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

Question

En un modelo multinivel, ¿cuáles son la práctica y la interpretación relacionados con las implicaciones de la estimación, o no, la estimación de efectos aleatorios de correlación de los parámetros? La razón práctica para preguntar esto, es que en el lmer marco en R, no hay ningún método implementado para la estimación de los valores de p a través de la MCMC técnicas cuando se hacen estimaciones en el modelo de las correlaciones entre los parámetros.

Por ejemplo, buscando en este ejemplo (partes en que se cita a continuación), ¿cuáles son las implicaciones prácticas de M2 frente M3. Obviamente, en un caso P5 no será estimado y en el otro se va.

Preguntas

1. Por razones prácticas (el deseo de obtener un p-valor a través de la MCMC técnicas) que uno puede querer ajustar un modelo sin las correlaciones entre los efectos aleatorios incluso si P5 es sustancialmente distinto de cero. Si uno hace esto, y, a continuación, las estimaciones de los valores de p a través de la MCMC técnica, son los resultados interpretables? (Sé que @Ben Bolker ha dicho anteriormente que "la combinación de pruebas de significación con MCMC es un poco incoherente, estadísticamente, aunque entiendo la necesidad de hacerlo (obtención de intervalos de confianza es más compatible)", por lo que si va a hacer dormir mejor por la noche pretender que dijo que los intervalos de confianza.)
2. Si falla uno de ellos para estimar P5, es que el mismo que afirmar que es 0?
3. Si P5 realmente no es cero, entonces ¿en qué son los valores estimados de P1-P4 afectados?
4. Si P5 realmente no es cero, entonces ¿de qué manera son las estimaciones de error para P1-P4 afectados?
5. Si P5 realmente no es cero, entonces ¿de qué maneras son interpretaciones de un modelo de no incluir P5 defectos?

El endeudamiento de @Mike Lawrence respuesta del (los más conocedores de lo que estoy libre para reemplazar por completo el modelo de notación, no estoy completamente seguro de que puede hacerlo con razonable fidelidad):

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

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

Los parámetros que pueden ser estimados:

P1: global interceptar

P2: efecto Aleatorio intercepta para V2 (es decir, para cada nivel de V2, que el nivel de la intercepción de la desviación de la global interceptar)

P3: Una sola estimación global del efecto (pendiente) de la V3

P4: el efecto de La V3 dentro de cada nivel de V2 (más específicamente, el grado en que la V3 efecto dentro de un nivel dado se desvía del efecto global de la V3), mientras que el cumplimiento de una correlación cero entre la intersección de las desviaciones y V3 efecto de las desviaciones en los niveles de V2.

P5: la correlación entre La intercepción de las desviaciones y V3 desviaciones a través de los niveles de V2

Respuestas derivadas de una lo suficientemente grande y amplio de simulación, junto con el código de acompañamiento en R usando lmer sería aceptable.

Answer 1

Considerar la sleepstudy de datos, incluidos en lme4. Bates se explica esto en su libro sobre lme4. En el capítulo 3, se considera que dos de los modelos para los datos.

$$M0: \textrm{Reacción} \sim 1 + \textrm{Días} + (1|\textrm{Asunto}) +(0+\textrm{Días}|\textrm{Asunto}) $$

y

$$MA: \textrm{Reacción} \sim 1 + \textrm{Días} + (\textrm{Días}|\textrm{Asunto}) $$

En el estudio participaron 18 sujetos, estudió durante un período de 10 privados de sueño los días. Los tiempos de reacción se calcula en la línea de base y en los días sucesivos. Hay 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. Modelo permite la posibilidad de una interacción entre el azar intercepto y de la pendiente efectos: imaginar, decir, que las personas con buenos tiempos de reacción sufren de manera más aguda a partir de los efectos de la privación del sueño. Esto implicaría una correlación positiva en el de efectos aleatorios.

En Bates ejemplo, no hubo ninguna correlación aparente de la Celosía de la trama y que no hay diferencia significativa entre los modelos. Sin embargo, para investigar la pregunta planteada anteriormente, me decidí a tomar los valores ajustados de la sleepstudy, poner encima de la correlación y la mirada en el rendimiento de los dos modelos.

Como se puede ver en la imagen, largos tiempos de reacción se asocia con una mayor pérdida de rendimiento. La correlación se utiliza para la simulación fue de 0,58

He simulado 1000 muestras, utilizando el simulador método en lme4, basado en los valores ajustados de mis datos artificiales. Me caben M0 y Ma a cada uno y miramos los resultados. El conjunto de datos original había 180 observaciones (10 para cada uno de los 18 sujetos), y los datos simulados tiene la misma estructura.

La conclusión es que hay muy poca diferencia.

1. Fija los parámetros tienen exactamente los mismos valores en ambos modelos.
2. Los efectos aleatorios son ligeramente diferentes. Hay 18 interceptar y 18 de la pendiente de efectos aleatorios para cada una de las muestras simuladas. Para cada muestra, estos efectos se ven obligados a añadir a 0, lo que significa que la diferencia de medias entre los dos modelos es (artificialmente) 0. Pero las varianzas y covarianzas diferentes. La mediana de la covarianza bajo MA fue de 104, contra 84 bajo M0 (valor real, 112). Las desviaciones de las pendientes y las intersecciones fueron mayores bajo M0 de MA, presumiblemente para obtener el extra de margen necesario en la ausencia de una covarianza parámetro.
3. El ANOVA método para lmer da una estadística F para la comparación de la Pendiente del modelo a un modelo con sólo un azar de intersección (ningún efecto debido a la privación del sueño). Claramente, este valor fue muy grandes bajo ambos modelos, pero era normalmente (pero no siempre) más bajo MA (media de 62 vs media de 55).
4. La covarianza y la varianza de los efectos fijos son diferentes.
5. Alrededor de la mitad del tiempo, sabe que MA es correcta. La mediana del valor de p para la comparación de M0 a MA es 0.0442. A pesar de la presencia de una significativa correlación y 180 equilibrada observaciones, el modelo correcto sería elegido sólo alrededor de la mitad del tiempo.
6. Los valores predichos difieren en los dos modelos, pero muy ligeramente. La diferencia media entre las predicciones es 0, con desviación estándar de 2.7. La sd de la predicción de los valores propios es 60.9

Entonces, ¿por qué sucede esto? @gung suponer, razonablemente, que el hecho de no incluir la posibilidad de una correlación de fuerzas de los efectos al azar para ser correlacionadas. Quizás debería, pero en esta implementación, el de efectos aleatorios, se les permite estar correlacionados, lo que significa que los datos son capaces de tirar de los parámetros en la dirección correcta, sin importar el modelo. La falsedad del mal modelo se muestra en la probabilidad, por lo que se puede (a veces) distinguir los dos modelos en ese nivel. El modelo de efectos mixtos es básicamente ajuste de regresiones lineales para cada sujeto, influenciado por lo que el modelo piensa que debería ser. El modelo equivocado de las fuerzas del fit de menos plausible valores que usted recibe bajo el modelo adecuado. Pero los parámetros, al final del día, se rigen por el ajuste a los datos reales.

Aquí está mi un poco torpe código. La idea era para que se ajuste el estudio del sueño de datos y, a continuación, construir un conjunto de datos simulados con los mismos parámetros, pero una mayor correlación de los efectos aleatorios. Ese conjunto de datos se alimentó a simular.lmer() para simular 1000 muestras, cada una de las cuales estaba a la altura de las dos maneras. Una vez que se había apareado equipada objetos, pude sacar diferentes características de la adaptación y la comparación de ellos, utilizando la prueba 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)
    }