Cómo estimar componentes de varianza con lmer para modelos con efectos aleatorios y compararlos con los resultados de lme

14

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?

KL_STKBK
fuente
El factor de tratamiento tiene 40 veces más grados de libertad que el factor fuente (¿seudoreplicación?). Esto indudablemente está disminuyendo el valor P del tratamiento.

Respuestas:

1

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:

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

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.

Bill Denney
fuente
1
no debería ser 1-pchisq(logLik(model1) - logLik(model2), 1)?
user81411