Skip to content

Resid() and estimable() functions with lmer

2 messages · Gustaf Granath, Douglas Bates

Original
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
Douglas Bates
#
On Nov 14, 2007 5:58 AM, Gustaf Granath <Gustaf.Granath at ebc.uu.se> wrote:
The model can be fit, it is the resid function that is not yet
implemented.  The reason it is not yet implemented is because there
are many different kinds of residuals for generalized linear mixed
models.  It requires considerably more design and coding than does
obtaining the residuals for a linear mixed model, and even those are
difficult to define for the general model (Do you incorporate the
random effects or not?  If you do, what do you use in the case of
multiple grouping factors?)

Right now I am thinking about other issues in the form of the model so
I can incorporate both generalized linear mixed models and nonlinear
mixed models.  Residuals for generalized linear mixed models are
further down the priority list.
I don't know what the estimable function expects the structure of an
lmer or glmer model to be but it would not surprise me if it did not
coincide with the current version.    The structure is continuously
evolving, not because I want to make things difficult for others but
because my understanding of the model and the computational methods
evolves.  On something as complicated as a generalized nonlinear mixed
model with crossed grouping factors it is not surprising that one
doesn't get the design down on the first try.  At least this one
doesn't.  I hope the structure of the classes and the algorithms will
stabilize soon (say, by the end of the year).

Even then it will be difficult to implement something like estimable
without doing finite difference calculations for generalized linear
mixed models.  A matrix that determines the marginal variance
covariance for the fixed-effects estimates (well, not quite but I
don't want to get into the technicalities) is calculated for linear
mixed models but not for generalized linear mixed models.