[R-sig-ME] Resid() and estimable() functions with lmer
Douglas Bates
bates at stat.wisc.edu
Wed Nov 14 14:17:45 CET 2007
On Nov 14, 2007 5:58 AM, Gustaf Granath <Gustaf.Granath at ebc.uu.se> wrote:
> 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
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.
> 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.
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.
> #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
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
More information about the R-sig-mixed-models
mailing list