[R] Resid() and estimable() functions with lmer

[Adriana Puentes]

Hi all,

Two questions:

1. Is there a way to evaluate models from lmer() with a poisson  
distribution? I get the following error message:

library(lme4)
lmer(tot.fruit~infl.treat+def.treat+(1|initial.size),family=poisson)->model
plot(fitted(model),resid(model))
Error: 'resid' is not implemented yet

Are there any other options?

2. Why doesn't the function estimable() in gmodels work with lmer()  
using a poisson distribution? I know that my coefficients are right  
because it works fine if I run the model without poisson distribution.

#CODE (with latest R version and package updates)#
#Without Poisson distribution#

library(lme4)
library(gmodels)
lmer(tot.fruit~infl.treat+def.treat+(1|initial.size))->model1
summary(model1)

Linear mixed-effects model fit by REML
Formula: tot.fruit ~ infl.treat + def.treat + (1 | initial.size)
   AIC  BIC logLik MLdeviance REMLdeviance
  2449 2471  -1218       2447         2437
Random effects:
  Groups       Name        Variance Std.Dev.
  initial.size (Intercept) 71.585   8.4608
  Residual                 49.856   7.0608
number of obs: 323, groups: initial.size, 292

Fixed effects:
             Estimate Std. Error t value
(Intercept)  12.8846     1.3028   9.890
infl.treat1  -0.4738     1.1819  -0.401
def.treat2   -3.5522     1.6022  -2.217
def.treat3   -2.1757     1.6461  -1.322
def.treat4   -2.1613     1.7003  -1.271

Correlation of Fixed Effects:
             (Intr) infl.1 df.tr2 df.tr3
infl.treat1 -0.413
def.treat2  -0.616  0.013
def.treat3  -0.641  0.002  0.493
def.treat4  -0.638  0.028  0.469  0.524

#Coefficients#

mean.infl.treat0=c(1,0,1/4,1/4,1/4)
mean.infl.treat1=c(1,1,1/4,1/4,1/4)
mean.def.treat1=c(1,1/2,0,0,0)
mean.def.treat2=c(1,1/2,1,0,0)
mean.def.treat3=c(1,1/2,0,1,0)
mean.def.treat4=c(1,1/2,0,0,1)
means=rbind(mean.infl.treat0,mean.infl.treat1,mean.def.treat1,mean.def.treat2,mean.def.treat3,mean.def.treat4)
estimable(model1,means,sim.lmer=TRUE,n.sim=1000)

                       Estimate Std. Error p value
(1 0 0.25 0.25 0.25) 10.897825  0.8356260       0
(1 1 0.25 0.25 0.25) 10.421633  0.9364839       0
(1 0.5 0 0 0)        12.577408  1.2055979       0
(1 0.5 1 0 0)         9.097063  1.1769760       0
(1 0.5 0 1 0)        10.523485  1.2238635       0
(1 0.5 0 0 1)        10.440959  1.2031459       0

#With Poisson distribution#

lmer(tot.fruit~infl.treat+def.treat+(1|initial.size),family=poisson)->model2
summary(model2)

Generalized linear mixed model fit using Laplace
Formula: tot.fruit ~ infl.treat + def.treat + (1 | initial.size)
  Family: poisson(log link)
   AIC  BIC logLik deviance
  1073 1096 -530.5     1061
Random effects:
  Groups       Name        Variance Std.Dev.
  initial.size (Intercept) 0.84201  0.91761
number of obs: 323, groups: initial.size, 292

Estimated scale (compare to  1 )  1.130529

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)  2.08865    0.10478  19.934  < 2e-16 ***
infl.treat1  0.10615    0.09412   1.128  0.25941
def.treat2  -0.37354    0.11398  -3.277  0.00105 **
def.treat3  -0.18566    0.12679  -1.464  0.14310
def.treat4  -0.10624    0.13451  -0.790  0.42962
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1

Correlation of Fixed Effects:
             (Intr) infl.1 df.tr2 df.tr3
infl.treat1 -0.418
def.treat2  -0.554  0.016
def.treat3  -0.599 -0.001  0.467
def.treat4  -0.623  0.055  0.460  0.548

estimable(model2,lsmeans,sim.lmer=TRUE,n.sim=1000)
Error in FUN(newX[, i], ...) :
   `param' has no names and does not match number of coefficients of  
model. Unable to construct coefficient vector

#End of CODE#

Any ideas as to why this is? It does say in the help of ?estimable  
that it should work for lmer objects.

Thanks in advance for any help,

Adriana Puentes

-- 
Adriana Puentes (MSc.)
PhD student
Plant Ecology
Evolutionary Biology Centre (EBC)
Uppsala University
Villavägen 14
75236, Uppsala,
Sweden

Webpage: http://www.vaxtbio.uu.se/resfold/puentes.htm
Work: (+46)18-471-2882
Mobile: (+46)(0)76-853-5861
Fax: (+46)18-553-419