En un modelo multinivel, ¿cuáles son las implicaciones prácticas de estimar los parámetros de correlación de efectos aleatorios versus no estimar?

En un modelo multinivel, ¿cuáles son las implicaciones prácticas y relacionadas con la interpretación de estimar los parámetros de correlación de efectos aleatorios versus no estimar? La razón práctica para preguntar esto es que en el marco de Imer en R, no hay un método implementado para estimar los valores de p mediante técnicas MCMC cuando las estimaciones se realizan en el modelo de las correlaciones entre parámetros.

Por ejemplo, mirando este ejemplo (porciones citadas a continuación), ¿cuáles son las implicaciones prácticas de M2 ​​versus M3? Obviamente, en un caso P5 no se estimará y en el otro sí.

Preguntas

1. Por razones prácticas (el deseo de obtener un valor p a través de las técnicas MCMC), uno podría querer ajustar un modelo sin correlaciones entre los efectos aleatorios, incluso si P5 es sustancialmente distinto de cero. Si uno hace esto, y luego estima los valores p a través de la técnica MCMC, ¿son los resultados interpretables? (Sé que @Ben Bolker ha mencionado anteriormente que "combinar pruebas de significación con MCMC es un poco incoherente, estadísticamente, aunque entiendo la necesidad de hacerlo (obtener intervalos de confianza es más compatible)" , así que si te hará dormir mejor por la noche finjo que dije intervalos de confianza.)
2. Si uno no puede estimar P5, ¿es lo mismo que afirmar que es 0?
3. Si P5 realmente no es cero, ¿de qué manera se ven afectados los valores estimados de P1-P4?
4. Si P5 realmente no es cero, ¿de qué manera se ven afectadas las estimaciones de error para P1-P4?
5. Si P5 realmente no es cero, ¿de qué manera las interpretaciones de un modelo que no incluye P5 son defectuosas?

Tomando prestado de la respuesta de @Mike Lawrence (los que tienen más conocimiento que yo son libres de reemplazar esto con la notación de modelo completo, no estoy completamente seguro de poder hacerlo con una fidelidad razonable):

M2: V1 ~ (1|V2) + V3 + (0+V3|V2)(Estimaciones P1 - P4)

M3: V1 ~ (1+V3|V2) + V3(Estimaciones P1-P5)

Parámetros que pueden estimarse:

P1 : una intercepción global

P2 : intercepciones de efecto aleatorio para V2 (es decir, para cada nivel de V2, la desviación de la intercepción de ese nivel de la intercepción global)

P3 : una estimación global única para el efecto (pendiente) de V3

P4 : El efecto de V3 dentro de cada nivel de V2 (más específicamente, el grado en que el efecto V3 dentro de un nivel dado se desvía del efecto global de V3), al tiempo que impone una correlación cero entre las desviaciones de intercepción y las desviaciones del efecto V3 entre niveles de V2.

P5 : La correlación entre las desviaciones de intercepción y las desviaciones de V3 entre los niveles de V2

Serían aceptables las respuestas derivadas de una simulación suficientemente grande y amplia junto con el código adjunto en R usando lmer.

Considere los datos del estudio del sueño, incluidos en lme4. Bates analiza esto en su libro en línea sobre lme4. En el capítulo 3, considera dos modelos para los datos.

y

En el estudio participaron 18 sujetos, estudiados durante un período de 10 días privados de sueño. Los tiempos de reacción se calcularon al inicio y en los días posteriores. Existe un claro efecto entre el tiempo de reacción y la duración de la privación del sueño. También hay diferencias significativas entre los sujetos. El modelo A permite la posibilidad de una interacción entre la intercepción aleatoria y los efectos de la pendiente: imagine, por ejemplo, que las personas con tiempos de reacción pobres sufren de manera más aguda los efectos de la privación del sueño. Esto implicaría una correlación positiva en los efectos aleatorios.

En el ejemplo de Bates, no hubo una correlación aparente del gráfico Lattice y no hubo diferencias significativas entre los modelos. Sin embargo, para investigar la pregunta planteada anteriormente, decidí tomar los valores ajustados del estudio del sueño, aumentar la correlación y observar el rendimiento de los dos modelos.

Como puede ver en la imagen, los largos tiempos de reacción están asociados con una mayor pérdida de rendimiento. La correlación utilizada para la simulación fue 0.58

Simulé 1000 muestras, usando el método de simulación en lme4, basado en los valores ajustados de mis datos artificiales. Encajé a M0 y Ma en cada uno y miré los resultados. El conjunto de datos original tenía 180 observaciones (10 para cada uno de los 18 sujetos), y los datos simulados tienen la misma estructura.

La conclusión es que hay muy poca diferencia.

1. Los parámetros fijos tienen exactamente los mismos valores en ambos modelos.
2. Los efectos aleatorios son ligeramente diferentes. Hay 18 efectos aleatorios de intercepción y 18 pendientes para cada muestra simulada. Para cada muestra, estos efectos se ven obligados a sumar 0, lo que significa que la diferencia media entre los dos modelos es (artificialmente) 0. Pero las variaciones y covarianzas difieren. La mediana de covarianza bajo MA fue 104, contra 84 bajo M0 (valor real, 112). Las variaciones de pendientes e intersecciones fueron mayores bajo M0 que MA, presumiblemente para obtener el margen de maniobra adicional necesario en ausencia de un parámetro de covarianza libre.
3. El método ANOVA para lmer proporciona una estadística F para comparar el modelo de pendiente con un modelo con solo una intercepción aleatoria (sin efecto debido a la privación del sueño). Claramente, este valor era muy grande en ambos modelos, pero era típicamente (pero no siempre) mayor en MA (media 62 frente a media de 55).
4. La covarianza y la varianza de los efectos fijos son diferentes.
5. Aproximadamente la mitad del tiempo, sabe que la MA es correcta. La mediana del valor p para comparar M0 con MA es 0.0442. A pesar de la presencia de una correlación significativa y 180 observaciones balanceadas, el modelo correcto se elegiría solo la mitad del tiempo.
6. Los valores pronosticados difieren bajo los dos modelos, pero muy ligeramente. La diferencia media entre las predicciones es 0, con un SD de 2.7. El SD ​​de los valores predichos mismos es 60.9

Entonces, ¿por qué pasa ésto? @gung supuso, razonablemente, que no incluir la posibilidad de una correlación obliga a que los efectos aleatorios no estén correlacionados. Quizás debería; pero en esta implementación, los efectos aleatorios pueden correlacionarse, lo que significa que los datos pueden llevar los parámetros en la dirección correcta, independientemente del modelo. La incorrección del modelo incorrecto se muestra en la probabilidad, por lo que puede (a veces) distinguir los dos modelos en ese nivel. El modelo de efectos mixtos básicamente ajusta regresiones lineales a cada sujeto, influenciado por lo que el modelo cree que deberían ser. El modelo incorrecto fuerza el ajuste de valores menos plausibles de los que obtiene con el modelo correcto. Pero los parámetros, al final del día, se rigen por el ajuste a los datos reales.

Aquí está mi código algo torpe. La idea era ajustar los datos del estudio del sueño y luego construir un conjunto de datos simulados con los mismos parámetros, pero una correlación mayor para los efectos aleatorios. Ese conjunto de datos se alimentó para simular.lmer () para simular 1000 muestras, cada una de las cuales se ajustó en ambos sentidos. Una vez que emparejé los objetos ajustados, pude extraer diferentes características del ajuste y compararlos, usando pruebas t, o lo que sea.

    # Fit a model to the sleep study data, allowing non-zero correlation
fm01 <- lmer(Reaction ~ 1 + Days +(1+Days|Subject), data=sleepstudy, REML=FALSE)
# Now use this to build a similar data set with a correlation = 0.9
# Here is the covariance function for the random effects
# The variances come from the sleep study. The covariance is chosen to give a larger correlation
sigma.Subjects <- matrix(c(565.5,122,122,32.68),2,2) 
# Simulate 18 pairs of random effects
ranef.sim <- mvrnorm(18,mu=c(0,0),Sigma=sigma.Subjects)
# Pull out the pattern of days and subjects.
XXM <- model.frame(fm01) 
n <- nrow(XXM) # Sample size
# Add an intercept to the model matrix.
XX.f <- cbind(rep(1,n),XXM[,2])
# Calculate the fixed effects, using the parameters from the sleep study. 
yhat <- XX.f %*%  fixef(fm01 )
# Simulate a random intercept for each subject
intercept.r <- rep(ranef.sim[,1], each=10) 
# Now build the random slopes
slope.r <- XXM[,2]*rep(ranef.sim[,2],each=10)
# Add the slopes to the random intercepts and fixed effects
yhat2 <- yhat+intercept.r+slope.r
# And finally, add some noise, using the variance from the sleep study
y <- yhat2 + rnorm(n,mean=0,sd=sigma(fm01))
# Here is new "sleep study" data, with a stronger correlation.
new.data <- data.frame(Reaction=y,Days=XXM$Days,Subject=XXM$Subject)
# Fit the new data with its correct model
fm.sim <- lmer(Reaction ~ 1 + Days +(1+Days|Subject), data=new.data, REML=FALSE)
# Have a look at it
xyplot(Reaction ~ Days | Subject, data=new.data, layout=c(6,3), type=c("p","r"))
# Now simulate 1000 new data sets like new.data and fit each one
# using the right model and zero correlation model.
# For each simulation, output a list containing the fit from each and
# the ANOVA comparing them.
n.sim <- 1000
    sim.data <- vector(mode="list",)
    tempReaction <- simulate(fm.sim, nsim=n.sim)
    tempdata <- model.frame(fm.sim)
    for (i in 1:n.sim){
        tempdata$Reaction <- tempReaction[,i]
			output0 <- lmer(Reaction ~ 1 + Days +(1|Subject)+(0+Days|Subject), data = tempdata, REML=FALSE)
			output1 <- lmer(Reaction ~ 1 + Days +(Days|Subject), data=tempdata, REML=FALSE)
			temp <- anova(output0,output1)
			pval <- temp$`Pr(>Chisq)`[2]
        sim.data[[i]] <- list(model0=output0,modelA=output1, pvalue=pval)
    }

Placidia ya ha proporcionado una respuesta exhaustiva utilizando datos simulados basados ​​en el sleepstudyconjunto de datos. Aquí hay otra respuesta (menos rigurosa) que también usa los sleepstudydatos.

Vemos que uno puede afectar la correlación estimada entre la intersección aleatoria y la pendiente aleatoria al "desplazar" la variable predictora aleatoria. Mira los resultados de los modelos fm1y fm2abajo:

library(lmer)

#Fit Models
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
k <- 3 # Shift "Days" by an arbitrary amount
fm2 <- lmer(Reaction ~ I(Days + k) + (I(Days + k)| Subject), sleepstudy)

fm1 # Model Output
# Linear mixed model fit by REML ['lmerMod']
# Formula: Reaction ~ Days + (Days | Subject)
# Data: sleepstudy
# REML criterion at convergence: 1743.628
# Random effects:
#   Groups   Name        Std.Dev. Corr
# Subject  (Intercept) 24.740       
# Days         5.922   0.07
# Residual             25.592       
# Number of obs: 180, groups:  Subject, 18
# Fixed Effects:
#   (Intercept)         Days  
# 251.41        10.47

fm2 # Model Output
# Linear mixed model fit by REML ['lmerMod']
# Formula: Reaction ~ I(Days + k) + (I(Days + k) | Subject)
# Data: sleepstudy
# REML criterion at convergence: 1743.628
# Random effects:
#   Groups   Name        Std.Dev. Corr 
# Subject  (Intercept) 29.498        
# I(Days + k)  5.922   -0.55
# Residual             25.592        
# Number of obs: 180, groups:  Subject, 18
# Fixed Effects:
#   (Intercept)  I(Days + k)  
# 220.00        10.47

# Random effects from both models
cbind(ranef(fm1)$Subject,ranef(fm2)$Subject)
# (Intercept)        Days (Intercept) I(Days + k)
# 308   2.2585654   9.1989719 -25.3383538   9.1989727
# 309 -40.3985769  -8.6197032 -14.5394628  -8.6197043
# 310 -38.9602458  -5.4488799 -22.6136027  -5.4488807
# 330  23.6904985  -4.8143313  38.1334933  -4.8143315
# 331  22.2602027  -3.0698946  31.4698868  -3.0698946
# 332   9.0395259  -0.2721707   9.8560377  -0.2721706
# 333  16.8404311  -0.2236244  17.5113040  -0.2236243
# 334  -7.2325792   1.0745761 -10.4563076   1.0745761
# 335  -0.3336958 -10.7521591  31.9227854 -10.7521600
# 337  34.8903508   8.6282840   9.0054946   8.6282850
# 349 -25.2101104   1.1734142 -28.7303527   1.1734141
# 350 -13.0699567   6.6142050 -32.9125736   6.6142054
# 351   4.5778352  -3.0152572  13.6236077  -3.0152574
# 352  20.8635924   3.5360133  10.2555505   3.5360138
# 369   3.2754530   0.8722166   0.6588028   0.8722167
# 370 -25.6128694   4.8224646 -40.0802641   4.8224648
# 371   0.8070397  -0.9881551   3.7715053  -0.9881552
# 372  12.3145393   1.2840297   8.4624492   1.2840300

A partir del resultado del modelo, vemos que la correlación de varianza aleatoria ha cambiado. Sin embargo, las pendientes (fijas y aleatorias) se mantuvieron igual, al igual que la estimación de la varianza residual. Las estimaciones de intersecciones (fijas y aleatorias) cambiaron en respuesta a la variable desplazada.

La covarianza de la pendiente de intercepción aleatoria descorrelacionada para LMM se analiza en las notas de clase del Dr. Jack Weiss aquí . Weiss señala que reducir la correlación de varianza de esta manera a veces puede ayudar con la convergencia del modelo, entre otras cosas.

El ejemplo anterior varía la correlación aleatoria (parámetro "P5"). Dirigiendo parcialmente el Q3 del OP, vemos en el resultado anterior que:

#   Parameter           Status
=================================
P1  Fixed Intercept     Affected
P2  Random Intercepts   Affected
P3  Fixed Slope         Not Affected
P4  Random Slopes       Not Affected
P5  Random Correlation  Affected