Modèles mixtes

Description

On s'intéresse à la croissance des pommes sur des donnes relevées sur difftérents pommiers. Les variables sont:

La question est de savoir en quoi la croissance des pommes diffère d'un arbre à l'autre.

Les programmes SAS correspondants se trouvent à la page http://www.growthmodel.org/mixedmodels/mixedmodels.htm.
Calcul du modèle
Lire les données.
Ne conserver que celle pour lesquel Diam>0:

ap=ap0[ap0$Diam>0,]
Une figure par arbre
attach(ap);
library(nlme);
plot(groupedData(Diam~Time | Tree),xlab="Time (Week/2)", ylab="Apple Diameter (Inches)");
Un modèle linéaire par pomme. Attention: la pomme 14 n'a qu'une mesure.
Il s'agit de comparer les différents coefficients beta0 et beta1 obtenus. La commande by permet d'appliquer la même fonction a des sous-ensembles des données (un peu comme tapply):

bet0=by(ap, factor(AppleID), function(x) lm(Diam~Time, data=x)$coef[1]);
bet1=by(ap, factor(AppleID), function(x) lm(Diam~Time, data=x)$coef[2]);

Calcul des moyennes et variances correspondantes. En raison de la pomme 14, une paire (bet0,bet1) est NA (on peut le vérifier en tapant is.na(bet1)):

cat("beta0=",mean(bet0,na.rm=T),"(+-)",sd(bet0,na.rm=T));
cat("beta1=",mean(bet1,na.rm=T),"(+-)",sd(bet1,na.rm=T));

(ou bien mean(bet1[!is.na(bet1)]...).
Un modèle à effets aléatoires
On construit un modèle où mixte où l'effet de groupe est la pomme.

    y=b0+b1.temps+b0(pomme)+b1(pomme).temps+u  :

moda=lme(Diam~Time,random=~Time|AppleID)

Expliciter les valeurs de b0, b1, de l'écart-type de b0(.) et de b1(.), la corrélation entre b0(.) et b1(.).
Comparer avec les valeurs précédement calculées.
Comparer l'estimation de lme avec celle de lmer donné par

library(lme4)
modla=lmer(Diam~Time+(Time|AppleID))
Importance de l'arbre?
On peut aussi tracer les coefficients (toujours calculés sur chaque pomme) arbre par arbre. Ici, mdls0 est un vecteur qui associe a chaque pomme son arbre.

mdls0=by(ap, factor(AppleID), function(x) factor(x[1,]$Tree));
par(mfrow=c(1,2));
boxplot(as.vector(bet0)~as.vector(mdls0));
boxplot(as.vector(bet1)~as.vector(mdls0));
par(mfrow=c(1,1))

Une différence entre les arbres semble-t-elle apparaître?

On considère le modèle

   y=b0+b1.tps+b0(arb)+b1(arb).tps+b0(pom)+b1(pom).tps+u

Les effets «arbre» et «pomme» sont emboîtés.
On l'estime avec:

Tim=Time*200; # Echec de convergence sans cette anti-normalisation!
modat=lme(Diam~Tim,random=~Tim|Tree/AppleID)
   #...random=list(~Tim|Tree,~Tim|Apple)) #pareil
   #...random=~Tim|Tree/Apple) #pareil
   #...random=list(Tree=~Tim,Apple =~Tim)) #pareil

On va comparer modat et moda, ainsi que modat et modt (arbre seul), et aussi moda et mod=lm(Diam~Time).
Attention, on les réestimera auparavant au maximum de vraisemblance en utilisant l'option method="ML".
On commencera par comparer ces trois méthodes differentes pour le test de modat contre moda:

  1.  Test du rapport de vraisemblance:
  x2=2*logLik(modat)-2*logLik(moda)
  d=attr(logLik(modat),"df")-attr(logLik(moda),"df")
  pchisq(x2, df=d, lower.tail=F)

  2.  Utiliser "lmtest":
  library(lmtest)
  lrtest(moda,modat)

  3.  Utiliser "anova":
  anova(moda,modat)

Conclure
Avec lme4
Refaire les tests. Attention l'option pour le maximum de vraisemblance est ici REML=F

library(lme4)
modlat=lmer(Diam~Tim+(Tim|Tree)+(Tim|AppleID),REML=F)
modla=lmer(...)
modlt=lmer(...)

Obtient-on les mêmes résultats en estimation, en test?