La obtención de un nonsence estimación de fijo pendiente después de haber incluido al azar de la pendiente en el modelo de

Question

Supongamos que tengo pseudo-datos, donde se incluyen tres temas y sus tomada analgésico dosis y sentí el dolor en los respectivos cuatro ensayos: (los números no hacen sentido médico, pero quiero que el caso extremo)

library(ggplot2)
library(lme4)

set.seed(1337)
x1 <- 1:4
y1 <- 100 + (rnorm(4) - 8 * x1)
x2 <- 26:29
y2 <- 500 + (rnorm(4) - 7 * x2)
x3 <- 51:54
y3 <- 1000 + (rnorm(4) - 10 * x3)
df_test <- data.frame(subject_id = factor(rep(c(1,2,3), each = 4)),
                      dosage = c(x1, x2, x3),
                      pain = round(c(y1, y2, y3), 1))
df_test

##    subject_id dosage  pain
## 1        1      1     92.2
## 2        1      2     82.6
## 3        1      3     75.7
## 4        1      4     69.6
## 5        2     26    317.3
## 6        2     27    313.0
## 7        2     28    304.9
## 8        2     29    299.1
## 9        3     51    491.9
## 10       3     52    479.6
## 11       3     53    471.0
## 12       3     54    458.3

Puesto que cada sujeto se repitió medido, me gustaría saber,"la relación entre la toma de la dosis y sentí el dolor" por la creación de un modelo mixto donde el sujeto es la agrupación factor de efectos aleatorios.

Dos modelos que he intentado se enumeran a continuación:

Modelo 1:

summary(lmer(pain ~ dosage + (1 | subject_id), data = df_test))

## Linear mixed model fit by REML ['lmerMod']
## Formula: pain ~ dosage + (1 | subject_id)
##    Data: df_test
## 
## REML criterion at convergence: 77.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.60297 -0.19477  0.04491  0.25937  1.53768 
## 
## Random effects:
##  Groups     Name        Variance  Std.Dev.
##  subject_id (Intercept) 161876.29 402.339 
##  Residual                    8.45   2.907 
## Number of obs: 12, groups:  subject_id, 3
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 512.2454   233.2030   2.197
## dosage       -8.1568     0.7489 -10.891
## 
## Correlation of Fixed Effects:
##        (Intr)
## dosage -0.088

Modelo 2:

summary(lmer(pain ~ dosage + (1+dosage | subject_id), data = df_test))

## Linear mixed model fit by REML ['lmerMod']
## Formula: pain ~ dosage + (1 + dosage | subject_id)
##    Data: df_test
## 
## REML criterion at convergence: 101.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.34815 -0.60873 -0.01502  0.61189  1.36244 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev. Corr 
##  subject_id (Intercept) 2526.259 50.262        
##             dosage         1.381  1.175   -1.00
##  Residual                402.622 20.065        
## Number of obs: 12, groups:  subject_id, 3
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  99.4286    34.2614   2.902
## dosage        7.2037     0.7812   9.221
## 
## Correlation of Fixed Effects:
##        (Intr)
## dosage -0.939

He leído las Preguntas acerca de cómo los efectos aleatorios son especificados en lmer. Mi comprensión es que, el Modelo 1 significa que cada sujeto tiene su propia intersección de un dolor o de la dosis, mientras que el Modelo 2, además, permite que cada sujeto tiene su propia pendiente de un dolor o de la dosis. Estoy en lo cierto?

Sin embargo, todavía estoy sorprendido por la enorme diferencia entre los coeficientes de la dosis de efectos fijos (-8.1568 en el Modelo 1 y 7.2037 en el Modelo 2). Yo pensaba que el Modelo 2 es el más adecuado para mi pregunta ya que los sujetos pueden tener varias reacciones a la droga, pero el coeficiente positivo es contradictorio con la tendencia de cada respectivo tema, haciendo la interpretación más razonable.

¿Cómo debo elegir el "correcto" modelo? Y ¿cómo es que los resultados son tan diferentes entre estos dos modelos - ¿cuáles son las razones matemáticas?

Answer 1

El número de sujetos en este ejemplo es bastante baja ($n = 3$). Esto puede conducir a una serie de problemas al tratar de estimación de efectos aleatorios (intercepta o pendientes) porque tenemos suficiente información. Para citar directamente de GLMM preguntas frecuentes sobre el "debo tratar factor xxx como fijos o aleatorios?" tema: "efectos aleatorios estimación está tratando de estimar entre un bloque de varianza". Para estimar con precisión la varianza, con sólo tres temas, que va a ser difícil y puede llevar a resultados incorrectos.

A ese respecto, lo que estamos viendo aquí como un "disparate estimación de fijo pendiente" es en realidad un mínimo local en la estimación ML. Si se trató de un diferente solver (por ejemplo, Nelder-Mead) obtendríamos la "respuesta correcta" (algún valor para el dosage coeficiente de la pendiente cerca de $-8$):

(lmer(pain ~ dosage + (1 +dosage | subject_id), data = df_test, 
      control = lmerControl(optimizer = "Nelder_Mead")))
Linear mixed model fit by REML ['lmerMod']
Formula: pain ~ dosage + (1 + dosage | subject_id)
   Data: df_test
REML criterion at convergence: 68.0749
Random effects:
 Groups     Name        Std.Dev. Corr 
 subject_id (Intercept) 475.649       
            dosage        2.365  -0.81
 Residual                 1.210       
Number of obs: 12, groups:  subject_id, 3
Fixed Effects:
(Intercept)       dosage  
    542.986       -8.215  

Podemos ver que el REML criterio a la convergencia ahora es 68.1 cuando el uso de Nelder-Mead mientras 101.3 cuando se utiliza BOBYQA. Configuración de optCtrl = list(iprint=3)) dentro lmerControl muestra que tanto las rutinas de inicio con un REML criterio de valor de ~106. El BOBYQA procedimiento no mejora mucho en esta partida estimar, pero la Nelder-Mead procedimiento es capaz de conseguir una significativa mejor solución.

Como @ameba señala correctamente, el hecho de que la opción solver bobyqa termina sin ninguna advertencia, es preocupante, pero no quiero ser demasiado rápidos para culpar bobyqa. Esta optimización de los procedimientos hace una aproximación cuadrática en su interior la solución de los pasos y no está diseñado para trabajar con muestras muy pequeñas.