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

Question

Realicé un experimento en el que crié a diferentes familias procedentes de dos poblaciones de origen distintas. A cada familia se le asignó uno de los dos tratamientos. Tras el experimento, medí varios rasgos en cada individuo. Para comprobar el efecto del tratamiento o de 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 va 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 utilizar fácilmente las sumas de los cuadrados (SS) para calcular la varianza explicada por cada factor. Pero para un modelo mixto (con estimación ML), no hay SS, por lo que pensé que podría utilizar el Tratamiento y la Fuente como efectos aleatorios también 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 que he utilizado lmer del paquete lme4:

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

Donde extraigo las varianzas del modelo utilizando 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. Utilizo entonces estos valores para calcular el porcentaje real de variación tomando la suma como la variación total.

El problema que tengo es la interpretación de los resultados del modelo inicial lme (con el tratamiento y la fuente como efectos fijos) y el modelo aleatorio para estimar los componentes de la varianza (con el tratamiento y la fuente como efecto aleatorio). En la mayoría de los casos encuentro que el porcentaje de varianza explicado por cada factor no se corresponde con la significación del efecto fijo.

Por ejemplo, para el rasgo HD, El lme inicial sugiere una tendencia para la interacción así como una significación para el Tratamiento. Utilizando un procedimiento hacia atrás, encuentro que el Tratamiento tiene una tendencia casi significativa. Sin embargo, al estimar los componentes de la varianza, encuentro que la Fuente tiene la mayor varianza, constituyendo 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

De ahí que mi pregunta sea, ¿es correcto lo que estoy haciendo? ¿O debería utilizar otra forma de estimar la cantidad de varianza explicada por cada factor (es decir, el tratamiento, la fuente y su interacción). Por ejemplo, ¿sería el tamaño del efecto una forma más apropiada?

Answer 1

Una forma habitual 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 chi-cuadrado:

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

Como la forma en que se calculan las probabilidades entre las funciones puede ser ligeramente diferente, normalmente sólo las compararé entre el mismo método.