¿Puede incluirse una variable en un modelo mixto como efecto fijo y como efecto aleatorio al mismo tiempo?

Question

Estoy construyendo un modelo mixto, donde los mismos 4 grupos son probados repetidamente 5 semanas en la medida M. Quiero probar el efecto del grupo en la medida M. ¿Puedo construir un modelo como el siguiente (digamos con R):

   lmer(M ~ Time + Group + (1|Group), data = mydata)

donde group ¿es un efecto fijo y aleatorio al mismo tiempo? Imponer dos efectos de un grupo parece conceptualmente contradictorio. Pero, ¿y si quiero probar el efecto de grupo?

Answer 1

Primero para asegurarse de que Group se define como un factor. Su modelo actual es $$M_{ij}=\alpha_i+\beta_1\mathrm{Time}_{ij}+u_i+e_{ij},$$ donde $i$ denota Group et $j$ los puntos de tiempo. Si ejecuta su modelo y prueba con ranef() se dará cuenta de que es difícil distinguir $\alpha_i$ et $u_i$ en la estimación, así $u_i$ casi igual a 0.

Dos posibles modelos alternativos son:

1. Modelo de efectos aleatorios, $$M_{ij}=\beta_0+\beta_1\mathrm{Time}_{ij}+u_i+e_{ij},$$ donde $\beta_0$ es la intercepción media, $u_i$ es la desviación individual aleatoria del intercepto medio de cada grupo y se supone que sigue una distribución.
2. Modelo de efectos fijos (en el contexto de la econometría), $$M_{ij}=\alpha_i+\beta_1\mathrm{Time}_{ij}+e_{ij},$$ donde $\alpha_i$ es el intercepto individual fijo.

---

Actualizaciones:

Utilizo el sleepstudy conjunto de datos como ilustración. Si la variable de agrupación (un factor) se incluye como covariable (Modelo fm2 más adelante), tanto los efectos aleatorios como su varianza tienden a ser cero. La explicación intuitiva es que, $\alpha_i$ et $u_i$ modelan básicamente la misma cantidad (el intercepto específico del grupo), aunque uno se supone fijo y otro aleatorio. La mayor parte de la variabilidad es absorbida primero por los interceptos fijos ( $\alpha_i$ ), por lo que los interceptos aleatorios $u_i$ tienden a ser todos cero. El código y los resultados se indican a continuación.

> fm1 <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy, REML = F)
> re1 = as.matrix(ranef(fm1)$Subject)
> summary(as.vector(re1))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-77.570  -7.460   5.701   0.000  15.850  71.920 
> VarCorr(fm1)
 Groups   Name        Std.Dev.
 Subject  (Intercept) 36.012  
 Residual             30.895  
> sd(re1)
[1] 35.76336

> fm2 <- lmer(Reaction ~ Days + Subject + (1 | Subject), sleepstudy, REML = F)
> re2 = as.matrix(ranef(fm2)$Subject)
> summary(as.vector(re2))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0       0       0       0       0       0 
> VarCorr(fm2)
 Groups   Name        Std.Dev.
 Subject  (Intercept)  0.00   
 Residual             29.31   
> sd(re2)
[1] 0

---

Actualizaciones 2:

Anteriormente utilicé ML en lmer porque encontré que la estimación de la varianza REML para el intercepto aleatorio parece estar fuera del valor real en este caso extremo. Sé que puede no tener sentido incluir tanto los interceptos fijos como los aleatorios específicos del grupo en un solo modelo, pero puede ser un ejemplo interesante. Tenga en cuenta que la única diferencia entre este modelo y el modelo de intercepción aleatoria es una matriz de diseño más grande para los efectos fijos.

Le site lmer El ejemplo siguiente utiliza REML, y el modelo es el mismo que el modelo fm2 arriba. Los efectos aleatorios estimados son todos cercanos a cero, y su desviación estándar es muy cercana a cero. Sé que la desviación estándar de los efectos aleatorios estimados no es una estimación perfecta de la varianza de los efectos aleatorios, pero los dos deberían corresponder de alguna manera. Pero la estimación de la varianza REML para los efectos aleatorios es de 33, que está muy lejos de cero.

> fm3 <- lmer(Reaction ~ Days + Subject + (1 | Subject), sleepstudy, REML = T, verbose = T)
> re3 = as.matrix(ranef(fm3)$Subject)
> summary(as.vector(re3))
      Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
-3.998e-12 -9.220e-13 -6.549e-13 -7.717e-13 -4.567e-13  6.893e-13 
> VarCorr(fm3)
 Groups   Name        Std.Dev.
 Subject  (Intercept) 33.050  
 Residual             30.991  
> sd(re3)
[1] 9.391908e-13

También he probado en Stata y se vuelve más interesante. El mixed utiliza un algoritmo EM pero no puede converger y por lo tanto da una estimación muy grande. A mi entender, REML y ML no deberían diferir tanto en este caso. Puede haber algunos problemas numéricos. Dado que las estimaciones dependen de las iteraciones, pensaré más en esto cuando tenga más tiempo.

. mixed reaction days i.subject || subject:, reml

Performing EM optimization: 

Performing gradient-based optimization: 

could not calculate numerical derivatives -- discontinuous region with missing values encountered
could not calculate numerical derivatives -- discontinuous region with missing values encountered

Computing standard errors:
standard-error calculation failed

Mixed-effects REML regression                   Number of obs     =        180
Group variable: subject                         Number of groups  =         18

                                                Obs per group:
                                                              min =         10
                                                              avg =       10.0
                                                              max =         10

                                                Wald chi2(18)     =     169.64
Log restricted-likelihood = -805.65036          Prob > chi2       =     0.0000

...

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
subject: Identity            |
                  var(_cons) |   105231.5          .             .           .
-----------------------------+------------------------------------------------
               var(Residual) |   960.4566          .             .           .
------------------------------------------------------------------------------
LR test vs. linear model: chibar2(01) = 2.3e-13       Prob >= chibar2 = 1.0000

Warning: convergence not achieved; estimates are based on iterated EM