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
