¿Qué hacer con la correlación de efectos aleatorios que equivale a 1 o -1?

La ocurrencia no tan infrecuente cuando se trata de modelos mixtos máximos complejos (estimar todos los efectos aleatorios posibles para datos y modelo dados) es perfecta (+1 o -1) o una correlación casi perfecta entre algunos efectos aleatorios. Para el propósito de la discusión, observemos el siguiente modelo y resumen del modelo

Model: Y ~ X*Cond + (X*Cond|subj)

# 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  

Random effects:
 Groups  Name             Variance Std.Dev. Corr             
 subject (Intercept)      0.85052  0.9222                    
         X                0.08427  0.2903   -1.00            
         CondB            0.54367  0.7373   -0.37  0.37      
         X:CondB          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 ***
CondB            -0.19707    0.06382   -3.09  0.00202 ** 
X:CondB           0.22809    0.05356    4.26 2.06e-05 ***

La supuesta razón detrás de estas correlaciones perfectas es que hemos creado un modelo que es demasiado complejo para los datos que tenemos. El consejo común que se da en estas situaciones es (por ejemplo, Matuschek et al., 2017; papel ) para fijar los coeficientes sobreparamizados a 0, porque estos modelos degenerados tienden a disminuir la potencia. Si observamos un cambio marcado en los efectos fijos en un modelo reducido, deberíamos aceptarlo; Si no hay cambio, entonces no hay problema en aceptar el original.

Sin embargo, supongamos que no estamos interesados ​​solo en los efectos fijos controlados por RE (efectos aleatorios), sino también en la estructura de RE. En el caso dado, sería teóricamente correcto asumir que Intercepty la pendiente Xtiene una correlación negativa distinta de cero. Siguen varias preguntas:

1. ¿Qué hacer en tales situaciones? ¿Deberíamos informar la correlación perfecta y decir que nuestros datos no son "suficientemente buenos" para estimar la correlación "real"? ¿O deberíamos informar el modelo de correlación 0? ¿O deberíamos intentar establecer alguna otra correlación a 0 con la esperanza de que la "importante" ya no sea perfecta? No creo que haya respuestas 100% correctas aquí, me gustaría escuchar sus opiniones.

2. ¿Cómo escribir el código que fija la correlación de 2 efectos aleatorios específicos a 0, sin influir en las correlaciones entre otros parámetros?

Matrices de covarianza de efectos aleatorios singulares

Obtener una estimación de correlación de efecto aleatorio de +1 o -1 significa que el algoritmo de optimización alcanza "un límite": las correlaciones no pueden ser superiores a +1 o inferiores a -1. Incluso si no hay errores o advertencias de convergencia explícitos, esto potencialmente indica algunos problemas con la convergencia porque no esperamos que las correlaciones verdaderas se encuentren en el límite. Como dijiste, esto generalmente significa que no hay suficientes datos para estimar todos los parámetros de manera confiable. Matuschek y col. 2017 dice que en esta situación el poder puede verse comprometido.

Otra forma de alcanzar un límite es obtener una estimación de varianza de 0: ¿Por qué obtengo una varianza cero de un efecto aleatorio en mi modelo mixto, a pesar de alguna variación en los datos?

$4\times 4$$4\times 4$(preprint no publicado) recomienda utilizar el análisis de componentes principales (PCA) para verificar si la matriz de covarianza obtenida es singular. Si es así, sugieren tratar esta situación de la misma manera que las situaciones singulares anteriores.

¿Entonces lo que hay que hacer?

Si no hay suficientes datos para estimar todos los parámetros de un modelo de manera confiable, entonces deberíamos considerar simplificar el modelo. Tomando su modelo de ejemplo X*Cond + (X*Cond|subj), hay varias formas posibles de simplificarlo:

1. Elimine uno de los efectos aleatorios, generalmente la correlación de orden superior:

   X*Cond + (X+Cond|subj)

2. Deshágase de todos los parámetros de correlación:

   X*Cond + (X*Cond||subj)

   Actualización: como señala @Henrik, la ||sintaxis solo eliminará las correlaciones si todas las variables a la izquierda son numéricas. Si Condestán involucradas variables categóricas (como ), uno debería usar su conveniente afexpaquete (o soluciones manuales engorrosas). Vea su respuesta para más detalles.

3. Deshágase de algunos de los parámetros de correlación dividiendo el término en varios, por ejemplo:

   X*Cond + (X+Cond|subj) + (0+X:Cond|subj)

4. Restrinja la matriz de covarianza de alguna manera específica, por ejemplo, estableciendo una correlación específica (la que alcanza el límite) a cero, como sugiere. No hay una forma integrada de lme4lograr esto. Consulte la respuesta de @BenBolker en SO para obtener una demostración de cómo lograr esto a través de un pirateo inteligente.

Al contrario de lo que dijiste, no creo que Matuschek et al. 2017 recomienda específicamente # 4. La esencia de Matuschek et al. 2017 y Bates et al. 2015 parece ser que uno comienza con el modelo máximo a la Barr et al. 2013 y luego disminuye la complejidad hasta que la matriz de covarianza sea de rango completo. (Además, a menudo recomendarían reducir aún más la complejidad para aumentar la potencia). Actualización: en contraste, Barr et al. recomendar reducir SOLAMENTE la complejidad si el modelo no converge; están dispuestos a tolerar matrices de covarianza singulares. Ver la respuesta de @ Henrik.

Si uno está de acuerdo con Bates / Matuschek, entonces creo que está bien probar diferentes formas de disminuir la complejidad para encontrar el que hace el trabajo mientras hace "el menor daño". Mirando mi lista anterior, la matriz de covarianza original tiene 10 parámetros; # 1 tiene 6 parámetros, # 2 tiene 4 parámetros, # 3 tiene 7 parámetros. Es imposible decir qué modelo eliminará las correlaciones perfectas sin ajustarlas.

Pero, ¿qué pasa si estás interesado en este parámetro?

La discusión anterior trata la matriz de covarianza de efectos aleatorios como un parámetro molesto. Plantea una pregunta interesante sobre qué hacer si está específicamente interesado en un parámetro de correlación que debe "abandonar" para obtener una solución significativa de rango completo.

Tenga en cuenta que fijar el parámetro de correlación en cero no necesariamente generará BLUP ( ranef) que no están correlacionados; de hecho, es posible que ni siquiera se vean tan afectados (ver la respuesta de @ Placidia para una demostración ). Entonces, una opción sería mirar las correlaciones de los BLUP e informar eso.

Otra opción, quizás menos atractiva, sería utilizar tratar subjectcomo un efecto fijo Y~X*cond*subj, obtener las estimaciones para cada sujeto y calcular la correlación entre ellos. Esto es equivalente a ejecutar Y~X*condregresiones separadas para cada sujeto por separado y obtener las estimaciones de correlación de ellos.

Consulte también la sección sobre modelos singulares en las preguntas frecuentes sobre el modelo mixto de Ben Bolker:

$\theta$

Estoy de acuerdo con todo lo dicho en la respuesta de ameba que proporciona un gran resumen de la discusión actual sobre este tema. Intentaré agregar algunos puntos adicionales y, de lo contrario, me referiré al folleto de mi reciente curso de modelo mixto que también resume estos puntos.

La supresión de los parámetros de correlación (opciones 2 y 3 en la respuesta de ameba) ||solo funciona para covariables numéricas lmery no para factores. Esto se discute en detalle con el código de Reinhold Kliegl .

Sin embargo, mi afexpaquete proporciona la funcionalidad para suprimir la correlación también entre los factores si el argumento expand_re = TRUEen la llamada a mixed()(ver también la función lmer_alt()). Esencialmente lo hace implementando el enfoque discutido por Reinhold Kliegl (es decir, transfigurando los factores en covariables numéricas y especificando la estructura de efectos aleatorios en esos).

Un simple ejemplo:

library("afex")
data("Machines", package = "MEMSS") # same data as in Kliegl code

# with correlation:
summary(lmer(score ~ Machine + (Machine  | Worker), data=Machines))
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr       
#  Worker   (Intercept) 16.6405  4.0793              
#           MachineB    34.5467  5.8776    0.48      
#           MachineC    13.6150  3.6899   -0.37  0.30
#  Residual              0.9246  0.9616              
# Number of obs: 54, groups:  Worker, 6

## crazy results:
summary(lmer(score ~ Machine + (Machine  || Worker), data=Machines))
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr     
#  Worker   (Intercept)  0.2576  0.5076            
#  Worker.1 MachineA    16.3829  4.0476            
#           MachineB    74.1381  8.6103   0.80     
#           MachineC    19.0099  4.3600   0.62 0.77
#  Residual              0.9246  0.9616            
# Number of obs: 54, groups:  Worker, 6

## as expected:
summary(lmer_alt(score ~ Machine + (Machine  || Worker), data=Machines))
# Random effects:
#  Groups   Name         Variance Std.Dev.
#  Worker   (Intercept)  16.600   4.0743  
#  Worker.1 re1.MachineB 34.684   5.8894  
#  Worker.2 re1.MachineC 13.301   3.6471  
#  Residual               0.926   0.9623  
# Number of obs: 54, groups:  Worker, 6

Para aquellos que no lo saben afex, la funcionalidad principal para los modelos mixtos es proporcionar valores p para los efectos fijos, por ejemplo:

(m1 <- mixed(score ~ Machine + (Machine  || Worker), data=Machines, expand_re = TRUE))
# Mixed Model Anova Table (Type 3 tests, KR-method)
# 
# Model: score ~ Machine + (Machine || Worker)
# Data: Machines
#    Effect      df        F p.value
# 1 Machine 2, 5.98 20.96 **    .002
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1

summary(m1)  
# [...]
# Random effects:
#  Groups   Name         Variance Std.Dev.
#  Worker   (Intercept)  27.4947  5.2435  
#  Worker.1 re1.Machine1  6.6794  2.5845  
#  Worker.2 re1.Machine2 13.8015  3.7150  
#  Residual               0.9265  0.9626  
# Number of obs: 54, groups:  Worker, 6
# [...]

Dale Barr de Barr et al. (2013) el documento es más cauteloso al recomendar reducir la estructura de efectos aleatorios que la presentada en la respuesta de ameba. En un reciente intercambio de Twitter escribió:

• "reducir el modelo introduce un riesgo desconocido de anticonservatividad, y debe hacerse con precaución, si es que lo hace". y
• "Mi principal preocupación es que las personas entiendan los riesgos asociados con la reducción del modelo y que minimizar este riesgo requiera un enfoque más conservador de lo que comúnmente se adopta (por ejemplo, cada pendiente probada en .05)".

Por lo tanto, se recomienda precaución.

Como uno de los revisores, también puedo proporcionar una idea de por qué nosotros, Bates et al. (2015) el documento permaneció inédito. Yo y los otros dos revisores (que firmamos, pero permanecerán sin nombre aquí) tuvimos algunas críticas con el enfoque de PCA (parece no tener principios y no hay evidencia de que sea superior en términos de poder). Además, creo que los tres criticaron que el documento no se centró en la cuestión de cómo especificar la estructura de efectos aleatorios, sino que también intenta introducir GAMM. Por lo tanto, el artículo de Bates et al (2015) se transformó en Matuschek et al. (2017) documento que aborda el tema de la estructura de efectos aleatorios con simulaciones y Baayen et al. (2017) documento que presenta GAMM.

Mi revisión completa de Bates et al. borrador se puede encontrar aquí . IIRC, las otras revisiones tenían puntos principales similares.

Yo también he tenido este problema cuando uso la estimación de máxima verosimilitud: solo uso el algoritmo Goldstein IGLS implementado a través del software MLwiN y no LME4 en R. Sin embargo, en todos y cada uno de los casos el problema se resolvió cuando cambié a la estimación MCMC usando el mismo software. Incluso he tenido una correlación superior a 3 que se resolvió cuando cambié la estimación. Usando IGLS, la correlación se calcula después de la estimación como la covarianza dividida por el producto de la raíz cuadrada del producto de las varianzas asociadas, y esto no tiene en cuenta la incertidumbre en cada una de las estimaciones constituyentes.

El software IGLS no 'sabe' que la covarianza implica una correlación y solo calcula estimaciones de una función de varianza constante, lineal, cuadrática, etc. Por el contrario, el enfoque MCMC se basa en la suposición de muestras de una distribución normal multivariada que corresponde a variaciones y covarianzas con buenas propiedades y propagación de errores completos, de modo que la incertidumbre en la estimación de las covarianzas se tiene en cuenta en la estimación de las variaciones. y viceversa.

MLwin posee la cadena de estimación MCMC con estimaciones IGLS y la matriz de covarianza de varianza definida no negativa podría necesitar ser alterada cambiando la covarianza a cero desde el principio antes de comenzar el muestreo.

Para un ejemplo trabajado ver

Desarrollar modelos multinivel para analizar la contextualidad, la heterogeneidad y el cambio utilizando MLwiN 3, Volumen 1 (actualizado en septiembre de 2017); El volumen 2 también está en RGate

https://www.researchgate.net/publication/320197425_Vol1Training_manualRevisedSept2017

Apéndice del Capítulo 10