Regresión lineal multinivel / de efectos mixtos / jerárquica. ¿Modelo y consejos de interpretación?

Question

Quiero ejecutar y luego reportar un modelo multinivel sobre los datos de ChickWeight en R.

He realizado varios análisis sobre los que tengo algunas preguntas. Me han dicho que el AIC es una buena forma de comprobar el rendimiento predictivo fuera de muestra de un modelo.

He considerado los siguientes modelos, y el de abajo es el que tiene el AIC más bajo, lo cual, si no me equivoco, ¿es bueno?

library(tidyverse)
library(lme4)
data("ChickWeight", package = "datasets")

chickm <- lmer(weight ~ Time + (1 | Chick), data = ChickWeight)
chickm2 <- lmer(weight ~ Time + (Time | Chick), data = ChickWeight)
chickm3 <- lmer(weight ~ Time + Diet + (Time | Chick), data = ChickWeight)
chickm4 <- lmer(weight ~ Time * Diet + (Time | Chick), data = ChickWeight)

Quiero informar de chickm4 ya que el efecto de interacción parece mejorar el modelo. Sin embargo, chickm3 también sería aceptable, ya que me gustaría incluir Diet.

La siguiente pregunta es sobre el intercepto que producen los modelos. El intercepto real (peso del polluelo cuando los días = 0 aka nace, es ~ 41g). Sin embargo chickm3 da el intercepto como 26g y chickm4 como 33g. Se puede ver en el gráfico de abajo que las líneas lineales de mejor ajuste causan esto. Los interceptos de las líneas están todos por debajo de los puntos verdaderos.

library(tidyverse)

ggplot(ChickWeight, aes(x = Time, y = weight, colour = Diet)) + geom_point() +
  stat_smooth(method = 'lm', se = F) + theme_minimal()

¿Es esto un problema o un ejemplo de un mal modelo, o simplemente un inconveniente de los datos y/o de los modelos lineales en comparación con los no lineales? Tengo que informar de un lineal regresión, así que no puedo cambiar esto.

Pregunta final. chickm3 muestra El modelo no converge con max|grad| = 0.00225537 (tol = 0.002, componente 1)' cuando lo resumo(). ¿Qué significa esto y es malo? Mi modelo chickm4 no hace esto.

Answer 1

Un par de puntos:

• Estos modelos están anidados y, por lo tanto, también puede utilizar la prueba de hipótesis formal mediante la prueba F para compararlos. Estas pruebas F las proporciona el lmerTest paquete.
• En general, no es bueno sólo considerar medidas como el AIC o los valores p para decidir qué términos incluir en el modelo, especialmente los efectos fijos en este caso. También debería ver cómo de sustantiva es la diferencia entre el modelo aditivo chickm3 y el modelo de interacción chickm4 . Podría evaluarlo más fácilmente utilizando un gráfico de efectos - véase, por ejemplo, el efectos o ggeffects paquetes.
• De hecho, el efecto de que el intercepto del modelo no se corresponda con el valor medio que se tiene para Time igual a cero y Diet igual al nivel de referencia es porque se ajusta una recta, y la estimación de la pendiente también afecta a la estimación del intercepto. Se podría mitigar este efecto ajustando un modelo no lineal para la variable Time variable utilizando, por ejemplo, splines.
• Por último, en cuanto al mensaje de advertencia de convergencia, podría ignorarlo en este caso porque está muy cerca del valor de tolerancia.

Answer 2

He hecho algunos modelos de ejemplo que deberían estar anidados si los pruebas en orden. El modelo que creo mejor es chickm6 que:

• elimina el efecto inicial de la dieta en el tiempo 0, lo cual es intuitivo, ya que en el tiempo 0 (imagino que al nacer), la dieta no debería tener ningún efecto sobre el peso de los pollitos.
• ha añadido una spline lineal de dos piezas, con el nudo situado en el Tiempo 6 (algo determinado a ojo por la curva de loess). Esto permite que la respuesta de referencia modelizada en el tiempo 0 (39,8) se acerque mucho más al peso medio observado (41), al tiempo que el modelo sigue siendo simple, explicable y lineal. Por supuesto, se puede probar con más nudos.

library(stats)
library(tidyverse)
library(lme4)
data("ChickWeight", package = "datasets")

ChickWeight <- as_tibble(ChickWeight) %>%
  mutate(Time_6p = pmax(0, Time - 6))

chickm4 <- lmer(weight ~ Time * Diet + (Time | Chick), data = ChickWeight)
chickm5 <- lmer(weight ~ Time + Time_6p + Time:Diet + (Time | Chick), data = ChickWeight)
chickm6 <- lmer(weight ~ Time + Time_6p + Time:Diet + Time_6p:Diet + (Time | Chick), data = ChickWeight)
chickm7 <- lmer(weight ~ Time*Diet + Time_6p*Diet + (Time | Chick), data = ChickWeight)

anova(chickm4, chickm5, chickm6, chickm7) # not sure if it actually does a refit with ML

pred_df <- ChickWeight %>%
  distinct(Time, Time_6p, Diet) %>%
  mutate(pred4 = predict(chickm4, newdata = ., re.form = ~0)) %>%
  mutate(pred6 = predict(chickm6, newdata = ., re.form = ~0))

ChickWeight %>%
  mutate(Chick = forcats::fct_reorder(as.character(Chick), sample(row_number()))) %>% 
  ggplot(aes(x = Time)) + 
  geom_line(aes(y = weight, color = Chick)) + 
  geom_line(data = pred_df, aes(y = pred4, linetype = 'pred4'), size = 1) +
  geom_line(data = pred_df, aes(y = pred6, linetype = 'pred6'), size = 1) + 
  geom_smooth(aes(y = weight), method = 'loess') +
  scale_color_discrete(guide = FALSE) +
  facet_wrap(~Diet)