4th Question:
Why have you not posted this on the R-sig-mixed-models list, where it
clearly belongs, rather than here?
-- Bert
On Sat, Dec 29, 2012 at 12:47 PM, Diego Pujoni
<diegopujoni@gmail.com>wrote:
> Dear colleagues,
>
> I have a data from a repeated measures design that I'm analysing
through a
> mixed model. Nine independent sampling units (flasks with culture medium
> with algae) were randomly divided into 3 groups ("c",
"t1", "t2"). There is
> no need for inclusion of the random effect of the intercept, because the
> nine sample units are homogeneous among each other (samples taken from the
> same culture). The algal concentration was measured every two days for 10
> days. The goal is to test differences between treatments.
>
> I estimated a model with only the intercept and the interaction Group * Day
> to test which is the best:
>
> library(nlme)
> library(lme4)
> Day = rep(c(0,2,4,6,8,10),each=9)
> Group =
rep(c("c","c","c","t1","t1","t1","t2","t2","t2"),6)
> Individual = rep(1:9,6)
> X = c(0.71,0.72,0.71,0.72,0.72,0.72,0.70,0.69,0.70,0.72,0.72,
> 0.71,0.72,0.72,0.71,0.71,0.70,0.71,0.73,0.73,0.69,0.74,
> 0.69,0.73,0.67,0.71,0.69,0.71,0.71,0.72,0.70,0.71,0.70,
> 0.52,0.64,0.60,0.70,0.73,0.73,0.67,0.66,0.71,0.47,0.56,
> 0.54,0.65,0.73,0.73,0.67,0.71,0.58,0.44,0.52,0.58)
>
> xyplot(X~Day, groups=Group)
>
> LME = lme(X ~ 1, random = ~Day|Individual)
> Erro em lme.formula(X ~ 1, random = ~Day | Individual) :
> nlminb problem, convergence error code = 1
> message = iteration limit reached without convergence (10)
>
> LME1 = lme(X ~ Group*Day, random = ~Day|Individual)
> Erro em lme.formula(X ~ Group * Day, random = ~Day | Individual) :
> nlminb problem, convergence error code = 1
> message = iteration limit reached without convergence (10)
>
> LMER = lmer(X ~ 1 + (Day|Individual))
> LMER1 = lmer(X ~ Group*Day + (Day|Individual))
>
> AIC(LMER)
> [1] -179.0302
>
> AIC(LMER1)
> [1] -151.1938
>
> anova(LMER,LMER1)
> Data:
> Models:
> LMER: X ~ 1 + (Day | Individual)
> LMER1: X ~ Group * Day + (Day | Individual)
> Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
> LMER 5 -187.2 -177.26 98.602
> LMER1 10 -203.2 -183.31 111.600 25.996 5 8.939e-05 ***
>
> xyplot(fitted(LMER)~Day, groups=Group)
> xyplot(fitted(LMER1)~Day, groups=Group)
>
> 1st question: Why function lme4:lmer converge, but the nlme:lme
doesn't?
> The first is better than the second?
>
> 2nd question: Why does the "anova" give distinct values of AIC
for the two
> models. If we look the AIC value of each model, the best model is LMER, but
> the "anova" says that LMER1 is the best.
>
> 3rd question: Why the fitted values of the model with only the intercept
> (LMER) vary over time?
>
> Thank you very much for any help
>
>
> Diego PJ
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
--
Bert Gunter
Genentech Nonclinical Biostatistics
Internal Contact Info:
Phone: 467-7374
Website:
http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm
[[alternative HTML version deleted]]