[R-sig-ME] Resid() and estimable() functions with lmer

[Gustaf Granath]

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,

Gustaf (PhD student)


Dept of Plant Ecology
Evolutionary Biology Centre (EBC)
Uppsala University