-
Notifications
You must be signed in to change notification settings - Fork 6
/
Copy path401-paquete-glmmTMB.Rmd
237 lines (179 loc) · 8.69 KB
/
401-paquete-glmmTMB.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
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
# Paquete **glmmTMB** {#pac-glmmTMB}
El paquete **glmmTMB** de @R-glmmTMB se utiliza para estimar modelos glmm por medio de máxima verosimilitud a través de 'TMB' (Template Model Builder). Se supone que los efectos aleatorios son gaussianos en la escala del predictor lineal y se integran utilizando la aproximación de Laplace. Los gradientes se calculan mediante la diferenciación automática. Al visitar este [enlace](https://cran.r-project.org/web/packages/glmmTMB/) se encontrará la página de apoyo del paquete, allí se puede consultar el manual de referencia y las viñetas.
## Función `glmmTMB`
La función `glmmTMB` es la principal función del paquete **glmmTMB**. Esta función sirve para ajustar un glmm y su estructura es la siguiente:
```{r, eval=FALSE}
glmmTMB(formula,
data = NULL,
family = gaussian(),
ziformula = ~0,
dispformula = ~1,
weights = NULL,
offset = NULL,
contrasts = NULL,
na.action,
se = TRUE,
verbose = FALSE,
doFit = TRUE,
control = glmmTMBControl(),
REML = FALSE,
start = NULL,
map = NULL,
sparseX = NULL
)
```
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.
- `family`: argumento para indicar la distribución de la variable respuesta.
- `REML`: valor lógico que sirve para indicar si queremos estimaciones maximizando la [verosimilitud restringida](https://en.wikipedia.org/wiki/Restricted_maximum_likelihood) o la verosimilitud usual.
## 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 `glmmTMB` 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$.
<div style="-moz-box-shadow: 1px 1px 3px 2px #000000;
-webkit-box-shadow: 1px 1px 3px 2px #000000;
box-shadow: 1px 1px 3px 2px #000000;">
```{block2, type='rmdexercise'}
Solución.
```
</div>
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.
```{r, eval=TRUE, echo=TRUE}
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.
```{r}
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.
```{r}
library(ggplot2)
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 modelo se usa la función `glmmTMB` de la siguiente forma.
```{r, message = FALSE, eval = TRUE, echo = TRUE}
library(glmmTMB)
fit1 <- glmmTMB(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`.
```{r, eval=TRUE}
summary(fit1)
```
Según el resultado anterior $\hat{\boldsymbol{\Theta}}=(\hat{\beta}_0=2.2378, \hat{\beta}_1=-6.0264, \hat{\sigma}_y=3.931, \hat{\sigma}_{b0}=24.066)^\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$.
<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'}
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 \@ref(pac-nlme).
¿Hay alguna similitud?
```
</div>
## Ejemplo: recuperando los interceptos aleatorios {-}
¿Cómo se pueden obtener los interceptos aleatorios a partir del modelo ajustado en la sección anterior?
<div style="-moz-box-shadow: 1px 1px 3px 2px #000000;
-webkit-box-shadow: 1px 1px 3px 2px #000000;
box-shadow: 1px 1px 3px 2px #000000;">
```{block2, type='rmdexercise'}
Solución.
```
</div>
Para obtener los interceptos aleatorios se usa la función `ranef` del paquete **glmmTMB** de @R-glmmTMB. A continuación vamos a obtener los interceptos aleatorios y los vamos a comparar con los $b_0$ simulados.
```{r}
interceptos_aleatorios <- ranef(fit1)
cbind(interceptos_aleatorios$grupo, b0=unique(b0))
```
De la salida anterior vemos que los $\tilde{b}_0$ son cercanos a los valores reales de $b_0$.
<div style="-moz-box-shadow: 1px 1px 3px 2px #808080;
-webkit-box-shadow: 1px 1px 3px 2px #808080;
box-shadow: 1px 1px 3px 2px #808080;">
```{block2, type='rmdwarning'}
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$.
```
</div>
## 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 `glmmTMB` 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 &\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$.
<div style="-moz-box-shadow: 1px 1px 3px 2px #000000;
-webkit-box-shadow: 1px 1px 3px 2px #000000;
box-shadow: 1px 1px 3px 2px #000000;">
```{block2, type='rmdexercise'}
Solución.
```
</div>
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.
```{r}
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.
```{r}
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.
```{r}
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.
```{r}
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.
```{r, message = FALSE, eval = TRUE, echo = TRUE}
library(glmmTMB)
fit2 <- glmmTMB(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 la salida de `summary`.
```{r}
summary(fit2)
```
Según el resultado anterior $\hat{\boldsymbol{\Theta}}=(\hat{\beta}_0=4.4102, \hat{\beta}_1=-8.0789, \hat{\phi}=0.507, \hat{\sigma}_{bo}=3.012)^\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$.