Realicé un experimento donde crié diferentes familias provenientes de dos poblaciones de origen diferentes. A cada familia se le asignó uno de dos tratamientos. Después del experimento medí varios rasgos en cada individuo. Para probar un efecto del tratamiento o la fuente, así como su interacción, utilicé un modelo lineal de efectos mixtos con la familia como factor aleatorio, es decir
lme(fixed=Trait~Treatment*Source,random=~1|Family,method="ML")
hasta ahora todo bien, ahora tengo que calcular los componentes de la varianza relativa, es decir, el porcentaje de variación que se explica por el tratamiento o la fuente, así como la interacción.
Sin un efecto aleatorio, podría usar fácilmente las sumas de cuadrados (SS) para calcular la varianza explicada por cada factor. Pero para un modelo mixto (con estimación de ML), no hay SS, por lo tanto, pensé que también podría usar Tratamiento y Fuente como efectos aleatorios para estimar la varianza, es decir
lme(fixed=Trait~1,random=~(Treatment*Source)|Family, method="REML")
Sin embargo, en algunos casos, lme no converge, por lo tanto, utilicé lmer del paquete lme4:
lmer(Trait~1+(Treatment*Source|Family),data=DATA)
Donde extraigo las variaciones del modelo usando la función de resumen:
model<-lmer(Trait~1+(Treatment*Source|Family),data=regrexpdat)
results<-VarCorr(model)
variances<-results[,3]
Obtengo los mismos valores que con la función VarCorr. Entonces uso estos valores para calcular el porcentaje real de variación tomando la suma como la variación total.
Donde estoy luchando es con la interpretación de los resultados del modelo inicial de lme (con tratamiento y fuente como efectos fijos) y el modelo aleatorio para estimar los componentes de varianza (con tratamiento y fuente como efecto aleatorio). Encuentro en la mayoría de los casos que el porcentaje de varianza explicado por cada factor no corresponde a la importancia del efecto fijo.
Por ejemplo, para el rasgo HD, la película inicial sugiere una tendencia a la interacción, así como una importancia para el tratamiento. Usando un procedimiento hacia atrás, encuentro que el tratamiento tiene una tendencia cercana a significativa. Sin embargo, al estimar los componentes de la varianza, encuentro que Source tiene la varianza más alta, representando hasta el 26.7% de la varianza total.
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
Y el 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
Por lo tanto, mi pregunta es, ¿es correcto lo que estoy haciendo? O debería usar otra forma de estimar la cantidad de varianza explicada por cada factor (es decir, Tratamiento, Fuente y su interacción). Por ejemplo, ¿serían los tamaños del efecto una forma más apropiada de hacerlo?
Respuestas:
Una forma común de determinar la contribución relativa de cada factor a un modelo es eliminar el factor y comparar la probabilidad relativa con algo como una prueba de ji cuadrado:
Como la forma en que se calculan las probabilidades entre funciones puede ser ligeramente diferente, típicamente solo las compararé entre el mismo método.
fuente
1-pchisq(logLik(model1) - logLik(model2), 1)
?