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?
Réponses:
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é:
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.
la source
1-pchisq(logLik(model1) - logLik(model2), 1)
?