-
Notifications
You must be signed in to change notification settings - Fork 6
/
305-apli-nlme.Rmd
149 lines (111 loc) · 6.02 KB
/
305-apli-nlme.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
# Aplicación con `nlme` {#apli-nlme}
En este capítulo se mostrará como usar el paquete `nlme` para la aplicación de modelos mixtos con la base de datos `Oxboys` del mismo paquete.
```{r echo=FALSE, out.width="80%", fig.align='center'}
knitr::include_graphics("images/crecimiento_chicos.png")
```
A continuación la base de datos a utilizar.
```{r}
library(nlme)
head(Oxboys)
```
Esta base de datos sobre crecimiento contiene la información sobre altura (`heigth`), edad estandarizada (`age`) de un grupo de 26 jóvenes. Como la base de datos `Oxboys` es de la clase groupedData, es posible aplicar un `plot` directamente y el resultado se muestra continuación.
```{r plot_oxboys_nlme, fig.width=12}
plot(Oxboys)
```
<div style="-moz-box-shadow: 1px 1px 3px 2px #00ffff;
-webkit-box-shadow: 1px 1px 3px 2px #00ffff;
box-shadow: 1px 1px 3px 2px #00ffff;">
```{block2, type='rmdtip'}
Es posible convertir un data.frame para que tenga la clase groupedData, consulte la ayuda de la función `groupedData` del paquete `nlme` para más detalles.
```
</div>
De la figura anterior vemos que las curvas de crecimiento inician a diferente altura (intercepto) y que la pendiente del crecimiento no son todas iguales, por ejemplo, el individuo 21 creció más rápido que el individuo 3. Esto nos hace pensar que un modelo con intercepto y pendiente aleatoria podrían ser adecuados para modelar el crecimiento.
En las siguientes ecuaciones se resume el modelo matemático que interesa en esta situación.
\begin{align*}
Height_{ij} | b_0, b_1 &\sim N(\mu_{ij}, \sigma^2_{Heigth}) \\
\mu_{ij} &= \beta_0 + \beta_1 Age_{ij} + b_{0i} + b_{1i} Age_{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*}
El vector de parámetros para este modelo sería $\boldsymbol{\Theta}=(\beta_0, \beta_1, \sigma_{height}, \sigma_{b0}, \sigma_{b1}, \sigma_{b0b1})^\top$.
Para ajustar este modelo a los datos con el paquete `nlme` podemos usar el siguiente código.
```{r}
fit <- lme(height ~ age, random= ~ 1 + age | Subject, data=Oxboys, method="REML")
```
Para obtener la tabla de resumen usamos:
```{r}
summary(fit)
```
De la salida anterior se obtiene que $\hat{\boldsymbol{\Theta}}=(\hat{\beta}_0=149.37, \hat{\beta}_1=6.53, \hat{\sigma}_{height}=0.66, \hat{\sigma}_{b0}=8.08, \hat{\sigma}_{b1}=1.68, \hat{\sigma}_{b0b1}=8.71)^\top$. La estimación $\hat{\sigma}_{b0b1}$ no aparece directamente en el summary pero se obtiene utilizando la ecuación $Cor=Cov/(\sigma_1 \sigma_2)$ que relaciona [correlación](https://en.wikipedia.org/wiki/Correlation_and_dependence), covarianza y desviaciones de los efectos aleatorios.
Usando la información anterior se puede escribir el modelo ajustado de la siguiente manera.
\begin{align*}
Height_{ij} &\sim N(\hat{\mu}_{ij}, 0.66^2) \\
\hat{\mu}_{ij} &= 149.37 + 6.53 Age_{ij} + b_{0i} + b_{1i} Age_{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}
8.08^2 & 8.71 \\
8.71 & 1.68^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 así:
```{r}
ranef(fit)
```
Los valores de los efectos fijos estimados se pueden obtener así:
```{r}
fixef(fit)
```
Usando la información de los efectos fijo y aleatorios obtenidos antes, es posible escribir la ecuación del modelo para cada individuo. Los efectos fijos estimados fueron $\hat{\beta}_0 \approx 149.37$ y $\hat{\beta}_1\approx 6.53$. Para el sujeto # 10 se obtuvo $\tilde{b}_{0, i=10} \approx -19.10$ y $\tilde{b}_{1, i=10} \approx -2.79$, así la media del individuo # 10 se calcula así:
\begin{align*}
\hat{\mu}_{i=10, j} &= \hat{\beta}_0 + \hat{\beta}_0 \, Age_{i=10, j} + \tilde{b}_{0, i=10} + \tilde{b}_{1, i=10} \, Age_{i=10, j} \\
\hat{\mu}_{i=10, j} &= 149.37 + 6.53 \, Age_{i=10, j} - 19.10 - 2.79 \, Age_{i=10, j} \\
\hat{\mu}_{i=10, j} &= 130.27 + 3.74 \, Age_{i=10, j}
\end{align*}
Lo anterior se puede resumir en el siguiente modelo.
\begin{align*}
Height_{i=10, j} &\sim N(\hat{\mu}_{i=10, j}, \hat{\sigma}^2_{Height}) \\
\hat{\mu}_{i=10, j} &= 130.27 + 3.74 \, Age_{i=10, j} \\
\hat{\sigma}_{Height} &= 0.66
\end{align*}
La expresión anterior para cada individuo con los efectos finales (fijos y aleatorios) se puede obtener con R así:
```{r}
coef(fit)
```
En la presente aplicación es posible incluir la recta de regresión para cada individuo al diagrama de dispersión original. El código de R para obtener esto es el siguiente.
```{r evol_fit_nlme}
library(lattice)
xyplot(height ~ age | Subject, data=Oxboys, fit=fit,
strip=strip.custom(bg="white"),
pch=16, cex=0.7, col='indianred1',
panel = function(x, y, ..., fit, subscripts) {
panel.xyplot(x, y, ...)
ypred <- fitted(fit)[subscripts]
panel.lines(x, ypred, col="deepskyblue3", lwd=1)
},
ylab="Height (cm)", xlab="Centered age")
```
En la figura anterior se tienen las observaciones (crecimiento) representado por los puntos rojos, adicionalmente, aparece una recta de color azul que representa la recta de regresión para cada individuo. De la figura se observa que la linea logra explicar la evolución del crecimiento para cada individuo.
## Ejercicios {-}
1. Repita el ejercicio anterior considerando un modelo sólo con intercepto aleatorio. Dibuje las rectas de regresión para cada individuo. ¿Qué opina de este modelo?
2. Repita el ejercicio anterior considerando un modelo sólo con pendiente aleatoria. Dibuje las rectas de regresión para cada individuo. ¿Qué opina de este modelo?
3. Estime la estatura para el individuo # 3 cuando su edad centrada sea de 1.1.
4. Replique los ejemplo [de este documento](https://socialsciences.mcmaster.ca/jfox/Books/Companion-1E/appendix-mixed-models.pdf).