5 Paquete lme4

El paquete lme4 de Bates et al. (2023) es uno de los paquetes más completos para modelos mixtos. Al visitar este enlace se encontrará la página de apoyo del paquete, allí se puede consultar el manual de referencia y las viñetas.

5.1 Función `lmer`

La función lmer es la principal función del paquete lme4. Esta función sirve para ajustar un modelo mixto y su estructura es la siguiente:

lmer(formula, data = NULL, REML = TRUE, control = lmerControl(),
     start = NULL, verbose = 0L, subset, weights, na.action,
     offset, contrasts = NULL, devFunOnly = FALSE, ...)

Los principales argumentos de la función son:

formula: es una fórmula similar a la usada en el modelo lineal clásico. Un ejemplo de fórmula sería y ~ 1 + x1 + x2 + (1 + x2 | grupo) con la cual se indican los efectos fijos y los efectos aleatorios del modelo. Más abajo hay una tabla con más detalles sobre la fórmula.
data: marco de datos donde están las variables.
REML: valor lógico que sirve para indicar si queremos estimaciones maximizando la verosimilitud restringida o la verosimilitud usual.

La siguiente imagen corresponde a la tabla 2 de la viñeta Fitting Linear Mixed-Effects Models using lme4. En la tabla las dos primeras columnas muestran formas equivalentes de incluir las estructuras de modelos mixtos más comunes.

Ejemplo: modelo normal con intercepto aleatorio

En este ejemplo vamos a simular observaciones \(n_i=50\) observaciones para \(G=10\) grupos (en total 500 obs) que tengan la estructura mostrada abajo. El objetivo del ejemplo es ilustrar el uso de la función lmer para estimar los parámetros del modelo.

\[\begin{align*} y_{ij} | b_0 &\sim N(\mu_{ij}, \sigma^2_y) \\ \mu_{ij} &= 4 - 6 x_{ij} + b_{0i} \\ \sigma^2_y &= 16 \\ b_{0} &\sim N(0, \sigma^2_{b0}=625) \\ x_{ij} &\sim U(-5, 6) \end{align*}\]

El vector de parámetros de este modelo es \(\boldsymbol{\Theta}=(\beta_0=4, \beta_1=-6, \sigma_y=4, \sigma_{b0}=25)^\top\).

Solución.

El código para simular las 500 observaciones se muestra a continuación. Observe que se fijó la semilla en dos ocasiones para que el lector pueda replicar el ejemplo y obtener los mismos resultados.

ni <- 50
G <- 10
nobs <- ni * G
grupo <- factor(rep(x=1:G, each=ni))
obs <- rep(x=1:ni, times=G)
set.seed(1234567)
x <- runif(n=nobs, min=-5, max=6)
set.seed(1234567)
b0 <- rnorm(n=G, mean=0, sd=sqrt(625)) # Intercepto aleatorio
b0 <- rep(x=b0, each=ni)               # El mismo intercepto aleatorio pero repetido
media <- 4 - 6 * x + b0
set.seed(1234567)
y <- rnorm(n=nobs, mean=media, sd=sqrt(16))
datos <- data.frame(obs, grupo, b0, x, media, y)

Vamos a explorar los datos simulados.

library(rmarkdown)
paged_table(datos, options = list(rows.print = 6))

El siguiente paso es dibujar los datos para explorar si sería apropiado usar un modelo con intercepto aleatorio (obvio porque así se simularon los datos). El código para dibujar los datos se muestra abajo.

library(ggplot2)

## 
## Attaching package: 'ggplot2'

## The following object is masked from 'package:randomForest':
## 
##     margin

ggplot(datos, aes(x, y, color=grupo) ) + 
  geom_point() + 
  labs(colour="Grupo/Cluster")

En la figura anterior se observa un patrón claro, todas las 10 nubes de puntos tienen la misma pendiente pero diferente intercepto con el eje vertical, eso se debe a que en la simulación se incluyó un \(b_0\).

Para estimar los parámetros del modelos se usa la función mler de la siguiente forma.

library(lme4)
fit1 <- lmer(y ~ x + (1 | grupo), data = datos)

La función summary se puede usar sobre el objeto fit1 para obtener una tabla de resumen, a continuación se ilustra el uso y la salida de summary.

summary(fit1)

Según el resultado anterior \(\hat{\boldsymbol{\Theta}}=(\hat{\beta}_0=2.2378, \hat{\beta}_1=-6.0264, \hat{\sigma}_y=3.9352, \hat{\sigma}_{b0}=25.369)^\top\) mientras que el vector real de parámetros es \(\boldsymbol{\Theta}=(\beta_0=4, \beta_1=-6, \sigma_y=4, \sigma_{b0}=25)^\top\).

Compare los resultados de la tabla anterior obtenida con la función lmer anterior con los resultados obtenidos con la función lme de capítulo 7.

¿Hay alguna similitud?

Ejemplo: recuperando los interceptos aleatorios

¿Cómo se pueden obtener los interceptos aleatorios a partir del modelo ajustado en la sección anterior?

Solución.

Para obtener los interceptos aleatorios se usa la función ranef del paquete lme4 de Bates et al. (2023). A continuación vamos a obtener los interceptos aleatorios y los vamos a comparar con los \(b_0\) simulados.

interceptos_aleatorios <- ranef(fit1)
cbind(interceptos_aleatorios$grupo, b0=unique(b0))

##    (Intercept)         b0
## 1     5.388795   3.917594
## 2    35.790443  34.345280
## 3    19.885504  18.266756
## 4   -31.893039 -33.770023
## 5     1.074942  -0.212874
## 6    10.680482   8.024547
## 7   -43.320003 -44.453710
## 8    25.363102  22.737596
## 9   -20.648195 -22.985108
## 10   -2.322031  -3.942871

De la salida anterior vemos que los \(\tilde{b}_0\) son cercanos a los valores reales de \(b_0\).

La comparación anterior solo es posible cuando usamos datos simulados. Cuando se usan datos de un fenómeno real no se tienen los valores de \(b_0\).

5.2 Función `glmer`

La función glmer es la función del paquete lme4 para ajustar modelo lineales generalizados mixtos y su estructura es la siguiente:

glmer(formula, data = NULL, family = gaussian, control = glmerControl(), 
      start = NULL, verbose = 0L, nAGQ = 1L, subset, weights, na.action, 
      offset, contrasts = NULL, mustart, etastart, devFunOnly = FALSE)

Los principales argumentos de la función son:

formula: es una fórmula similar a la usada en el modelo lineal clásico. Un ejemplo de fórmula sería y ~ 1 + x1 + x2 + (1 + x2) con la cual se indican los efectos fijos y los efectos aleatorios del modelo. Más abajo hay una tabla con más detalles sobre la fórmula.
data: marco de datos donde están las variables.
family: argumento para seleccionar la distribución de la variable respuesta. Para más detalles de las distribuciones y funciones de enlace se recomienda ver la ayuda de la función family escribiendo esto en la consola ?family.

Ejemplo: modelo gamma con intercepto aleatorio

En este ejemplo vamos a simular observaciones \(n_i=20\) observaciones para \(G=10\) grupos (en total 200 obs) que tengan la estructura mostrada abajo. El objetivo del ejemplo es ilustrar el uso de la función glmer para estimar los parámetros del modelo.

\[\begin{align*} y_{ij} | b_0 &\sim Gamma(\mu_{ij}, \phi) \\ \log(\mu_{ij}) &= 2 - 8 x_{ij} + b_{0i} \\ \phi &= 0.5 \\ b_{0} &\sim N(0, \sigma^2_{b0}=9) \\ x_{ij} &\sim U(0, 1) \end{align*}\]

El vector de parámetros de este modelo es \(\boldsymbol{\Theta}=(\beta_0=2, \beta_1=-8, \phi=0.5, \sigma_{b0}=3)^\top\).

Solución.

La función rgamma_glm que se muestra a continuación es una modificación de la función rgamma para tener la parametrización usada en los glm.

rgamma_glm <- function(n, mu, phi) {
  x <- rgamma(n=n, shape=1/phi, scale=mu*phi)
  return(x)
}

A continuación el código para simular datos del modelo de interés.

ni <- 20
G <- 10
nobs <- ni * G
grupo <- factor(rep(x=1:G, each=ni))
obs <- rep(x=1:ni, times=G)
set.seed(123456)
x <- runif(n=nobs, min=0, max=1)
set.seed(123456)
b0 <- rnorm(n=G, mean=0, sd=sqrt(9)) # Intercepto aleatorio
b0 <- rep(x=b0, each=ni)             # El mismo intercepto aleatorio pero repetido
media <- exp(2 - 8 * x + b0)
set.seed(123456)
y <- rgamma_glm(n=nobs, mu=media, phi=0.5)
datos <- data.frame(obs, grupo, b0, x, media, y)

Vamos a explorar los datos simulados.

library(rmarkdown)
paged_table(datos, options = list(rows.print = 6, cols.print=6))

El siguiente paso es explorar los datos simulados. El código para dibujar los datos se muestra abajo.

library(ggplot2)
ggplot(datos, aes(x, y, color=grupo) ) + 
  geom_point() + 
  labs(colour="Grupo/Cluster")

Para estimar los parámetros del modelos se usa la función glmer de la siguiente forma.

library(lme4)
fit2 <- glmer(y ~ x + (1 | grupo), family=Gamma(link="log"), data = datos)

La función summary se puede usar sobre el objeto fit2 para obtener una tabla de resumen, a continuación se muestra la salida de summary.

summary(fit2)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: y ~ x + (1 | grupo)
##    Data: datos
## 
##      AIC      BIC   logLik deviance df.resid 
##    591.6    604.8   -291.8    583.6      196 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.4164 -0.7151 -0.2155  0.5631  3.0609 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  grupo    (Intercept) 4.3265   2.0800  
##  Residual             0.4728   0.6876  
## Number of obs: 200, groups:  grupo, 10
## 
## Fixed effects:
##             Estimate Std. Error  t value Pr(>|z|)    
## (Intercept)   4.3965     0.9566    4.596 4.31e-06 ***
## x            -8.0784     0.0411 -196.549  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##   (Intr)
## x -0.021

Según el resultado anterior \(\hat{\boldsymbol{\Theta}}=(\hat{\beta}_0=4.3965, \hat{\beta}_1=-8.0784, \hat{\phi}=0.4728, \hat{\sigma}_{bo}=4.3265)^\top\) mientras que el vector real de parámetros es \(\boldsymbol{\Theta}=(\beta_0=2, \beta_1=-8, \phi=0.50, \sigma_{b0}=3)^\top\).

El valor de \(\hat{\phi}\) se encuentra en la línea Residual y columna Variance del bloque Random effects que está en el summary.

Ejemplo: modelo inversa gaussiana con intercepto y pendiente aleatoria

\[\begin{align*} y_{ij} | b_0, b_1 &\sim InvGauss(\mu_{ij}, \phi) \\ \log(\mu_{ij}) &= 2 - 8 x_{ij} + b_{0i} + b_{1i} x_{ij} \\ \phi &= 0.5 \\ \left ( \begin{matrix} b_{0} \\ b_{1} \end{matrix} \right ) &\sim N\left ( \left [ \begin{matrix} 0 \\ 0 \end{matrix} \right ], \left [ \begin{matrix} \sigma^2_{b0}=1 & \sigma_{b01}=0.5 \\ \sigma_{b01} & \sigma^2_{b1}=1 \end{matrix} \right ] \right ) \\ x_{ij} &\sim U(0, 1) \end{align*}\]

El vector de parámetros de este modelo es \(\boldsymbol{\Theta}=(\beta_0=2, \beta_1=-8, \phi=0.5, \sigma_{b0}^2=1, \sigma_{b1}^2=1, \sigma_{b01}=0.5)^\top\).

Solución.

Para simular valores de una distribución inversa gaussiana vamos a usar la función rinvgauss del paquete statmod de (R-statmod?).

A continuación el código para simular datos del modelo de interés.

ni <- 20
G <- 10
nobs <- ni * G
grupo <- factor(rep(x=1:G, each=ni))
obs <- rep(x=1:ni, times=G)
set.seed(12345)
x <- runif(n=nobs, min=0, max=1)
set.seed(12345)

Sigma <- matrix(c(1, 0.5, # Matriz de var-cov
                  0.5, 1), ncol=2, nrow=2)

b <- MASS::mvrnorm(n=G, mu=rep(0, 2), Sigma=Sigma)   # Simulando b0 y b1
b <- apply(b, MARGIN=2, function(c) rep(c, each=ni)) # Replicando
b0 <- as.vector(b[, 1]) # Separando los b0
b1 <- as.vector(b[, 2]) # Separando los b1

media <- exp(2 - 8 * x + b0 + b1 * x)
phi <- 0.5

set.seed(12345)
y <- statmod::rinvgauss(n=nobs, mean=media, dispersion=phi)
datos <- data.frame(obs, grupo, b0, x, media, y)

Vamos a explorar los datos simulados.

library(rmarkdown)
paged_table(datos, options = list(rows.print = 6, cols.print=6))

El siguiente paso es explorar los datos simulados. El código para dibujar los datos se muestra abajo.

library(ggplot2)
ggplot(datos, aes(x, y, color=grupo) ) + 
  geom_point() + 
  labs(colour="Grupo/Cluster")

Para estimar los parámetros del modelos se usa la función glmer de la siguiente forma.

library(lme4)
fit3 <- glmer(y ~ x + (1 + x | grupo), 
              family=inverse.gaussian(link="log"), data = datos)

La función summary se puede usar sobre el objeto fit3 para obtener una tabla de resumen, a continuación se muestra la salida de summary.

summary(fit3)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: inverse.gaussian  ( log )
## Formula: y ~ x + (1 + x | grupo)
##    Data: datos
## 
##      AIC      BIC   logLik deviance df.resid 
##  -1060.0  -1040.2    536.0  -1072.0      194 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7657 -0.5806 -0.1294  0.4044  6.5944 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  grupo    (Intercept) 0.6118   0.7822       
##           x           0.3740   0.6116   0.33
##  Residual             0.7266   0.8524       
## Number of obs: 200, groups:  grupo, 10
## 
## Fixed effects:
##             Estimate Std. Error t value Pr(>|z|)    
## (Intercept)   1.7907     0.2947   6.077 1.22e-09 ***
## x            -8.0237     0.2345 -34.218  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##   (Intr)
## x 0.268

Según el resultado anterior \(\hat{\boldsymbol{\Theta}}=(\hat{\beta}_0=1.7907, \hat{\beta}_1=-8.0237, \hat{\phi}=0.7266, \hat{\sigma}_{b0}^2=0.6118, \hat{\sigma}_{b1}^2=0.3740, \hat{\sigma}_{b01}=0.1579)^\top\) mientras que el vector real de parámetros es \(\boldsymbol{\Theta}=(\beta_0=2, \beta_1=-8, \phi=0.5, \sigma_{b0}^2=1, \sigma_{b1}^2=1, \sigma_{b01}=0.5)^\top\).

El parámetro estimado \(\hat{\sigma}_{b01}\) se obtiene utilizando la ecuación de correlación (\(\rho\)) que relaciona la covarianza y desviaciones de los efectos aleatorios: \(\rho_{b0b1} = \sigma_{b0b1}/(\sigma_{b0} \times \sigma_{b1})\).

El valor de \(\hat{\phi}\) se encuentra en la línea Residual y columna Variance del bloque Random effects que está en el summary.

References

Bates, Douglas, Martin Maechler, Ben Bolker, and Steven Walker. 2023. Lme4: Linear Mixed-Effects Models Using Eigen and S4. https://github.com/lme4/lme4/.

5 Paquete lme4

5.1 Función lmer

Ejemplo: modelo normal con intercepto aleatorio

Ejemplo: recuperando los interceptos aleatorios

5.2 Función glmer

Ejemplo: modelo gamma con intercepto aleatorio

Ejemplo: modelo inversa gaussiana con intercepto y pendiente aleatoria

References

5.1 Función `lmer`

5.2 Función `glmer`