¿Cómo simular modelos de efectos aleatorios?

Question

Simular modelos lineales es bastante sencillo:

set.seed(42)    
years <- rnorm(100, 12, 8)
work_hours <- rnorm(100, 8, 2)
income <- 2*years + 0.5*work_hours + 2000 + rnorm(100, 0, 10)
plot(work_hours, income2)
lmmodel <- lm(income ~ years + work_hours)
summary(lmmodel)

O modelos logísticos:

set.seed(42)
x1 <-  rnorm(100) 
x2 <-  rnorm(100)
z <- 1 + 2*x1 + 3*x2   
pr <- 1/(1+exp(-z)) 

y = rbinom(100,1,pr) 

df <- data.frame(y=y,x1=x1,x2=x2)
logitmodel <- glm( y~x1+x2,data=df,family="binomial")
summary(logitmodel)

Entonces, ¿cómo se simulan los modelos de efectos aleatorios? Hay muchos "sabores" con esta clase de modelos. Mirando el [libro] de Faraway[1] los hay:

• Bloques como efectos aleatorios
• Parcelas divididas
• Efectos anidados
• Efectos cruzados
• Modelos multinivel
• Medidas repetidas
• Datos longitudinales/de panel
• Modelos de efectos mixtos para respuestas no normales

¿Cómo podría simularlos para poder jugar con ellos?

[1]: Ampliación del modelo lineal con R - John Faraway

Answer 1

Así es como simulo los efectos aleatorios. Voy a demostrar para la regresión lineal, pero extenderlo a un GLM diferente debe ser sencillo.

Empecemos con un modelo de intercepción aleatoria. El modelo suele escribirse como

$$ y = XB + Z\gamma $$

Dónde $Z$ es un indicador del grupo y $\gamma_i$ se distribuye normalmente con media 0 y cierta varianza. La simulación de este modelo es la siguiente...

groups<-1:5
N <- 250
g <- factor(sample(groups, replace = T, size = N), levels = groups)
x <- rnorm(N)
X <- model.matrix(~x)
Z <- model.matrix(~g-1)

beta <- c(10, 2)
gamma <- rnorm(length(groups), 0, 0.25)

y = X %*% beta + Z%*% gamma + rnorm(N, 0, 0.3)

Adaptemos un modelo mixto y veamos si recuperamos algunas de estas estimaciones

library(lme4)

model = lmer(y~x + (1|g), data = d)

summary(model)
inear mixed model fit by REML ['lmerMod']
Formula: y ~ x + (1 | g)
   Data: d

REML criterion at convergence: 136.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.85114 -0.65429 -0.00888  0.65268  2.63459 

Random effects:
 Groups   Name        Variance Std.Dev.
 g        (Intercept) 0.05771  0.2402  
 Residual             0.09173  0.3029  
Number of obs: 250, groups:  g, 5

Fixed effects:
            Estimate Std. Error t value
(Intercept)  9.95696    0.10914   91.23
x            2.00198    0.01993  100.45

Correlation of Fixed Effects:
  (Intr)
x -0.008

Los efectos fijos tienen buena pinta, y la desviación típica del grupo (0,25) se estima con bastante precisión, al igual que la desviación típica residual.

Los modelos de pendiente aleatoria son similares. Bajo el supuesto de que cada pendiente proviene de una distribución normal, entonces podemos escribir la pendiente como

$$ y = Bx + \beta_i x$$

Aquí $B$ es la media de la población y $\beta_i$ es el efecto del grupo i. He aquí una simulación

library(tidyverse)

groups<-1:5
N <- 250
g <- sample(groups, replace = T, size = N)
x <- rnorm(N)
X <- model.matrix(~x)

B <- c(10, 2)
beta <- rnorm(length(groups), 0, 2)

y = X %*% B + x*beta[g] + rnorm(N, 0, 0.3)

y un modelo...

library(lme4)

d = tibble(y, x, g)

model = lmer(y ~ x + (x|g), data = d)

summary(model)

Linear mixed model fit by REML ['lmerMod']
Formula: y ~ x + (x | g)
   Data: d

REML criterion at convergence: 158.9

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.95141 -0.65904  0.02218  0.61932  2.66614 

Random effects:
 Groups   Name        Variance  Std.Dev. Corr
 g        (Intercept) 2.021e-05 0.004496     
          x           3.416e+00 1.848314 1.00
 Residual             9.416e-02 0.306856     
Number of obs: 250, groups:  g, 5

Fixed effects:
            Estimate Std. Error t value
(Intercept) 10.00883    0.01984  504.47
x            2.05913    0.82682    2.49

Correlation of Fixed Effects:
  (Intr)
x 0.099

Estos son los coeficientes de los 5 grupos

coef(model)
$g
  (Intercept)         x
1    10.00135 -1.015180
2    10.01335  3.919787
3    10.00934  2.270760
4    10.01081  2.873636
5    10.00928  2.246626

y compararlos con los valores verdaderos

B[2] + beta

-0.9406479  3.9195119  2.2976457  2.8536623  2.3539863

```

Answer 2

Basta con escribir una fórmula (algebraica) para el modelo y simular a partir de esa descripción. Pondré un ejemplo muy sencillo, un modelo con múltiples observaciones de los mismos sujetos, con una estructura de covarianza intercambiable. Dicha estructura puede representarse con un intercepto aleatorio para cada sujeto. También una covariable a nivel de sujeto: $$ y_{ij}=\mu + \alpha x_i + \epsilon_i + \epsilon_{ij} $$ para $i=1,2,\dotsc,n$ y $j=1,\dotsc,k$ dentro de cada asignatura. Se trata, pues, de un modelo equilibrado. El mismo principio se utiliza para las situaciones desequilibradas, pero eso da más programación. A continuación, debemos especificar los valores de los parámetros fijos y las distribuciones de los efectos aleatorios $\epsilon_i, \epsilon_{ij}$ . Un simple código R es:

N <- 20 # Number subjects
k <- 4  # Number obs within subject
set.seed(7*11*13) # My public seed

id <- as.factor(1:N)
x <-  runif(N, 1, 5)
idran <- rnorm(N, 0, 1)
obsran <- rnorm(N*k, 0, 2)
mu <- 10.
alpha <- 1.

X <- rep(x, each=k)
Y <- mu + alpha*X + rep(idran, each=k) + obsran

Un gráfico de estos datos simulados es:

Para situaciones más complejas ayudaría algún paquete preprogramado, existe un paquete simstudy en CRAN que puede ayudar. Véase también Matrices de modelos de efectos mixtos y https://stackoverflow.com/questions/30896540/extract-raw-model-matrix-of-random-effects-from-lmer-objects-lme4-r , https://stackoverflow.com/questions/55199251/how-to-create-a-simulation-of-a-small-data-set-in-r .