Modelos mixtos con R
Christian Castro
2024-10-22
1 Introducción
Los modelos mixtos son útiles cuando la estructura de datos que estudiamos es jerárquica o posee algún tipo de agrupación. Ejemplos de datos agrupados son, por ejemplo, estudiantes dentro de escuelas, pacientes dentro de hospitales o mediciones repetidas de los mismos individuos.
Los modelos mixtos combinan efectos fijos y efectos aleatorios. Los primeros son aquellos constantes para todos los individuos del set de datos, mientras los segundos son los efectos que varían entre los grupos.
Los modelos mixtos ofrecen varias ventajas sobre simplemente separar los datos en grupos:
Los modelos mixtos permiten modelar tanto la variabilidad dentro de los grupos (intra) como entre los grupos (inter).
Al utilizar toda la información disponible y modelar las correlaciones dentro de los grupos, los modelos mixtos pueden proporcionar estimaciones más precisas y eficientes que los análisis separados por grupos.
Los modelos mixtos permiten incluir efectos aleatorios, que pueden capturar la variabilidad no explicada por los efectos fijos. Esto es útil cuando hay factores no medidos que podrían influir en las observaciones dentro de los grupos.
Separar los datos en grupos y analizarlos por separado puede llevar a conclusiones específicas para cada grupo, pero los modelos mixtos permiten hacer inferencias más generales que se aplican a la población completa, considerando la variabilidad entre los grupos.
Los modelos mixtos son robustos frente a datos desbalanceados, donde el número de observaciones puede variar significativamente entre los grupos. Esto puede ser problemático en análisis separados, pero los modelos mixtos pueden manejar esta variabilidad de manera más efectiva.
Los modelos mixtos permiten modelar interacciones complejas entre efectos fijos y aleatorios, proporcionando una comprensión más profunda de cómo diferentes factores influyen en la variable de respuesta.
1.1 Primer ejemplo.
Consideremos un conjunto de datos que proporciona una simulación de cómo las horas de estudio y los efectos específicos de cada escuela pueden influir en las puntuaciones de los estudiantes en una prueba.
Se consideran 1000 estudiantes agrupados en 10 escuelas.
school: Representa las diferentes escuelas en el estudio.
hours: La cantidad de horas que cada estudiante ha dedicado al estudio. Estas horas se generan aleatoriamente siguiendo una distribución normal con una media de 5 horas y una desviación estándar de 2 horas.
score: La puntuación obtenida por cada estudiante en una prueba. Esta puntuación se calcula en función de varios factores:
- Un valor base de 50 puntos.
- Un incremento de 10 puntos por cada hora de estudio.
- Un efecto aleatorio específico de cada escuela, que puede aumentar o disminuir la puntuación.
- Un término de error aleatorio para añadir variabilidad adicional.1.1.1 Apliquemos un modelo lineal simple.
##
## Call:
## lm(formula = score ~ hours, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.6854 -5.5945 -0.1313 5.5825 24.7978
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 50.7563 0.6927 73.28 <2e-16 ***
## hours 9.8136 0.1281 76.62 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.026 on 998 degrees of freedom
## Multiple R-squared: 0.8547, Adjusted R-squared: 0.8546
## F-statistic: 5871 on 1 and 998 DF, p-value: < 2.2e-16
1.1.2 Apliquemos un modelo mixto.
# Modelo lineal mixto
model_lmm <- lmer(score ~ hours + (1 | school), data = data)
# Resumen del modelo lineal mixto
summary(model_lmm)## Linear mixed model fit by REML ['lmerMod']
## Formula: score ~ hours + (1 | school)
## Data: data
##
## REML criterion at convergence: 6125.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0352 -0.6848 0.0118 0.7021 3.3473
##
## Random effects:
## Groups Name Variance Std.Dev.
## school (Intercept) 43.16 6.57
## Residual 25.50 5.05
## Number of obs: 1000, groups: school, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 50.77737 2.12307 23.92
## hours 9.80945 0.08094 121.19
##
## Correlation of Fixed Effects:
## (Intr)
## hours -0.192REML (Restricted Maximum Likelihood): Método utilizado para estimar los parámetros del modelo mixto.
Formula score ~ 1 + (1 | school) indica que estamos modelando score con un intercepto fijo y un intercepto aleatorio para school.
Recordemos que los efectos fijos son aquellos que son constantes para todos los individuos en el estudio, mientras que los efectos aleatorios varían entre los grupos o individuos.
1.1.3 Obtengamos los efectos aleatorios por escuela:
## school intercept
## 1 1 -4.23430052
## 2 2 -4.94361505
## 3 3 -0.08656364
## 4 4 -0.65645784
## 5 5 -11.82852535
## 6 6 5.54007363
## 7 7 1.75963135
## 8 8 12.28860502
## 9 9 3.93576307
## 10 10 -1.77461066
Supongamos que el intercepto fijo es 50. Si una escuela tiene un intercepto aleatorio de +5, esto significa que los estudiantes de esa escuela, en promedio, tienen un puntaje de 55. Si otra escuela tiene un intercepto aleatorio de -3, los estudiantes de esa escuela, en promedio, tienen un puntaje de 47.
2 Modelos mixtos
Los modelos de efectos mixtos, o simplemente mixtos, generalmente se refieren a una mezcla de efectos fijos y aleatorios.
El término “efectos fijos” es quizás un término pobre, pero no por ello menos adecuado, para los efectos principales típicos que se pueden observar en un modelo de regresión lineal, es decir, la parte no aleatoria de un modelo mixto. En algunos contextos, se los denomina “efecto promedio de la población”. Los efectos aleatorios son simplemente aquellos específicos de una unidad de observación, cualquiera sea su definición.
En términos de estructura, los datos pueden tener una o varias fuentes de agrupamiento, y ese agrupamiento puede ser jerárquico, de modo que los agrupamientos estén anidados dentro de otros agrupamientos. Un ejemplo serían las pruebas de aptitud escolar administradas varias veces a los estudiantes (observaciones repetidas anidadas dentro de los estudiantes, estudiantes anidados dentro de las escuelas, escuelas anidadas dentro de los distritos). En otros casos, no hay una estructura de anidamiento. Un ejemplo sería un experimento de tiempo de reacción donde los participantes realizan el mismo conjunto de tareas. Si bien las observaciones están anidadas dentro de un individuo, las observaciones también están agrupadas según el tipo de tarea. Algunos usan los términos anidado y cruzado para distinguir entre estos escenarios. Además, el agrupamiento puede ser equilibrado o no. Podríamos esperar más equilibrio en estudios de naturaleza experimental, pero definitivamente no en otros casos, por ejemplo, cuando el agrupamiento es algo así como una unidad geográfica y las observaciones son personas.
En lo que sigue, veremos modelos de efectos mixtos en todas estas situaciones de datos. En general, nuestro enfoque será el mismo, ya que dicha agrupación es en realidad más una propiedad de los datos que del modelo. Sin embargo, es importante tener una idea de la flexibilidad de los modelos mixtos para manejar una variedad de situaciones de datos. También vale la pena señalar que pueden surgir otras estructuras, como temporales o espaciales. Las mencionaremos en algunos lugares de este documento, pero no serán el foco.
Existen varios tipos de modelos mixtos que se pueden utilizar dependiendo de la naturaleza de los datos y el objetivo del análisis. Aquí tienes una lista de los tipos más comunes:
- Modelos Lineales Mixtos (LMM) Estos modelos combinan efectos fijos y aleatorios y son adecuados para datos continuos.
modelo_lmm <- lmer(y ~ x1 + x2 + (1 | grupo), data = datos)
- Modelos Lineales Generalizados Mixtos (GLMM) Extienden los LMM para datos que no siguen una distribución normal, como datos binarios o de conteo.
Ejemplo para datos binarios:
modelo_glmm_bin <- glmer(y_bin ~ x1 + x2 + (1 | grupo), data = datos, family = binomial)
Ejemplo para datos de conteo:
modelo_glmm_pois <- glmer(y_count ~ x1 + x2 + (1 | grupo), data = datos, family = poisson)
- Modelos Lineales Mixtos Multinivel Permiten incluir múltiples niveles de efectos aleatorios, útiles para datos anidados.
Ejemplo:
modelo_multinivel <- lmer(y ~ x1 + x2 + (1 | escuela/clase), data = datos)
- Modelos de Supervivencia Mixtos Para datos de tiempo hasta el evento, aunque no se manejan directamente con lme4, se pueden usar paquetes como coxme.
Ejemplo:
library(coxme)
modelo_coxme <- coxme(Surv(tiempo, evento) ~ x1 + x2 + (1 | grupo), data = datos)
- Modelos No Lineales Mixtos Para ajustar modelos no lineales con efectos aleatorios, se puede usar el paquete nlme.
Ejemplo:
library(nlme)
modelo_nlme <- nlme(y ~ f(x1, x2), data = datos, fixed = x1 + x2 ~ 1, random = x1 ~ 1 | grupo)
- Modelos Mixtos Multinomial Para datos categóricos con más de dos categorías, se pueden usar modelos mixtos multinomiales.
Ejemplo:
library(mme)
# Ajuste de un modelo multinomial mixto
modelo_multinomial <- mme::mme(y ~ x1 + x2, random = ~ 1 | grupo, data = datos)
- Modelos Lineales Mixtos (LMM)
2.1 1.1 Modelos de Efectos Cruzados
modelo_cruzado <- lmer(y ~ x1 + x2 + (1 | estudiante) + (1 | profesor), data = datos)2.2 1.2 Modelos de Intercepto y Pendiente Aleatorios
modelo_intercepto_pendiente <- lmer(y ~ x1 + x2 + (1 + x1 + x2 | grupo), data = datos)
2.3 1.3 Modelos de Pendiente Aleatoria
modelo_pendiente <- lmer(y ~ x1 + x2 + (1 + x1 | grupo), data = datos)
2.4 Modelo de intersecciones aleatorias.
En el contexto de los Modelos Lineales Mixtos (LMMs), los modelos de intersecciones aleatorias se utilizan para manejar la variabilidad entre grupos. Aquí hay algunos puntos clave sobre estos modelos:
Intersecciones Aleatorias: Se utilizan para capturar la variabilidad entre diferentes niveles de un factor aleatorio. Por ejemplo, en un estudio con múltiples bloques o regiones, cada bloque puede tener su propia intersección aleatoria.
Estructura Jerárquica: Los modelos pueden tener efectos aleatorios anidados (e.g., plantas dentro de bloques) o cruzados (e.g., plantas medidas en diferentes condiciones).
Estimación: Se utilizan métodos como la Máxima Verosimilitud (ML) y la Máxima Verosimilitud Restringida (REML) para estimar los parámetros del modelo.
Distribución Normal: Tanto los residuos como los efectos aleatorios deben seguir una distribución normal para que el modelo sea válido.
Un modelo de intersecciones aleatorias es un tipo de modelo lineal mixto que incluye efectos aleatorios para capturar la variabilidad entre diferentes grupos o sujetos. En estos modelos, se asume que cada grupo o sujeto tiene su propio intercepto, que varía aleatoriamente alrededor de un valor promedio común.
Componentes del Modelo de Intersecciones Aleatorias
Efectos Fijos: Representan la relación promedio entre las variables predictoras y la variable respuesta en toda la población.
Efectos Aleatorios: Capturan la variabilidad específica de cada grupo o sujeto. En el caso de intersecciones aleatorias, solo el intercepto varía entre los grupos, mientras que las pendientes se mantienen constantes.
2.4.1 Ejemplo: GPA del estudiante
Para nuestro primer ejemplo, evaluaremos los factores que predicen el promedio de calificaciones (GPA) de la universidad. Cada uno de los 200 estudiantes es evaluado en seis ocasiones (cada semestre durante los primeros tres años), por lo que tenemos observaciones agrupadas dentro de los estudiantes. Tenemos otras variables, como la situación laboral, el sexo y el GPA de la escuela secundaria. Algunas estarán en forma tanto numérica como etiquetada.
2.4.2 El modelo de regresión estándar
Ahora veamos el modelo subyacente. Podemos mostrarlo de un par de maneras diferentes. Primero, comenzamos con una regresión estándar para orientarnos.
\[\mathrm{gpa} = b_{\mathrm{intercept}} + b_{\mathrm{occ}}\cdot \mathrm{occasion} + \epsilon\]
Disponemos de coeficientes (\(b_*\)) para la intersección y el efecto del tiempo. Se supone un error (\(\epsilon\)) distribuído normalmente con media 0 y desviación estándar \(\sigma\).
\[\epsilon \sim \mathcal{N}(0, \sigma)\]
Una forma alternativa de escribir el modelo que pone énfasis en el proceso de generación de datos subyacente para \(\mathrm{gpa}\) se puede mostrar de la siguiente manera.
\[\mathrm{gpa} \sim \mathcal{N}(\mu, \sigma)\] \[\mu = b_{\mathrm{intercept}} + b_{\mathrm{occ}}\cdot \mathrm{occasion}\]
Más técnicamente, las variables GPA y \(\mu\) tienen un subíndice implícito para denotar cada observación, pero también puedes pensarlo como un modelo para un solo individuo en un solo punto temporal.
2.5 El modelo mixto
2.5.1 Representación inicial
Ahora mostramos una forma de representar un modelo mixto que incluye un efecto único para cada estudiante. Considere el siguiente modelo para un solo estudiante. Esto demuestra que el efecto específico del estudiante, es decir, la desviación en el promedio de calificaciones solo para ese estudiante siendo quien es, puede verse como una fuente adicional de varianza.
\[\mathrm{gpa} = b_{\mathrm{intercept}} + b_{\mathrm{occ}}\cdot \mathrm{occasion} + (\mathrm{effect}_{\mathrm{student}} + \epsilon)\]
We would (usually) assume the following for the student effects.
\[\mathrm{effect}_{\mathrm{student}} \sim \mathcal{N}(0, \tau)\] La diferencia principal entre este modelo mixto y una regresión estándar es el efecto del estudiante. En promedio, este efecto es cero, pero varía de un estudiante a otro con una desviación estándar (\(\tau\)).
Entonces, los efectos del estudiante son aleatorios y, específicamente, están distribuidos normalmente con una media de cero y una desviación estándar estimada (\(\tau\)). En otras palabras, conceptualmente, la única diferencia entre este modelo mixto y una regresión estándar es el efecto del estudiante, que en promedio es nulo, pero típicamente varía de un estudiante a otro en una cantidad que en promedio es (\(\tau\)).
Si lo reorganizamos, podemos enfocarnos en los coeficientes del modelo, en lugar de verlo como una fuente adicional de error.
\[\mathrm{gpa} = (b_{\mathrm{intercept}} + \mathrm{effect}_{\mathrm{student}}) + b_{\mathrm{occ}}\cdot \mathrm{occasion} + \epsilon\] O más sucintamente:
\[\mathrm{gpa} = b_{\mathrm{int\_student}} + b_{\mathrm{occ}}\cdot \mathrm{occasion} + \epsilon\]
De esta manera, tendremos interceptos específicos para cada estudiante, ya que cada persona tendrá su propio efecto único añadido al intercepto general, resultando en un intercepto diferente para cada persona.
\[b_{\mathrm{int\_student}} \sim \mathcal{N}(b_{\mathrm{intercept}}, \tau)\]
Ahora vemos los interceptos como distribuidos normalmente con una media del intercepto general y una desviación estándar. Por lo tanto, esto a menudo se llama un modelo de interceptos aleatorios.
2.5.2 Como un modelo multinivel
Otra forma de mostrar el modelo mixto es comúnmente vista en la literatura de modelos multinivel. Se muestra de manera más explícita como un modelo de regresión en dos partes, una a nivel de observación y otra a nivel de estudiante.
\[\mathrm{gpa} = b_{\mathrm{int\_student}} + b_{\mathrm{occ}}\cdot \mathrm{occasion} + \epsilon\]
\[b_{\mathrm{int\_student}} = b_{\mathrm{intercept}} + \mathrm{effect}_{\mathrm{student}}\]
Sin embargo, después de “insertar” la parte del segundo nivel en la primera, es idéntico al modelo anterior.
Observa cómo no tenemos un efecto específico del estudiante para la ocasión. En este contexto, se dice que la ocasión es un efecto fijo únicamente, y no hay un componente aleatorio. Sin embargo, esto definitivamente no tiene por qué ser así, como veremos más adelante.
2.6 Aplicación
2.6.1 Visualización inicial
Graficamos el promedio de calificaciones en función de la ocasión (es decir, el semestre) para tener una idea de la variabilidad en los puntos de partida y las tendencias.
Todas las trayectorias de los estudiantes se muestran como trayectorias atenuadas. La tendencia general, estimada mediante la regresión, se muestra en rojo. Dos cosas se destacan. Una es que los estudiantes tienen mucha variabilidad al comenzar. En segundo lugar, si bien la tendencia general en el GPA es ascendente con el tiempo como esperaríamos, los estudiantes individuales pueden variar en esa trayectoria.
2.6.2 Regresión estándar
Primero, veremos la regresión y sólo el indicador de tiempo como covariable, que trataremos como numérico para simplificar.
##
## Call:
## lm(formula = gpa ~ occasion, data = gpa)
##
## Coefficients:
## (Intercept) occasion
## 2.5992 0.1063load('data/gpa.RData')
gpa_lm = lm(gpa ~ occasion, data = gpa)
## summary(gpa_lm)| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 2.599 | 0.018 | 145.7 | 0 |
| occasion | 0.106 | 0.006 | 18.04 | 0 |
| Observations | Residual Std. Error | \(R^2\) | Adjusted \(R^2\) |
|---|---|---|---|
| 1200 | 0.3487 | 0.2136 | 0.2129 |
Lo anterior nos indica que, al comenzar, es decir, cuando la ocasión es cero, el promedio de GPA, denotado por la intersección, es 2.6. Además, a medida que avanzamos de un semestre a otro, podemos esperar que el GPA aumente en aproximadamente 0.11 puntos. Esto estaría bien, excepto que ignoramos la agrupación. Un efecto secundario de hacerlo es que nuestros errores estándar están sesgados y, por lo tanto, las afirmaciones sobre la significación estadística basadas en ellos estarían equivocadas. Sin embargo, lo más importante es que simplemente no podemos explorar el efecto del estudiante, que sería interesante por sí mismo.
Un enfoque alternativo que podríamos adoptar sería ejecutar regresiones separadas para cada estudiante. Sin embargo, esto tiene muchas desventajas: no es fácil resumirlo cuando hay muchos grupos, normalmente habría muy pocos datos dentro de cada conglomerado para hacerlo (como en este caso) y los modelos están sobrecontextualizados, lo que significa que ignoran lo que los estudiantes tienen en común. Compararemos este enfoque con el modelo mixto más adelante.
2.6.3 3 Ejecución de un modelo mixto
A continuación, ejecutamos un modelo mixto que permitirá un efecto
específico del estudiante. Este tipo de modelo se puede llevar a cabo
fácilmente en R, específicamente con el paquete lme4. A continuación, el código se verá igual que el
que usaste para la regresión con lm, pero con
un componente adicional que especifica el efecto del grupo, es decir, el
estudiante. El (1|estudiante) significa que estamos
permitiendo que la intersección, representada por 1, varíe
según el estudiante. Con el modelo mixto, obtenemos los mismos
resultados que la regresión, pero como veremos, tendremos más de qué
hablar.
library(lme4)
gpa_mixed = lmer(gpa ~ occasion + (1 | student), data = gpa)
## summary(gpa_mixed)| term | value | se | t | p_value | lower_2.5 | upper_97.5 |
|---|---|---|---|---|---|---|
| Intercept | 2.599 | 0.022 | 119.800 | 0 | 2.557 | 2.642 |
| occasion | 0.106 | 0.004 | 26.096 | 0 | 0.098 | 0.114 |
| group | effect | variance | sd |
|---|---|---|---|
| student | Intercept | 0.064 | 0.252 |
| Residual | NA | 0.058 | 0.241 |
En primer lugar, vemos que los coeficientes (o en este contexto, los efectos fijos) para el intercepto y el tiempo son los mismos que vimos con la regresión estándar [^lmlmercoef], como sería su interpretación. Los errores estándar, por otra parte, son diferentes aquí, aunque al final nuestra conclusión sería la misma en lo que respecta a la significación estadística. Nótese específicamente que el error estándar para el intercepto ha aumentado. Conceptualmente, se puede pensar que permitir interceptos aleatorios por persona nos permite obtener información sobre el individuo, al tiempo que reconocemos la incertidumbre con respecto al promedio general que estábamos subestimando antes [^sewithin].
Si bien tenemos coeficientes y errores estándar, es posible que haya notado que lme4 no proporciona valores p. Hay varias razones para esto, a saber, que con los modelos mixtos estamos tratando esencialmente con diferentes tamaños de muestra, los \(N_c\) dentro de los conglomerados, que pueden variar de un conglomerado a otro (¡e incluso ser una sola observación!), y N observaciones totales, lo que nos coloca en una especie de situación difusa con respecto a las distribuciones de referencia, los grados de libertad del denominador y cómo aproximarse a una “mejor” solución. Otros programas proporcionan valores p automáticamente como si no hubiera ningún problema y sin decirle qué enfoque usan para calcularlos (hay varios). Además, esas aproximaciones pueden ser muy deficientes en algunos escenarios o hacer suposiciones que pueden no ser apropiadas para la situación[^fuzzyp].
Sin embargo, es más sencillo obtener intervalos de confianza y podemos hacerlo con lme4 de la siguiente manera[^confint].
confint(gpa_mixed)## Computing profile confidence intervals ...| group | effect | variance | sd | sd_2.5 | sd_97.5 | var_prop |
|---|---|---|---|---|---|---|
| student | Intercept | 0.064 | 0.252 | 0.225 | 0.282 | 0.523 |
| Residual | NA | 0.058 | 0.241 | 0.231 | 0.252 | 0.477 |
2.6.3.1 Componentes de la varianza
Una cosa nueva en comparación con la salida de una regresión estándar es la desviación estándar/varianza estimada del efecto del estudiante (\(\tau\)/\(\tau^2\) en nuestra representación de la fórmula anterior). Esto nos indica cuánto, en promedio, varía el GPA al pasar de un estudiante a otro. En otras palabras, incluso después de hacer una predicción basada en el punto temporal, cada estudiante tiene su propia desviación única, y ese valor (en términos de la desviación estándar) es la desviación promedio estimada entre los estudiantes. Es importante notar que las calificaciones varían debido al estudiante más del doble de lo que varían debido a un cambio de semestre. Este es un aspecto interpretativo importante que no está disponible en un modelo de regresión estándar.
Another way to interpret the variance output is to note percentage of the student variance out of the total, or 0.064 / 0.122 = 52%. In this setting, this value is also called the intraclass correlation, because it is also an estimate of the within cluster correlation, as we’ll see later.
2.6.3.2 Estimates of the random effects
After running the model, we can actually get estimates of the student effects[^blup]. I show two ways for the first five students, both as random effect and as random intercept (i.e. intercept + random effect).
ranef(gpa_mixed)$student %>% head(5)
# showing mixedup::extract_random_effects(gpa_mixed)| group_var | effect | group | value | se | lower_2.5 | upper_97.5 |
|---|---|---|---|---|---|---|
| student | Intercept | 1 | -0.071 | 0.092 | -0.251 | 0.109 |
| student | Intercept | 2 | -0.216 | 0.092 | -0.395 | -0.036 |
| student | Intercept | 3 | 0.088 | 0.092 | -0.091 | 0.268 |
| student | Intercept | 4 | -0.187 | 0.092 | -0.366 | -0.007 |
| student | Intercept | 5 | 0.030 | 0.092 | -0.149 | 0.210 |
coef(gpa_mixed)$student %>% head(5)| group_var | effect | group | value | se | lower_2.5 | upper_97.5 |
|---|---|---|---|---|---|---|
| student | Intercept | 1 | 2.528 | 0.095 | 2.343 | 2.713 |
| student | Intercept | 2 | 2.383 | 0.095 | 2.198 | 2.568 |
| student | Intercept | 3 | 2.687 | 0.095 | 2.502 | 2.872 |
| student | Intercept | 4 | 2.412 | 0.095 | 2.227 | 2.597 |
| student | Intercept | 5 | 2.629 | 0.095 | 2.444 | 2.814 |
Note that we did not allow occasion to vary, so it is a constant, i.e. fixed, effect for all students.
Often, we are keenly interested in these effects, and want some sense of uncertainty regarding them. With lme4 this typically would be done via bootstrapping, specifically with the bootMer function within lme4. However, for some users this may be a bit of a more complex undertaking. The merTools package provides an easier way to get this with the predictInterval function[^predinterval]. Or you can go straight to the plot of them.
library(merTools)
predictInterval(gpa_mixed) # for various model predictions, possibly with new data
REsim(gpa_mixed) # mean, median and sd of the random effect estimates
plotREsim(REsim(gpa_mixed)) # plot the interval estimatesEl siguiente gráfico muestra los efectos aleatorios estimados para cada estudiante y su estimación de intervalo (una versión modificada del gráfico producido por esa última línea de código [^mertoolsplotlabels]). Recuerde que los efectos aleatorios se distribuyen normalmente con una media de cero, que se muestra mediante la línea horizontal. Los intervalos que no incluyen cero están en negrita. En este caso, dichos estudiantes tienen un resultado inicial relativamente más alto o más bajo en comparación con un estudiante típico.
##
## Adjuntando el paquete: 'visibly'## The following object is masked from 'package:maditr':
##
## :=
2.6.3.3 Prediction
Let’s now examine standard predictions vs. cluster-specific predictions. As with most R models, we can use the predict function on the model object.
predict(gpa_mixed, re.form=NA) %>% head()## 1 2 3 4 5 6
## 2.599214 2.705529 2.811843 2.918157 3.024471 3.130786In the above code we specified not to use the random effects
re.form=NA, and as such, our predictions for the
observations are pretty much what we’d get from the standard linear
model.
predict_no_re = predict(gpa_mixed, re.form=NA)
predict_lm = predict(gpa_lm)But each person has their unique intercept, so let’s see how the predictions differ when we incorporate that information.
predict_with_re = predict(gpa_mixed)Depending on the estimated student effect, students will start above or below the estimated intercept for all students. The following visualizes the unconditional prediction vs. the conditional prediction that incorporates the random intercept for the first two students.
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levelsWe can see that the predictions from the mixed model are shifted because of having a different intercept. For these students, the shift reflects their relatively poorer start.
2.7 Cluster Level Covariates
Note our depiction of a mixed model as a multilevel model.
\[\mathrm{gpa} = b_{\mathrm{int\_student}} + b_{\mathrm{occ}}\cdot \mathrm{occasion} + \epsilon\]
\[b_{\mathrm{int\_student}} = b_{\mathrm{intercept}} + \mathrm{effect}_{\mathrm{student}}\] If we add student a student level covariate, e.g sex, to the model, we then have the following.
\[b_{\mathrm{int\_student}} = b_{\mathrm{intercept}} + b_{\mathrm{sex}}\cdot \mathrm{sex} + \mathrm{effect}_{\mathrm{student}}\]
Which, after plugging in, we still have the same model as before, just with an additional predictor.
\[\mathrm{gpa} = b_{\mathrm{intercept}} + b_{\mathrm{occ}}\cdot \mathrm{occasion}+ b_{\mathrm{sex}}\cdot \mathrm{sex} + (\mathrm{effect}_{\mathrm{student}} + \epsilon)\]
Por lo tanto, al final, agregar covariables a nivel de clúster no tiene ningún efecto inusual en cómo pensamos acerca del modelo[^mlevel]. Simplemente las agregamos a nuestro conjunto de variables predictoras. Tenga en cuenta también que podemos crear covariables a nivel de clúster como medias de grupo o algún otro resumen de las variables a nivel de observación. Esto es especialmente común cuando los clústeres representan unidades geográficas y las observaciones son personas. Por ejemplo, podríamos tener el ingreso como una covariable a nivel de persona y usar la mediana para representar la riqueza general de la región geográfica.
2.8 Resumen de los conceptos básicos de los modelos mixtos
Los modelos mixtos nos permiten tener en cuenta la estructura observada en los datos. Si solo se utilizaran para eso, tendríamos una inferencia más precisa en relación con lo que obtendríamos si ignoráramos esa estructura. Sin embargo, ¡obtenemos mucho más! Entendemos mejor las fuentes de variabilidad en la variable objetivo. También obtenemos estimaciones específicas de los parámetros del modelo para cada grupo, lo que nos permite comprender exactamente cómo difieren los grupos entre sí. Además, esto a su vez permite una predicción específica para cada grupo y, por lo tanto, una predicción mucho más precisa, suponiendo que haya una varianza apreciable debido a la agrupación. En resumen, los modelos mixtos tienen mucho que ganar, incluso en los entornos más simples.
2.9 Exercises for Starting Out
2.9.1 Sleep
For this exercise, we’ll use the sleep study data from the lme4 package. The following describes it.
The average reaction time per day for subjects in a sleep deprivation study. On day 0 the subjects had their normal amount of sleep. Starting that night they were restricted to 3 hours of sleep per night. The observations represent the average reaction time (in milliseconds) on a series of tests given each day to each subject.
After loading the package, the data can be loaded as follows. I show the first few observations.
datatable(sleepstudy, options = list(
pageLength = 10,
autoWidth = TRUE,
initComplete = JS(
"function(settings, json) {",
"$('table.dataTable').css({'font-size': '10px', 'width': '80%'});",
"}"
)
))Ejecuta una regresión con Reaction como la variable objetivo y Days como el predictor.
Ejecuta un modelo mixto con un intercepto aleatorio para Subject.
Interpreta los componentes de varianza y los efectos fijos.
2.9.2 Añadiendo la covariable a nivel de clúster
Vuelve a ejecutar el modelo mixto con los [GPA data][Mixed model]
añadiendo la covariable a nivel de clúster sex, o GPA de la
escuela secundaria (highgpa), o ambos. Interpreta todos los aspectos de
los resultados.
What happened to the student variance after adding cluster level covariates to the model?
2.9.3 Simulating a mixed model
The following represents a simple way to simulate a random intercepts model. Note each object what each object is, and make sure the code make sense to you. Then run it.
set.seed(1234) # this will allow you to exactly duplicate your result
Ngroups = 100
NperGroup = 3
N = Ngroups * NperGroup
groups = factor(rep(1:Ngroups, each = NperGroup))
u = rnorm(Ngroups, sd = .5)
e = rnorm(N, sd = .25)
x = rnorm(N)
y = 2 + .5 * x + u[groups] + e
d = data.frame(x, y, groups)Which of the above represent the fixed and random effects? Now run the following.
model = lmer(y ~ x + (1|groups), data=d)
summary(model)
confint(model)
library(ggplot2)
ggplot(aes(x, y), data=d) +
geom_point()Do the results seem in keeping with what you expect?
In what follows we’ll change various aspects of the data, then rerun the model after each change, then summarize and get confidence intervals as before. For each note specifically at least one thing that changed in the results.
- First calculate or simply eyeball the intraclass correlation coefficient \(\frac{\textrm{random effect variance}}{\textrm{residual + random effect variance}}\). In addition, create a density plot of the random effects as follows.