Cómo tratar el ajuste singular en los modelos mixtos

Question

Digamos que tenemos un modelo

mod <- Y ~ X*Condition + (X*Condition|subject)

# Y = logit variable  
# X = continuous variable  
# Condition = values A and B, dummy coded; the design is repeated 
#             so all participants go through both Conditions  
# subject = random effects for different subjects 

summary(model)
Random effects:
 Groups  Name             Variance Std.Dev. Corr             
 subject (Intercept)      0.85052  0.9222                    
         X                0.08427  0.2903   -1.00            
         ConditionB       0.54367  0.7373   -0.37  0.37      
         X:ConditionB     0.14812  0.3849    0.26 -0.26 -0.56
Number of obs: 39401, groups:  subject, 219

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)       2.49686    0.06909   36.14  < 2e-16 ***
X                -1.03854    0.03812  -27.24  < 2e-16 ***
ConditionB       -0.19707    0.06382   -3.09  0.00202 ** 
X:ConditionB      0.22809    0.05356    4.26 2.06e-05 ***

Aquí observamos un ajuste singular, porque la correlación entre el intercepto y los efectos aleatorios x es -1. Ahora, según este útil enlace Una forma de tratar este modelo es eliminar los efectos aleatorios de orden superior (por ejemplo, X:ConditionB) y ver si eso supone una diferencia a la hora de probar la singularidad. La otra es utilizar el enfoque bayesiano, por ejemplo, el blme paquete para evitar la singularidad.

¿Cuál es el método preferido y por qué?

Lo pregunto porque utilizar la primera o la segunda conduce a resultados diferentes - en el primer caso, eliminaré el efecto aleatorio X:CondiciónB y no podré estimar la correlación entre los efectos aleatorios X y X:CondiciónB. Por otro lado, utilizando blme me permite mantener X:ConditionB y estimar la correlación dada. No veo ninguna razón por la que deba siquiera utilizar las estimaciones no bayesianas y eliminar los efectos aleatorios cuando se producen ajustes singulares cuando puedo estimar todo con el enfoque bayesiano.

¿Puede alguien explicarme las ventajas y los problemas de utilizar uno u otro método para tratar los ataques singulares?

Gracias.

Answer 1

Cuando se obtiene un ajuste singular, esto suele indicar que el modelo está sobreajustado, es decir, que la estructura de efectos aleatorios es demasiado compleja para ser soportada por los datos, lo que naturalmente lleva a aconsejar la eliminación de la parte más compleja de la estructura de efectos aleatorios (normalmente las pendientes aleatorias). La ventaja de este enfoque es que conduce a un modelo más parsimonioso que no está sobreajustado.

Sin embargo, antes de hacer nada, ¿tiene una buena razón para querer X , Condition y su interacción, todo para variar según el tema en primer lugar? ¿Sugiere esto la teoría de cómo se generan los datos?

Si desea ajustar el modelo con la estructura de efectos aleatorios máxima, y lme4 obtiene un ajuste singular, entonces el ajuste del mismo modelo en un marco bayesiano podría muy bien informarle por qué lme4 tuvo problemas, inspeccionando los gráficos de rastreo y lo bien que convergen las distintas estimaciones de los parámetros. La ventaja de adoptar el enfoque bayesiano es que al hacerlo se puede descubrir un problema con el modelo original, es decir, la razón por la que la estructura de efectos aleatorios máximos no está respaldada por los datos) o podría descubrir por qué lme4 es incapaz de ajustarse al modelo. Me he encontrado con situaciones en las que un modelo bayesiano no converge bien, a menos que se utilicen priores informativos, lo que puede estar bien o no.

En resumen, ambos enfoques tienen mérito.

Sin embargo, yo siempre partiría de un punto en el que el modelo inicial es parsimonioso y se basa en el conocimiento experto del dominio para determinar la estructura de efectos aleatorios más adecuada. Especificar las variables de agrupación es relativamente fácil, pero las pendientes aleatorias no suelen tienen a incluir. Inclúyalos sólo si tienen un sentido teórico sólido Y están respaldados por los datos.

Editar: Se menciona en los comentarios que hay razones teóricas sólidas para ajustar la estructura de efectos aleatorios máximos. Por lo tanto, una forma relativamente fácil de proceder con un modelo bayesiano equivalente es cambiar la llamada a glmer con stan_glmer de la rstanarm paquete - está diseñado para ser plug and play. Cuenta con valores previos por defecto, por lo que se puede obtener rápidamente un modelo ajustado. El paquete también tiene muchas herramientas para evaluar la convergencia. Si se comprueba que todos los parámetros convergen a valores plausibles, todo va bien. Sin embargo, puede haber una serie de problemas: por ejemplo, una varianza estimada en o por debajo de cero, o una estimación que continúa a la deriva. El sitio web mc-stan.org contiene una gran cantidad de información y un foro de usuarios.

Answer 2

¡Este es un hilo muy interesante, con respuestas y comentarios interesantes! Ya que esto no se ha mencionado todavía, quería señalar que tenemos muy pocos datos para cada sujeto (según tengo entendido). En efecto, cada sujeto sólo tiene dos valores para cada una de las variables de respuesta Y, la variable categórica Condición y la variable continua X. En concreto, sabemos que los dos valores de Condición son A y B.

Si siguiéramos la modelización de regresión en dos etapas en lugar de la modelización de efectos mixtos, ni siquiera podríamos ajustar un modelo de regresión lineal a los datos de un sujeto específico, como se ilustra en el ejemplo de juguete siguiente para uno de los sujetos:

y <- c(4, 7)
condition <- c("A", "B")
condition <- factor(condition)
x <- c(0.2, 0.4)

m <- lm(y ~ condition*x)
summary(m)

El resultado de este modelo específico sería:

Call:
lm(formula = y ~ condition * x)

Residuals:
ALL 2 residuals are 0: no residual degrees of freedom!

Coefficients: (2 not defined because of singularities)
         Estimate Std. Error t value Pr(>|t|)
(Intercept)         4         NA      NA       NA
conditionB          3         NA      NA       NA
x                  NA         NA      NA       NA
conditionB:x       NA         NA      NA       NA

Residual standard error: NaN on 0 degrees of freedom
Multiple R-squared:      1,     Adjusted R-squared:    NaN 
F-statistic:   NaN on 1 and 0 DF,  p-value: NA

Obsérvese que el ajuste del modelo adolece de singularidades, ya que estamos tratando de estimar 4 coeficientes de regresión más la desviación estándar del error utilizando sólo 2 observaciones.

Las singularidades persistirían incluso si observáramos a este sujeto dos veces -en lugar de una- bajo cada condición. Sin embargo, si observáramos al sujeto 3 veces bajo cada condición, nos libraríamos de las singularidades:

y <- c(4, 7, 3, 5, 1, 2)
condition <- c("A", "B", "A","B","A","B")
condition <- factor(condition)
x <- c(0.2, 0.4, 0.1, 0.3, 0.3, 0.5)

m2 <- lm(y ~ condition*x)
summary(m2)

Aquí está la salida R correspondiente a este segundo ejemplo, del que han desaparecido las singularidades:

>     summary(m2)

Call:
lm(formula = y ~ condition * x)

Residuals:
    1       2       3       4       5       6 
1.3333  2.3333 -0.6667 -1.1667 -0.6667 -1.1667 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)     4.667      3.555   1.313    0.320
conditionB      6.000      7.601   0.789    0.513
x             -10.000     16.457  -0.608    0.605
conditionB:x   -5.000     23.274  -0.215    0.850

Residual standard error: 2.327 on 2 degrees of freedom
Multiple R-squared:  0.5357,    Adjusted R-squared:  -0.1607 
F-statistic: 0.7692 on 3 and 2 DF,  p-value: 0.6079

Por supuesto, el modelo de efectos mixtos no se ajusta a modelos de regresión lineal no relacionados y separados para cada sujeto, sino que se ajusta a modelos "relacionados" cuyos interceptos y/o pendientes se desvían aleatoriamente en torno a un intercepto y/o pendiente típicos, de forma que las desviaciones aleatorias del intercepto y/o la pendiente típicos siguen una distribución normal con media cero y alguna desviación estándar desconocida.

Aun así, mi intuición sugiere que el modelo de efectos mixtos tiene dificultades con la pequeña cantidad de observaciones -sólo 2- disponibles para cada sujeto. Cuanto más se cargue el modelo con pendientes aleatorias, más dificultades tendrá probablemente. Sospecho que, si cada sujeto aportara 6 observaciones en lugar de 2 (es decir, 3 por condición), ya no tendría problemas para acomodar todas las pendientes aleatorias.

Me parece que esto podría ser (?) un caso en el que el diseño actual del estudio no apoya las ambiciones de modelización compleja - para apoyar esas ambiciones, se necesitarían más observaciones en cada condición para cada sujeto (o al menos para algunos de los sujetos?). Esto es sólo mi intuición, así que espero que otros puedan añadir sus ideas a mis observaciones anteriores. Gracias de antemano.