Comment estimer les composantes de la variance avec lmer pour les modèles à effets aléatoires et les comparer avec les résultats lme

14

J'ai réalisé une expérience où j'ai élevé différentes familles provenant de deux populations sources différentes. Chaque famille a reçu l'un des deux traitements. Après l'expérience, j'ai mesuré plusieurs traits sur chaque individu. Pour tester un effet du traitement ou de la source ainsi que leur interaction, j'ai utilisé un modèle à effet mixte linéaire avec la famille comme facteur aléatoire, c'est-à-dire

lme(fixed=Trait~Treatment*Source,random=~1|Family,method="ML")

jusqu'ici tout va bien, maintenant je dois calculer les composantes de la variance relative, c'est-à-dire le pourcentage de variation qui est expliqué par le traitement ou la source ainsi que l'interaction.

Sans effet aléatoire, je pourrais facilement utiliser les sommes des carrés (SS) pour calculer la variance expliquée par chaque facteur. Mais pour un modèle mixte (avec estimation ML), il n'y a pas de SS, donc j'ai pensé que je pouvais aussi utiliser Traitement et Source comme effets aléatoires pour estimer la variance, c'est-à-dire

lme(fixed=Trait~1,random=~(Treatment*Source)|Family, method="REML")

Cependant, dans certains cas, lme ne converge pas, j'ai donc utilisé lmer du package lme4:

lmer(Trait~1+(Treatment*Source|Family),data=DATA)

Où j'extrais les variances du modèle en utilisant la fonction de résumé:

model<-lmer(Trait~1+(Treatment*Source|Family),data=regrexpdat)
results<-VarCorr(model)
variances<-results[,3]

J'obtiens les mêmes valeurs qu'avec la fonction VarCorr. J'utilise ensuite ces valeurs pour calculer le pourcentage réel de variation en prenant la somme comme variation totale.

Là où je me bats, c'est avec l'interprétation des résultats du modèle lme initial (avec traitement et source comme effets fixes) et le modèle aléatoire pour estimer les composantes de la variance (avec traitement et source comme effet aléatoire). Je trouve dans la plupart des cas que le pourcentage de variance expliqué par chaque facteur ne correspond pas à la signification de l'effet fixe.

Par exemple pour le trait HD, la lme ​​initiale suggère une tendance à l'interaction ainsi qu'une signification pour le traitement. En utilisant une procédure en arrière, je trouve que le traitement a une tendance proche de significative. Cependant, en estimant les composantes de la variance, je trouve que la source a la variance la plus élevée, représentant jusqu'à 26,7% de la variance totale.

La lme:

anova(lme(fixed=HD~as.factor(Treatment)*as.factor(Source),random=~1|as.factor(Family),method="ML",data=test),type="m")
                                      numDF denDF  F-value p-value
(Intercept)                                1   426 0.044523  0.8330
as.factor(Treatment)                       1   426 5.935189  0.0153
as.factor(Source)                          1    11 0.042662  0.8401
as.factor(Treatment):as.factor(Source)     1   426 3.754112  0.0533

Et le lmer:

summary(lmer(HD~1+(as.factor(Treatment)*as.factor(Source)|Family),data=regrexpdat))
Linear mixed model fit by REML 
Formula: HD ~ 1 + (as.factor(Treatment) * as.factor(Source) | Family) 
   Data: regrexpdat 
    AIC    BIC logLik deviance REMLdev
 -103.5 -54.43  63.75   -132.5  -127.5
Random effects:
 Groups   Name                                      Variance  Std.Dev. Corr                 
 Family   (Intercept)                               0.0113276 0.106431                      
          as.factor(Treatment)                      0.0063710 0.079819  0.405               
          as.factor(Source)                         0.0235294 0.153393 -0.134 -0.157        
          as.factor(Treatment)L:as.factor(Source)   0.0076353 0.087380 -0.578 -0.589 -0.585 
 Residual                                           0.0394610 0.198648                      
Number of obs: 441, groups: Family, 13

Fixed effects:
            Estimate Std. Error t value
(Intercept) -0.02740    0.03237  -0.846

Par conséquent, ma question est, est-ce correct ce que je fais? Ou devrais-je utiliser une autre façon d'estimer la quantité de variance expliquée par chaque facteur (c.-à-d. Traitement, source et leur interaction). Par exemple, les tailles d'effet seraient-elles une manière plus appropriée de procéder?

KL_STKBK
la source
Le facteur de traitement a 40 fois plus de degrés de liberté que le facteur source (pseudo-réplication?). Cela fait sans aucun doute baisser la valeur P du traitement.

Réponses:

1

Une façon courante de déterminer la contribution relative de chaque facteur à un modèle consiste à supprimer le facteur et à comparer la probabilité relative avec quelque chose comme un test du chi carré:

pchisq(logLik(model1) - logLik(model2), 1)

Comme la façon dont les probabilités sont calculées entre les fonctions peut être légèrement différente, je ne les compare généralement qu'entre la même méthode.

Bill Denney
la source
1
ne devrait-il pas en être ainsi 1-pchisq(logLik(model1) - logLik(model2), 1)?
user81411