Comparar un modelo mixto (sujeto como efecto aleatorio) con un modelo lineal simple (sujeto como efecto fijo)

Estoy terminando un análisis de un gran conjunto de datos. Me gustaría tomar el modelo lineal utilizado en la primera parte del trabajo y volver a ajustarlo usando un modelo lineal mixto (LME). El LME sería muy similar con la excepción de que una de las variables utilizadas en el modelo se utilizaría como un efecto aleatorio. Estos datos provienen de muchas observaciones (> 1000) en un pequeño grupo de sujetos (~ 10) y sé que modelar el efecto del sujeto se realiza mejor como un efecto aleatorio (esta es una variable que quiero cambiar). El código R se vería así:

my_modelB <- lm(formula = A ~ B + C + D)    
lme_model <- lme(fixed=A ~ B + C, random=~1|D, data=my_data, method='REML')

Todo funciona bien y los resultados son muy similares. Sería bueno si pudiera usar algo como RLRsim o un AIC / BIC para comparar estos dos modelos y decidir cuál es el más apropiado. Mis colegas no quieren informar sobre la LME porque no hay una forma fácilmente accesible de elegir cuál es "mejor", aunque creo que la LME es el modelo más apropiado. ¿Alguna sugerencia?

r regression hypothesis-testing mixed-model lme4-nlme MudPhud
fuente

Respuestas:

Esto es para agregar a la respuesta de @ ocram porque es demasiado largo para publicar como comentario. Lo trataría A ~ B + Ccomo su modelo nulo para que pueda evaluar la importancia estadística de una Dintercepción aleatoria de nivel en una configuración de modelo anidado. Como señaló Oram, las condiciones de regularidad se violan cuando , y el estadístico de prueba de razón de probabilidad (LRT) no necesariamente se distribuirá asintóticamente . La solución que me enseñaron fue arrancar el LRT (cuya distribución de arranque probablemente no será ) paramétricamente y calcular un valor p de arranque como este: $H_0: \sigma^2 = 0$ $\chi^2$ $\chi^2$

library(lme4)
my_modelB <- lm(formula = A ~ B + C)
lme_model <- lmer(y ~ B + C + (1|D), data=my_data, REML=F)
lrt.observed <- as.numeric(2*(logLik(lme_model) - logLik(my_modelB)))
nsim <- 999
lrt.sim <- numeric(nsim)
for (i in 1:nsim) {
    y <- unlist(simulate(mymodlB))
    nullmod <- lm(y ~ B + C)
    altmod <- lmer(y ~ B + C + (1|D), data=my_data, REML=F)
    lrt.sim[i] <- as.numeric(2*(logLik(altmod) - logLik(nullmod)))
}
mean(lrt.sim > lrt.observed) #pvalue

La proporción de LRT bootstrapped más extrema que la LRT observada es el valor p.

bloqueado
fuente

Gracias por completar mi respuesta. Además, a veces las personas usan una mezcla de chi-cuadrados en lugar de una distribución de chi-cuadrado para la estadística de prueba.

ocram

@ocram +1 por su comentario sobre la decisión de tratar la variable como aleatoria o fija por separado del análisis. @MudPhud Si su PI no comprende el problema e insiste en un valor p, entonces tal vez solo le muestre el resultado de la prueba del efecto aleatorio (que incluiría de todos modos en la redacción).

cerrado el

Gracias por el codigo. Cuando me encontré con que el resultado es que ninguno de los LRTs bootstrap son mayores que la observada, por lo que esto significa que puedo adhieren a la película sin los efectos aleatorios o incluso la variable original tirado.

MudPhud

@MudPhud: ¿Recibió algún error? Intente escribir lrt.simpara asegurarse de que no todos sean ceros, en cuyo caso el culpable más probable sería que no tenga el paquete lme4instalado.

cerrado el

No son 0, solo son muy pequeños (~ 1e-6) en comparación con los observados (63.95).

MudPhud

No estoy totalmente seguro de averiguar qué modelo se ajusta cuando utiliza la función lme. (Supongo que se supone que el efecto aleatorio sigue una distribución normal con media cero?). Sin embargo, el modelo lineal es un caso especial del modelo mixto cuando la varianza del efecto aleatorio es cero. Aunque existen algunas dificultades técnicas (porque está en el límite del espacio de parámetros para la varianza), debería ser posible probar vs ... $0$ $H_0:variance = 0$ $H_1: variance > 0$

EDITAR

Para evitar confusiones: la prueba mencionada anteriormente se usa a veces para decidir si el efecto aleatorio es significativo ... pero no para decidir si se debe transformar o no en un efecto fijo.

ocram
fuente

La pregunta es: ¿hay pruebas para decidir si la variable debe modelarse como un efecto mixto o aleatorio? De lo contrario, podría hacer la prueba que describió y luego probarla con un chi-cuadrado dist (no estoy seguro de cuál sería la prueba adecuada).

MudPhud

@MudPhud: modelar una variable como un efecto fijo o aleatorio debería decidirse antes del análisis, cuando se planifica el estudio. Depende, en particular, del alcance de sus conclusiones. Los efectos aleatorios permiten una mayor generalización. También podría evitar algunas dificultades técnicas. Por ejemplo, las asintóticas pueden romperse cuando crece el número de parámetros, como es el caso cuando una variable categórica con muchos niveles se considera una variable fija.

ocram

Estoy de acuerdo, pero cuando traté de explicarle esto a mi PI, él simplemente se dio la vuelta y pidió un valor p de algún tipo. Quiero incluir este análisis en un manuscrito, pero él no lo incluirá si no hay una justificación más concreta.

MudPhud

@MudPhud: Que yo sepa, no existe un valor p para tal decisión. Si el interés se centra en el efecto de los niveles específicos elegidos, entonces debe considerarse como fijo. Si los niveles de factores disponibles se ven como una muestra aleatoria de una población más grande y se desean inferencias para la población más grande, el efecto debe ser aleatorio.

ocram