¿Debo eliminar los interceptos aleatorios de mi modelo?

Question

He recogido algunos datos sobre los tiempos de respuesta (Y) en dos condiciones variables (X1 y X2). Las condiciones son variables continuas, aunque las establezco en valores fijos de 1,2,3,4 y 5.

Tengo 10 sujetos y cada sujeto fue expuesto a cada combinación de X1 y X2.

Así que decidí utilizar un modelo lineal de efectos mixtos:

mod <- lmer(Y ~ X1 * X2 + (X1 * X2 | Subject), dt)

Sin embargo, hubo un ajuste singular. Así que seguí las indicaciones de este post:
Cómo simplificar una estructura aleatoria singular cuando las correlaciones declaradas no están cerca de +1/-1

que dice que hay que ejecutar un análisis de componentes principales en la matriz de varianza-covarianza de los efectos aleatorios, y luego ajustar un modelo de "correlación cero" y buscar una varianza cercana a cero.

Tras seguir estas instrucciones, comprobé que la varianza de los interceptos aleatorios era muy cercana a cero, por lo que eliminé los interceptos aleatorios con:

mod = lmer(Y ~ X1 * X2 + (0 + X1 * X2 | Subject), dt)

y el modelo convergió sin singularidad.

Pero esto parece muy extraño porque creo que necesito controlar las medidas repetidas dentro de cada sujeto con los interceptos aleatorios. ¿Debería ajustar efectos fijos para los sujetos en su lugar?

Answer 1

En primer lugar, siempre es una buena idea trazar los datos:

library(lme4)
library(tidyverse)
dt <- read.csv("https://raw.githubusercontent.com/camhsdoc/data/main/random_slopes.csv")

dt$Subject <- as.factor(dt$Subject)
ggplot(dt, aes(y = Y, x = X1, colour = Subject, group = Subject)) + geom_point() + geom_smooth(method = "lm", se = FALSE)

Así que podemos ver directamente que parece haber variación en los interceptos, aunque el rango de la x empiece en 1.

Así que probablemente haya algo más en estos datos.

Sin embargo, hubo un ajuste singular. Así que seguí la guía de este post que dice que hay que ejecutar un análisis de componentes principales en la matriz de varianza-covarianza de los efectos aleatorios.

...y si miramos la salida de esto vemos:

mod <- lmer(Y ~ X1 * X2 + (X1 * X2 | Subject), dt)
## boundary (singular) fit: see ?isSingular
summary(rePCA(mod))
## $Subject
## Importance of components:
##                          [,1]   [,2]   [,3]      [,4]
## Standard deviation     1.2748 1.1303 0.8397 2.245e-05
## Proportion of Variance 0.4505 0.3541 0.1954 0.000e+00
## Cumulative Proportion  0.4505 0.8046 1.0000 1.000e+00

Así que la varianza se explica por los tres primeros componentes, lo que indica que uno de ellos es innecesario. En realidad, este paso no es necesario, ya que esto es lo que define un ajuste singular, y ya sabemos que es singular.

y luego ajustar un modelo de "correlación cero" y buscar una varianza cercana a cero.

mod1 <- lmer(Y ~ X1 * X2 + (X1 * X2 || Subject), dt)
## boundary (singular) fit: see ?isSingular
summary(mod1)
## 
## Random effects:
##  Groups    Name        Variance  Std.Dev. 
##  Subject   (Intercept) 3.987e-09 6.314e-05
##  Subject.1 X1          1.098e+00 1.048e+00
##  Subject.2 X2          9.007e-01 9.490e-01
##  Subject.3 X1:X2       1.213e+00 1.101e+00
##  Residual              9.995e-01 9.998e-01
## Number of obs: 250, groups:  Subject, 10

y efectivamente vemos que la varianza de los interceptos aleatorios es efectivamente muy pequeña, lo que nos lleva al siguiente modelo:

mod2 <- lmer(Y ~ X1 * X2 + (0 + X1 * X2 || Subject), dt)
summary(mod2)
## 
## Random effects:
##  Groups    Name  Variance Std.Dev.
##  Subject   X1    1.0983   1.0480  
##  Subject.1 X2    0.9007   0.9490  
##  Subject.2 X1:X2 1.2128   1.1013  
##  Residual        0.9995   0.9998  
## Number of obs: 250, groups:  Subject, 10

lo cual, a primera vista, parece estar bien. Sin embargo, hemos visto en el gráfico anterior que hay cierta variación en los interceptos. En mi experiencia he encontrado a menudo problemas de ajuste singular cuando se incluye una interacción como pendiente aleatoria. Así que yo también consideraría este modelo:

mod3 <- lmer(Y ~ X1 * X2 + (X1 + X2 | Subject), dt)
summary(mod3)
##  
## 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr       
##  Subject  (Intercept) 100.671  10.033              
##           X1           10.310   3.211   -0.95      
##           X2           12.375   3.518   -0.97  0.91
##  Residual               5.572   2.360              
## Number of obs: 250, groups:  Subject, 10

Lo destacable de este modelo es que aunque el ajuste no es singular las correlaciones entre los efectos aleatorios son muy altas.

Una forma de avanzar a partir de aquí es utilizar un enfoque de probabilidad para elegir el modelo que "mejor" se ajuste:

anova(mod2, mod3)
## refitting model(s) with ML (instead of REML)
## Data: dt
## Models:
## mod2: Y ~ X1 * X2 + ((0 + X1 | Subject) + (0 + X2 | Subject) + (0 + 
## mod2:     X1:X2 | Subject))
## mod3: Y ~ X1 * X2 + (X1 + X2 | Subject)
##      npar     AIC     BIC  logLik deviance Chisq Df Pr(>Chisq)
## mod2    8  881.84  910.01 -432.92   865.84                    
## mod3   11 1261.61 1300.34 -619.80  1239.61     0  3          1
AIC(mod2); AIC(mod3)
## [1] 883.4465
## [1] 1261.92
BIC(mod2); BIC(mod3)
## [1] 911.6182
## [1] 1300.656

También podemos observar la precisión de la predicción interna, con el RMSE:

(predict(mod2) - dt$Y) ^2 %>% sum() %>% sqrt()
## [1] 14.80808
(predict(mod3) - dt$Y)^2 %>% sum() %>% sqrt()
## [1] 35.1342

En todos los casos, es preferible el modelo con interceptos aleatorios pero sin pendientes aleatorias para la interacción. Dado que ya tenemos razones para querer ajustar los interceptos aleatorios (medidas repetidas) y parece que hay variación en el gráfico, tenemos una razón bastante sólida para elegir este modelo.