-
Notifications
You must be signed in to change notification settings - Fork 6
/
Copy path303-apli-lme4.Rmd
214 lines (161 loc) · 10.8 KB
/
303-apli-lme4.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
# Aplicación con `lme4` {#apli-lme4}
En este capítulo se mostrará como usar el paquete `lme4` para la aplicación de modelos mixtos con la base de datos `sleepstudy` del mismo paquete.
```{r echo=FALSE, out.width="40%", fig.align='center'}
knitr::include_graphics("images/sleep_study.png")
```
A continuación la base de datos a utilizar.
```{r, message = FALSE}
library(lme4)
head(sleepstudy)
```
Esta base de datos sobre el tiempo de reacción promedio por día para un conjunto de individuos, en un estudio de privación del sueño, contiene la información sobre el tiempo de reacción promedio (`Reaction`), el número de días de privación del sueño (`Days`), donde el día 0 corresponde al día en el que los indiviuos tenían su cantidad normal de sueño, y el número del individuo (en total 18) sobre el que se realizó la observación (`Subject`). A partir del día 0, hubo una restricción en cada individuo a 3 horas de sueño por noche.
```{r plot_sleepstudy_lme4, fig.align="center"}
library(ggplot2)
ggplot(data = sleepstudy, aes(x = Days, y = Reaction, color = Subject)) +
geom_point() +
theme_bw() +
facet_wrap(~ Subject) + labs(y = "Reaction time") +
theme(legend.position = "none")
```
De la figura anterior vemos que el tiempo de reacción promedio, tanto en el día 0 como en los siguientes días de prueba (del día 1 al día 9), son distintos en cada uno de los individuos. Esta situación conlleva a probar la hipótesis de que el tiempo de reacción promedio en una serie de pruebas varía según los individuos. Esto es, ajustar un modelo donde el intercepto y la pendiente se consideran como efectos aleatorios.
Un modelo lineal mixto que describe la anterior situación se puede escribir como:
\begin{align*}
Reaction_{ij} | b_0, b_1 &\sim N(\mu_{ij}, \sigma^2_{Reaction}) \\
\mu_{ij} &= \beta_0 + \beta_1 Days_{ij} + b_{0i} + b_{1i} Days_{ij} \\
\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} & \sigma_{b01} \\
\sigma_{b01} & \sigma^2_{b1}
\end{matrix} \right ]
\right )
\end{align*}
Aquí, los individuos ($i$) varían en el tiempo de reacción promedio tanto en su intercepto ($b_{0i}$) como en su pendiente ($b_{1i}$), que en conjunto componen la varianza total en dicho tiempo atribuible a la variación entre individuos. Esta contribución individual se cuantifica usando un modelo de intercepto y pendiente aleatoria con distribución normal ($N$). La variación entre individuos en intercepto y pendiente es $\sigma^2_{b0}$ y $\sigma^2_{b1}$, respectivamente. La covarianza entre el intercepto y la pendiente esta dada por $\sigma_{b01}$.
El vector de parámetros para este modelo sería $\boldsymbol{\Theta}=(\beta_0, \beta_1, \sigma_{reaction}, \sigma_{b0}, \sigma_{b1}, \sigma_{b0b1})^\top$.
Para ajustar el modelo de intercepto y pendiente aleatoria planteado usando el paquete `lme4` podemos usar el siguiente código:
```{r}
fit <- lmer(Reaction ~ Days + (Days | Subject), REML = TRUE, data = sleepstudy)
```
Para obtener la tabla de resumen usamos:
```{r}
summary(fit)
```
De la salida anterior se obtienen los siguientes parámetros estimados
$$
\hat{\boldsymbol{\Theta}} = (\hat{\beta_0}=251.40, \hat{\beta_1}=10.47, \hat{\sigma}_{reaction}=25.59, \hat{\sigma}_{b0}=24.74, \hat{\sigma}_{b1}=5.92, \hat{\sigma}_{b0b1}=10.25)^\top
$$
El último parámetro estimado se obtiene utilizando la ecuación de [correlación](https://en.wikipedia.org/wiki/Correlation_and_dependence) que relaciona la correlación, la covarianza y las desviaciones de los efectos aleatorios: $\rho_{b0b1} = \sigma_{b0b1}/(\sigma_{b0} \times \sigma_{b1})$
Usando la información anterior se puede escribir el modelo ajustado de la siguiente manera:
\begin{align*}
Reaction_{ij} &\sim N(\hat{\mu}_{ij}, 25.59^2) \\
\hat{\mu}_{ij} &= 251.40 + 10.47 Days_{ij} + \tilde{b}_{0i} + \tilde{b}_{1i} Days_{ij} \\
\left (
\begin{matrix}
b_{0} \\ b_{1}
\end{matrix}
\right ) &\sim
N\left ( \left [ \begin{matrix}
0 \\ 0
\end{matrix} \right ],
\left [ \begin{matrix}
24.74^2 & 10.25 \\
10.25 & 5.92^2
\end{matrix} \right ]
\right )
\end{align*}
Los elementos $b_{0i}$ y $b_{1i}$ se deben substituir por sus respectivas predicciones $\tilde{b}_{0i}$ y $\tilde{b}_{1i}$ y se pueden obtener del modelo ajustado de esta forma:
```{r Prediciones de los efectos aleatorios}
ranef(fit)
```
Y los valores de los efectos fijos estimados se pueden obtener así:
```{r Valores de los efectos fijos estimados}
fixef(fit)
```
Con base en la información anterior de efectos aleatorios y fijos, es posible escribir la ecuación del modelo para cada individuo. Para esto, se debe considerar los efectos fijos estimados ($\hat{\beta}_0 \approx 251.40$ y $\hat{\beta}_1\approx 10.47$) y los efectos aleatorios de cada uno de los individuos (por ejemplo para el individuo `308`, $\tilde{b}_{0, i=308} \approx 2.26$ y $\tilde{b}_{1, i=308} \approx 9.20$). Así, el valor medio del individuo `308` se calcula como:
\begin{align*}
\hat{\mu}_{i=308, j} &= \hat{\beta}_0 + \hat{\beta}_0 \, Days_{i=308, j} + \tilde{b}_{0, i=308} + \tilde{b}_{1, i=308} \, Days_{i=308, j} \\
\hat{\mu}_{i=308, j} &= 251.40 + 10.47 \, Days_{i=308, j} 2.26 + 9.20 \, Days_{i=308, j} \\
\hat{\mu}_{i=308, j} &= 253.66 + 19.67 \, Days_{i=308, j}
\end{align*}
Lo anterior se puede resumir en el siguiente modelo.
\begin{align*}
Reaction_{i=308, j} &\sim N(\hat{\mu}_{i=308, j}, \hat{\sigma}^2_{Reaction}) \\ \hat{\mu}_{i=308, j} &= 253.66 + 19.67 \, Days_{i=308, j} \\
\hat{\sigma}_{Reaction} &= 25.59
\end{align*}
Los efectos fijos y aleatorios de la expresión anterior para cada uno de los individuos se pueden obtener con `R` de la siguiente forma:
```{r Efectos fijos y aleatorios para cada uno de los individuos}
coef(fit)
```
A continuación usted podrá observar el diagrama de dispersión mostrado al inicio de este capitulo con la recta de regresión para cada individuo. El código de `R` para obtener esto se presenta a continuación:
```{r Gráfico con recta de regresión, fig.align = "center"}
fit <- lmer(Reaction ~ Days + (Days | Subject), REML = TRUE, data = sleepstudy)
sleepstudy$pred_inter_pend_aleatorio <- predict(fit)
ggplot(data = sleepstudy, aes(x = Days, y = pred_inter_pend_aleatorio, color = Subject)) +
geom_line() +
geom_point(aes(x = Days, y = Reaction, color = Subject)) +
geom_abline(intercept = 251.40, slope = 10.47, color = "black", linetype = "dashed", linewidth = 0.5) +
theme_bw() +
facet_wrap(~ Subject) +
theme(legend.position = "none") + labs(y = "Reaction time")
```
La figura anterior corresponde a un modelo de intercepto y pendiente aleatoria, en el que se permite que tanto los interceptos como las pendientes varíen según los individuos. Las líneas continuas corresponde a la recta de regresión ajustada a los datos. Los puntos representan las observaciones (tiempo de reacción promedio por día) medidas en cada uno de los individuos. La línea negra discontinua representa el valor medio global de la distribución de los efectos aleatorios.
A continuación usted encontrará una figura con cuatro páneles donde se muestran los resultados de cuatro modelos distintos, entre ellos el modelo mixto con intercepto y pendiente aleatoria del ejemplo anterior. Con base en estas figuras, responda las preguntas siguientes.
```{r Graficos de los distintos modelo, echo = FALSE, eval = TRUE, message = FALSE, fig.align="center"}
library('see')
fit_1 <- lm(Reaction ~ Days, data = sleepstudy)
sleepstudy$modelo_simple <- predict(fit_1)
Graf_1 <- ggplot(data = sleepstudy, aes(x = Days, y = modelo_simple, color = Subject)) +
geom_point(aes(x = Days, y = Reaction, color = Subject)) +
geom_abline(intercept = 251.40, slope = 10.47, color = "black", linetype = "solid", size = 0.5) +
theme_bw() +
theme(legend.position = "none") +
labs(x = "Días", y = "Tiempo de reacción", title = 'Modelo de regresión clásico') +
theme(axis.text = element_text(size = 8, color = "black", face = "bold"),
axis.title = element_text(size = 8, face = "bold"),
plot.title = element_text(size = 8, color = "black", face = "bold"))
fit_2 <- lmer(Reaction ~ Days + (1 | Subject), data = sleepstudy)
sleepstudy$pred_inter_aleatorio <- predict(fit_2)
Graf_2 <- ggplot(data = sleepstudy, aes(x = Days, y = pred_inter_aleatorio, color = Subject)) +
geom_line() +
geom_point(aes(x = Days, y = Reaction, color = Subject)) +
theme_bw() +
theme(legend.position = "none") +
labs(x = "Días", y = "Tiempo de reacción", title = 'Modelo mixto con intercepto aleatorio') +
theme(axis.text = element_text(size = 8, color = "black", face = "bold"),
axis.title = element_text(size = 8, face = "bold"),
plot.title = element_text(size = 8, color = "black", face = "bold"))
fit_3 <- lmer(Reaction ~ Days + (0 + Days | Subject), data = sleepstudy)
sleepstudy$pred_pend_aleatorio <- predict(fit_3)
Graf_3 <- ggplot(data = sleepstudy, aes(x = Days, y = pred_pend_aleatorio, color = Subject)) +
geom_line() +
geom_point(aes(x = Days, y = Reaction, color = Subject)) +
theme_bw() +
theme(legend.position = "none") +
labs(x = "Días", y = "Tiempo de reacción", title = 'Modelo mixto con pendiente aleatoria') +
theme(axis.text = element_text(size = 8, color = "black", face = "bold"),
axis.title = element_text(size = 8, face = "bold"),
plot.title = element_text(size = 8, color = "black", face = "bold"))
fit_4 <- lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy)
sleepstudy$pred_inter_pend_aleatorio <- predict(fit_4)
Graf_4 <- ggplot(data = sleepstudy, aes(x = Days, y = pred_inter_pend_aleatorio, color = Subject)) +
geom_line() +
geom_point(aes(x = Days, y = Reaction, color = Subject)) +
theme_bw() +
theme(legend.position = "none") +
labs(x = "Días", y = "Tiempo de reacción", title = 'Modelo mixto con intercepto y pendiente\n aleatoria') +
theme(axis.text = element_text(size = 8, color = "black", face = "bold"),
axis.title = element_text(size = 8, face = "bold"),
plot.title = element_text(size = 8, color = "black", face = "bold"))
plots(Graf_1, Graf_2, Graf_3, Graf_4, n_columns = 2,
tags = paste0("Panel ", 1:4))
```
## Ejercicios {-}
1. Ajuste el modelo con intercepto aleatorio mostrado en el panel 2. ¿Qué opina de este modelo?
2. Ajuste el modelo con pendiente aleatoria presentado en el panel 3. ¿Qué opina de este modelo?
3. Ajustar solo un intercepto aleatorio permite que los individuos varíen asumiendo que los mismos tienen una pendiente común (panel 2). Al ajustar solo una pendiente aleatoria (panel 3) permite que la pendiente de un predictor varíe en función de los individuos (la variable de agrupación). Con base esto y teniendo en cuenta el modelo de intercepto y pendiente aleatoria (panel 4), evalúe cual de estos estos modelos permite un mejor ajuste de los datos presentados en la base de datos `sleepstudy` del paquete `lme4`.