¿Por qué obtengo una variación cero de un efecto aleatorio en mi modelo mixto, a pesar de alguna variación en los datos?

Hemos ejecutado una regresión logística de efectos mixtos utilizando la siguiente sintaxis;

# fit model
fm0 <- glmer(GoalEncoding ~ 1 + Group + (1|Subject) + (1|Item), exp0,
             family = binomial(link="logit"))
# model output
summary(fm0)

Sujeto y elemento son los efectos aleatorios. Estamos obteniendo un resultado impar, que es el coeficiente y la desviación estándar para el término sujeto son ambos cero;

Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial  ( logit )
Formula: GoalEncoding ~ 1 + Group + (1 | Subject) + (1 | Item)
Data: exp0

AIC      BIC      logLik deviance df.resid 
449.8    465.3   -220.9    441.8      356 

Scaled residuals: 
Min     1Q Median     3Q    Max 
-2.115 -0.785 -0.376  0.805  2.663 

Random effects:
Groups  Name        Variance Std.Dev.
Subject (Intercept) 0.000    0.000   
Item    (Intercept) 0.801    0.895   
Number of obs: 360, groups:  Subject, 30; Item, 12

Fixed effects:
                Estimate Std. Error z value Pr(>|z|)    
 (Intercept)     -0.0275     0.2843    -0.1     0.92    
 GroupGeMo.EnMo   1.2060     0.2411     5.0  5.7e-07 ***
 ---
 Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Correlation of Fixed Effects:
             (Intr)
 GroupGM.EnM -0.002

Esto no debería estar sucediendo porque obviamente hay variación entre los sujetos. Cuando ejecutamos el mismo análisis en stata

xtmelogit goal group_num || _all:R.subject || _all:R.item

Note: factor variables specified; option laplace assumed

Refining starting values: 

Iteration 0:   log likelihood = -260.60631  
Iteration 1:   log likelihood = -252.13724  
Iteration 2:   log likelihood = -249.87663  

Performing gradient-based optimization: 

Iteration 0:   log likelihood = -249.87663  
Iteration 1:   log likelihood = -246.38421  
Iteration 2:   log likelihood =  -245.2231  
Iteration 3:   log likelihood = -240.28537  
Iteration 4:   log likelihood = -238.67047  
Iteration 5:   log likelihood = -238.65943  
Iteration 6:   log likelihood = -238.65942  

Mixed-effects logistic regression               Number of obs      =       450
Group variable: _all                            Number of groups   =         1

                                                Obs per group: min =       450
                                                               avg =     450.0
                                                               max =       450

Integration points =   1                        Wald chi2(1)       =     22.62
Log likelihood = -238.65942                     Prob > chi2        =    0.0000

------------------------------------------------------------------------------
        goal |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
   group_num |   1.186594    .249484     4.76   0.000     .6976147    1.675574
       _cons |  -3.419815   .8008212    -4.27   0.000    -4.989396   -1.850234
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
_all: Identity               |
               sd(R.subject) |   7.18e-07   .3783434             0           .
-----------------------------+------------------------------------------------
_all: Identity               |
                 sd(R.trial) |   2.462568   .6226966      1.500201    4.042286
------------------------------------------------------------------------------
LR test vs. logistic regression:     chi2(2) =   126.75   Prob > chi2 = 0.0000

Note: LR test is conservative and provided only for reference.
Note: log-likelihood calculations are based on the Laplacian approximation.

los resultados son los esperados con un coeficiente diferente de cero / se para el término sujeto.

Originalmente pensamos que esto podría tener algo que ver con la codificación del término Asunto, pero cambiar esto de una cadena a un entero no hizo ninguna diferencia.

Obviamente, el análisis no funciona correctamente, pero no podemos precisar la fuente de las dificultades. (Nota: alguien más en este foro ha tenido un problema similar, pero este hilo sigue sin respuesta ).

Esto se analiza con cierta extensión en https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html (busque "modelos singulares"); es común, especialmente cuando hay un pequeño número de grupos (aunque 30 no es particularmente pequeño en este contexto).

Una diferencia entre lme4muchos otros paquetes es que muchos paquetes, incluido lme4el predecesor nlme, manejan el hecho de que las estimaciones de varianza no deben ser negativas ajustando la varianza en la escala logarítmica: eso significa que las estimaciones de varianza no pueden ser exactamente cero, simplemente muy muy pequeña. lme4, por el contrario, utiliza la optimización restringida, por lo que puede devolver valores que son exactamente cero (consulte http://arxiv.org/abs/1406.5823 p. 24 para obtener más información). http://rpubs.com/bbolker/6226 da un ejemplo.

En particular, si observa de cerca los resultados de la varianza entre sujetos de Stata, tiene una estimación de 7.18e-07 (en relación con una intersección de -3.4) con una desviación estándar de Wald de .3783434 (¡esencialmente inútil en este caso!) Y un IC del 95% listado como "0"; esto es técnicamente "distinto de cero", pero es tan cercano a cero como informará el programa ...

Es bien conocido y teóricamente demostrable (por ejemplo, Stram y Lee Biometrics 1994) que la distribución nula para los componentes de la varianza es una mezcla de una masa puntual ('pico') en cero y una distribución chi-cuadrado alejada de cero. Como era de esperar (pero no sé si está comprobado / bien conocido), la distribución de muestreo de las estimaciones del componente de varianza a menudo tiene un pico en cero, incluso cuando el valor verdadero no es cero; consulte, por ejemplo, http://rpubs.com/ bbolker / 4187 para un ejemplo, o el último ejemplo en la ?bootMerpágina:

library(lme4)
library(boot)
## Check stored values from a longer (1000-replicate) run:
load(system.file("testdata","boo01L.RData",package="lme4"))
plot(boo01L,index=3) 

No creo que haya un problema. La lección del resultado del modelo es que, aunque existe una variación "obviamente" en el rendimiento del sujeto, el alcance de esta variación del sujeto puede explicarse total o virtualmente solo por el término de varianza residual solo. No hay suficiente variación adicional a nivel de sujeto para justificar la adición de un efecto aleatorio adicional a nivel de sujeto para explicar toda la variación observada.

Piénsalo de esta manera. Imagine que estamos simulando datos experimentales bajo este mismo paradigma. Establecimos los parámetros para que haya una variación residual ensayo por ensayo, pero 0 variación a nivel de sujeto (es decir, todos los sujetos tienen la misma "media verdadera" más error). Ahora, cada vez que simulamos datos de este conjunto de parámetros, por supuesto encontraremos que los sujetos no tienen exactamente el mismo rendimiento. Algunos terminan con puntajes bajos, algunos con puntajes altos. Pero todo esto se debe solo a la variación residual del nivel de prueba. "Sabemos" (en virtud de haber determinado los parámetros de simulación) que en realidad no hay ninguna variación a nivel de sujeto.