El presente documento realizará una exposición del documento Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1), 1–48. https://doi.org/10.18637/jss.v067.i01

A lo largo de la discusión sobre lme4, trabajaremos con un conjunto de datos sobre el tiempo promedio de reacción por día para sujetos en un estudio de privación de sueño. En el día 0, los sujetos tuvieron su cantidad normal de sueño. A partir de esa noche, se les restringió a 3 horas de sueño por noche. La variable de respuesta, “Reacción”, representa los tiempos de reacción promedio en milisegundos (ms) en una serie de pruebas realizadas cada día a cada sujeto.

Hemos elegido el ejemplo de sleepstudy porque es un ejemplo relativamente pequeño y simple para ilustrar la teoría y la práctica subyacente en lmer.

library(lme4)
## Loading required package: Matrix
str(sleepstudy)
## 'data.frame':    180 obs. of  3 variables:
##  $ Reaction: num  250 259 251 321 357 ...
##  $ Days    : num  0 1 2 3 4 5 6 7 8 9 ...
##  $ Subject : Factor w/ 18 levels "308","309","310",..: 1 1 1 1 1 1 1 1 1 1 ...
head(sleepstudy)
##   Reaction Days Subject
## 1 249.5600    0     308
## 2 258.7047    1     308
## 3 250.8006    2     308
## 4 321.4398    3     308
## 5 356.8519    4     308
## 6 414.6901    5     308

El tiempo de reacción de cada sujeto aumenta aproximadamente de forma lineal con el número de días de privación de sueño. Sin embargo, parece que los sujetos también varían en las pendientes e intercepciones de estas relaciones, lo que sugiere un modelo con pendientes e intercepciones aleatorias. Como veremos, dicho modelo se puede ajustar minimizando el criterio REML

fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
summary(fm1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Reaction ~ Days + (Days | Subject)
##    Data: sleepstudy
## 
## REML criterion at convergence: 1743.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9536 -0.4634  0.0231  0.4634  5.1793 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  Subject  (Intercept) 612.10   24.741       
##           Days         35.07    5.922   0.07
##  Residual             654.94   25.592       
## Number of obs: 180, groups:  Subject, 18
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  251.405      6.825  36.838
## Days          10.467      1.546   6.771
## 
## Correlation of Fixed Effects:
##      (Intr)
## Days -0.138

Las estimaciones de las desviaciones estándar de los efectos aleatorios para la intercepción y la pendiente son de 24,74 ms y 5,92 ms/día. Los coeficientes de efectos fijos, β, son de 251,4 ms y 10,47 ms/día para la intercepción y la pendiente. En este modelo, una interpretación de estos efectos fijos es que son los valores medios estimados de la población para la intercepción y la pendiente aleatoria.

La función lmer se compone de cuatro módulos en gran parte independientes. En el primer módulo, se analiza y convierte una fórmula de modelo mixto en los datos necesarios para especificar un modelo mixto lineal (Sección 2). El segundo módulo utiliza estos datos para construir una función en R que toma los parámetros de covarianza, θ, como argumentos y devuelve el doble negativo del logaritmo del perfil de la verosimilitud o el criterio REML (Sección 3). El tercer módulo optimiza esta función objetivo para producir estimaciones de máxima verosimilitud (ML) o REML de θ (Sección 4). Finalmente, el cuarto módulo proporciona utilidades para interpretar el modelo optimizado (Sección 5).

Para ilustrar esta modularidad, recreamos el objeto fm1 en una serie de cuatro pasos modulares.

  1. módulo de fórmula.
parsedFormula <- lFormula(formula = Reaction ~ Days + (Days | Subject),
 data = sleepstudy)
  1. objetivo de la función
 devianceFunction <- do.call(mkLmerDevfun, parsedFormula)
  1. obtimizacion del modelo
optimizerOutput <- optimizeLmer(devianceFunction)
  1. módulo de salida
mkMerMod(rho = environment(devianceFunction), opt = optimizerOutput,
 reTrms = parsedFormula$reTrms, fr = parsedFormula$fr)
## Linear mixed model fit by REML ['lmerMod']
## REML criterion at convergence: 1743.628
## Random effects:
##  Groups   Name        Std.Dev. Corr
##  Subject  (Intercept) 24.741       
##           Days         5.922   0.07
##  Residual             25.592       
## Number of obs: 180, groups:  Subject, 18
## Fixed Effects:
## (Intercept)         Days  
##      251.41        10.47

Seccion 2: Formula modelo mixto.

Como la mayoría de las funciones de ajuste de modelos en R, lmer toma como sus dos primeros argumentos una fórmula que especifica el modelo y los datos con los que evaluar la fórmula. Este segundo argumento, data, es opcional pero recomendado y suele ser el nombre de un marco de datos de R. En la función lm de R para ajustar modelos lineales, las fórmulas tienen la forma resp ~ expr, donde resp determina la variable de respuesta y expr es una expresión que especifica las columnas de la matriz del modelo. Las fórmulas para la función lmer contienen términos especiales de efectos aleatorios.

resp ~ FEexpr + (REexpr1 | factor1) + (REexpr2 | factor2) + ...

Donde FEexpr es una expresión que determina las columnas de la matriz del modelo de efectos fijos, X, y los términos de efectos aleatorios, (REexpr1 | factor1) y (REexpr2 | factor2), determinan tanto la matriz del modelo de efectos aleatorios, Z, como la estructura del factor de covarianza relativa, Λθ. En principio, una fórmula de modelo mixto puede contener un número arbitrario de términos de efectos aleatorios, pero en la práctica, el número de dichos términos suele ser bajo.

Antes de describir los detalles sobre cómo lme4 analiza las fórmulas de modelos mixtos, proporcionamos una explicación informal y luego algunos ejemplos. Nuestra discusión asume familiaridad con el paradigma de modelado estándar de R.

Cada término de efectos aleatorios tiene la forma (expr | factor). La expresión expr se evalúa como una fórmula de modelo lineal, lo que produce una matriz de modelo siguiendo las mismas reglas utilizadas en las funciones de modelado estándar de R (por ejemplo, lm o glm). La expresión factor se evalúa como un factor en R.

Una forma de entender el operador de barra vertical es como un tipo especial de interacción entre la matriz de modelo y el factor de agrupamiento. Esta interacción garantiza que las columnas de la matriz de modelo tengan efectos diferentes para cada nivel del factor de agrupamiento. Lo que hace que esto sea un tipo especial de interacción es que estos efectos se modelan como variables aleatorias no observadas, en lugar de parámetros fijos desconocidos. Se ha escrito mucho sobre las diferencias prácticas y filosóficas importantes entre estos dos tipos de interacciones. Por ejemplo, la implementación de efectos aleatorios de estas interacciones se puede utilizar para obtener estimaciones de contracción de los coeficientes de regresión o para tener en cuenta la falta de independencia en los residuos debido a la estructura de bloque o a mediciones repetidas (Para más información, consultar http://dx.doi.org/10.1214/009053604000001048)

La Tabla 1 proporciona varios ejemplos de los lados derechos de las fórmulas de modelos mixtos. El primer ejemplo, (1 | g), es la fórmula de modelo mixto más simple, donde cada nivel del factor de agrupación, g, tiene su propia intersección aleatoria. La media y la desviación estándar de estas intersecciones son parámetros a estimar. Nuestra descripción de este modelo incorpora cualquier media distinta de cero de los efectos aleatorios como parámetros de efectos fijos. Si se desea especificar que una intersección aleatoria tiene medias conocidas a priori, se puede utilizar la función offset como en el segundo modelo en la Tabla 1. Este modelo no contiene efectos fijos, o más precisamente, la matriz de modelo de efectos fijos, X, tiene columnas cero y β tiene longitud cero.

Tabla 1. Ejemplos de los lados derechos de las fórmulas de modelos de efectos mixtos. Los nombres de los factores de agrupación se denominan g, g1 y g2, y las covariables y desplazamientos conocidos a priori como x y o.
Formula Alternativa Descripción
(1 | g) 1 + (1 | g) Intercepto aleatorio con media fija
0 + offset(o) + (1 | g) -1 + offset(o) + (1 | g) Intercepto aleatorio con medias a priori.
(1 | g1/g2) (1 | g1) + (1 | g1:g2) Intercepto variable entre g1 y g2 dentro de g1.
(1 | g1) + (1 | g2) 1 + (1 | g1) + (1 | g2) Intercepto variable entre g1 y g2.
x + (x | g) 1 + x + (1 + x | g) Intercepto y pendiente correlacionados.
x + (x || g) 1 + x + (1 | g) + (0 + x | g) Intercepto y pendiente no correlacionados.

También podemos construir modelos con múltiples factores de agrupamiento. Por ejemplo, si las observaciones están agrupadas por g2, que está anidado dentro de g1, entonces la tercera fórmula en la Tabla 1 puede usarse para modelar la variación en la intercepción. Un objetivo común en el modelado mixto es tener en cuenta esa estructura anidada (o jerárquica). Sin embargo, uno de los aspectos más útiles de lme4 es que puede usarse para ajustar efectos aleatorios asociados con factores de agrupamiento no anidados. Por ejemplo, supongamos que los datos están agrupados cruzando completamente dos factores, g1 y g2; en ese caso, se puede emplear la cuarta fórmula en la Tabla 1. Estos modelos son comunes en la teoría de respuesta a los ítems, donde los factores de sujeto y de ítem están completamente cruzados. Además de las intercepciones variables, también podemos tener pendientes variables (por ejemplo, los datos de sleepstudy). El quinto ejemplo en la Tabla 1 proporciona un modelo donde tanto la intercepción como la pendiente varían entre los niveles del factor de agrupamiento.

Efectos aleatorios no correlacionados

Por defecto, lme4 asume que todos los coeficientes asociados al mismo término de efectos aleatorios están correlacionados. Para especificar una pendiente e intercepción no correlacionadas (por ejemplo), se puede utilizar la notación de doble barra, (x || g), o, de manera equivalente, utilizar múltiples términos de efectos aleatorios, x + (1 | g) + (0 + x | g), como se muestra en el último ejemplo de la Tabla 1. Por ejemplo, si se examinan los resultados del modelo fm1 de los datos de sleepstudy utilizando summary(fm1), se observará que la correlación estimada entre la pendiente para los Días y la intercepción es bastante baja (0.066) (más adelante para obtener más información sobre cómo extraer la matriz de covarianza de efectos aleatorios). Podemos usar la notación de doble barra para ajustar un modelo que excluya un parámetro de correlación:

fm2 <- lmer(Reaction ~ Days + (Days || Subject), sleepstudy)

Aunque los modelos mixtos donde se asume que las pendientes y las intercepciones aleatorias son independientes se utilizan comúnmente para reducir la complejidad de los modelos con pendientes aleatorias, tienen una leve desventaja sutil. Los modelos en los que se permite que las pendientes y las intercepciones tengan una correlación distinta de cero (por ejemplo, fm1) son invariantes a cambios aditivos en el predictor continuo (en este caso, los Días). Esta invarianza se descompone cuando se restringe la correlación a cero; cualquier cambio en el predictor conducirá necesariamente a un cambio en la correlación estimada, y en la verosimilitud y predicciones del modelo.

Por ejemplo, podemos eliminar la correlación en fm1 simplemente agregando una cantidad igual a la proporción de las desviaciones estándar entre sujetos estimadas multiplicadas por la correlación estimada (es decir, σslope/σintercept · ρslope:intercept) a la variable Días. El uso de modelos como fm2 debería idealmente restringirse a casos en los que el predictor se mida en una escala de razón (es decir, el punto cero en la escala tiene un significado, no solo es una ubicación definida por conveniencia o convención), como es el caso aquí.

LA implementación se llama lme4pureR y está disponible actualmente en Github https://github.com/lme4/lme4pureR/

La interfaz de usuario para el módulo de salida consiste en gran medida en un conjunto de métodos (Tabla 8) para objetos de clase ‘merMod’, que es la clase de objetos devueltos por lmer. Además de estos métodos, la función getME se puede utilizar para extraer varios objetos de un modelo mixto ajustado en lme4.

graficos

lme4 proporciona herramientas para generar la mayoría de las representaciones gráficas diagnósticas estándar (a excepción de aquellas que incorporan medidas de influencia, por ejemplo, gráficos de distancia de Cook y gráficos de palanca), siguiendo el modelo de los gráficos diagnósticos de nlme. En particular, el método de gráficos ‘plot’ en R base para modelos lineales (objetos de la clase lm) puede ser utilizado para crear los gráficos estándar de ajuste versus residuos

plot(fm1, type = c("p", "smooth"))

plot(fm1, sqrt(abs(resid(.))) ~ fitted(.), type = c("p", "smooth"))

En contraste con el método de gráfico para objetos ‘lm’, estos gráficos de escala-ubicación (scale-location) y Q-Q están basados en residuos crudos en lugar de residuos estandarizados.

Además de estos gráficos diagnósticos estándar, que examinan la validez de varias suposiciones (linealidad, homocedasticidad, normalidad) a nivel de los residuos, también se pueden utilizar los métodos dotplot y qqmath para los modos condicionales (es decir, los objetos ranef.mer generados por ranef(fit)) para verificar patrones interesantes y la normalidad en los modos condicionales. lme4 no proporciona diagnósticos de influencia, pero estos (y otros procedimientos diagnósticos útiles) están disponibles en los paquetes dependientes HLMdiag e influence.ME (Para más información, consultar https://doi.org/10.1017/CBO9780511790942)

Finalmente, la simulación predictiva posterior es una herramienta diagnóstica generalmente útil, adaptada de métodos bayesianos, para explorar la idoneidad del modelo. Los usuarios seleccionan alguna estadística resumida de interés. Luego calculan la estadística resumida para un conjunto de simulaciones y observan dónde se sitúan los datos observados dentro de la distribución simulada; si los datos observados son extremos, podríamos concluir que el modelo es una representación deficiente de la realidad. Por ejemplo, utilizando el ajuste del estudio del sueño y eligiendo el rango intercuartílico de los tiempos de reacción como la estadística resumida

iqrvec <- sapply(simulate(fm1, 1000), IQR)
obsval <- IQR(sleepstudy$Reaction)
post.pred.p <- mean(obsval >= c(obsval, iqrvec))

El valor p predictivo posterior (unilateral) de 0.795 indica que el modelo representa los datos adecuadamente, al menos para esta estadística resumida. A diferencia del caso bayesiano completo, el procedimiento descrito aquí no tiene en cuenta la incertidumbre en los parámetros estimados. No obstante, debería ser una aproximación razonable cuando la variación residual es grande.

Para ilustrar el método de actualización para objetos merMod, construimos un modelo solo con interceptos aleatorios. Esto se puede realizar de diversas formas, pero optamos por eliminar primero el término de efectos aleatorios (Days | Subject) y luego agregar un nuevo término con un intercepto aleatorio.

fm3 <- update(fm1, . ~ . - (Days | Subject) + (1 | Subject))
formula(fm3)
## Reaction ~ Days + (1 | Subject)

El método de resumen para objetos merMod ofrece información sobre el ajuste del modelo. Aquí consideramos la salida de summary(fm3) en cuatro pasos. Las primeras líneas de salida indican que el modelo fue ajustado por REML, así como el valor del criterio REML en la convergencia (que también se puede extraer usando REMLcrit(fm3)).

summary(fm3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Reaction ~ Days + (1 | Subject)
##    Data: sleepstudy
## 
## REML criterion at convergence: 1786.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2257 -0.5529  0.0109  0.5188  4.2506 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Subject  (Intercept) 1378.2   37.12   
##  Residual              960.5   30.99   
## Number of obs: 180, groups:  Subject, 18
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 251.4051     9.7467   25.79
## Days         10.4673     0.8042   13.02
## 
## Correlation of Fixed Effects:
##      (Intr)
## Days -0.371
REMLcrit(fm3)
## [1] 1786.465

Descomposición secuencial y comparación de modelos

Los objetos de la clase merMod tienen un método anova que devuelve estadísticas F correspondientes a la descomposición secuencial de las contribuciones de los términos de efectos fijos. Para ilustrar esta descomposición secuencial ANOVA, ajustamos un modelo con términos lineales y cuadráticos ortogonales para Days.

fm4 <- lmer(Reaction ~ polyDays[ , 1] + polyDays[ , 2] +
 (polyDays[ , 1] + polyDays[ , 2] | Subject),
 within(sleepstudy, polyDays <- poly(Days, 2)))

Las magnitudes relativas de las dos sumas de cuadrados indican que el término cuadrático explica mucha menos variación que el término lineal. Además, las magnitudes de las dos estadísticas F sugieren fuertemente la significancia del término lineal, pero no del término cuadrático. Es importante señalar que esta es solo una evaluación informal de la significancia ya que no hay valores p asociados con estas estadísticas F, un tema que se discute con más detalle en la siguiente subsección.

anova(fm1, fm2, fm3)
## refitting model(s) with ML (instead of REML)
## Data: sleepstudy
## Models:
## fm3: Reaction ~ Days + (1 | Subject)
## fm2: Reaction ~ Days + ((1 | Subject) + (0 + Days | Subject))
## fm1: Reaction ~ Days + (Days | Subject)
##     npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)    
## fm3    4 1802.1 1814.8 -897.04   1794.1                          
## fm2    5 1762.0 1778.0 -876.00   1752.0 42.0754  1  8.782e-11 ***
## fm1    6 1763.9 1783.1 -875.97   1751.9  0.0639  1     0.8004    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El resultado muestra estadísticas χ2 que representan la diferencia en desviación entre modelos sucesivos, así como valores p basados en comparaciones de pruebas de razón de verosimilitud. En este caso, la comparación secuencial muestra que agregar un efecto lineal de tiempo no correlacionado con la intersección conduce a una caída enorme y significativa en la desviación (∆desviación ≈ 42, p ≈ 10−10), mientras que la adición adicional de correlación entre la pendiente e intersección conduce a un cambio trivial y no significativo en la desviación (∆desviación ≈ 0.06). Para objetos de clase lmerMod, el comportamiento predeterminado es volver a ajustar los modelos con ML si se ajustan con REML = TRUE, lo cual es necesario para obtener respuestas sensatas al comparar modelos que difieren en sus efectos fijos; esto se puede controlar a través del argumento refit.

Una de las decisiones de diseño más controvertidas de lme4 ha sido omitir la presentación de valores p asociados con descomposiciones secuenciales de ANOVA de efectos fijos. La ausencia de resultados analíticos para las distribuciones nulas de estimaciones de parámetros en situaciones complejas (por ejemplo, diseños desequilibrados o parcialmente cruzados) es un problema de larga data en la inferencia de modelos mixtos. Aunque las distribuciones nulas (y las distribuciones muestrales de estimaciones no nulas) son asintóticamente normales, estas distribuciones no siguen una distribución t para muestras de tamaño finito, ni tampoco lo hacen las distribuciones nulas correspondientes de las diferencias en desviaciones escaladas F. Por lo tanto, los métodos aproximados para calcular los grados de libertad aproximados para las distribuciones t, o los grados de libertad del denominador para las estadísticas F, son a lo sumo soluciones ad hoc.