[R-br] Seleção de modelo misto

Diego Pujoni diegopujoni em gmail.com
Sexta Dezembro 14 12:07:05 BRST 2012


Olá pessoal,

tenho um experimento de medidas repetidas que estou analisando através de
um modelo misto. Nove unidades amostrais independentes (frascos com meio de
cultura com algas) foram alocadas aleatoriamente em 3 grupos
("c","t1","t2"). Não há a necessidade de inclusão de efeito aleatório de
intercepto, pois as nove unidades amostrais são homogêneas entre si
(retiradas da mesma amostra mãe). A concentração de algas foi mensurada a
cada dois dias durante 10 dias. O objetivo é testar diferenças entre
os grupos (i.e. tratamentos).

Estimei um modelo somente com o intercepto e com a interação Group*Day,
para testar qual é o melhor

library(nlme)
library(lme4)
Day = rep(c(0,2,4,6,8,10),each=9)
Group = as.factor(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)

1a pergunta: Por que a função nlme:lme não converge, mas a lme4:lmer
converge? A segunda é melhor que a primeira?
2a pergunta: Por que a função "anova" dá valores distintos de AIC para os
dois modelos. Se olharmos o valor do AIC de cada um, o modelo LMER é
melhor que o LMER1, mas a "anova" diz o contrário.
3a pergunta: Por que os valores ajustados do modelo somente com o
intercepto (LMER) variam ao longo do tempo?

Muito obrigado pela ajuda

                                               Diego PJ
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20121214/2b83c160/attachment.html>


Mais detalhes sobre a lista de discussão R-br