-
Notifications
You must be signed in to change notification settings - Fork 6
/
Copy path701-R2.Rmd
190 lines (134 loc) · 7.38 KB
/
701-R2.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
# Coeficiente $R^2$ {#R2}
En este capítulo vamos mostrar como se puede calcular la medida $R^2$ propuesta por @nakagawa2012 y @nakagawa2017 para modelos lineales generalizados mixtos.
La función `r.squaredGLMM()` del paquete **MuMIn** de @R-MuMIn permite calcular $R^2$ para varios tipos de modelos con efectos mixtos. La función puede calcular $R^2$ para modelos de la clase merMod, lme, glmmTMB, glmmADMB, glmmPQL, cpglm(m) y
(g)lm.
@nakagawa2017 define dos $R^2$ para LMM (linear mixed models) de la siguiente manera:
- $R^2_{m}$: representa la varianza explicada por los efectos fijos.
- $R^2_{c}$: representa la varianza explicada por todo el modelo, incluyendo ambos efectos fijos y aleatorios.
Las expresiones para obtener ambos $R^2$ son:
$$
R^2_{m} = \frac{\sigma^2_f}{\sigma^2_f + \sigma^2_\alpha + \sigma^2_\epsilon},
$$
$$
R^2_{c} = \frac{\sigma^2_f + \sigma^2_\alpha}{\sigma^2_f + \sigma^2_\alpha + \sigma^2_\epsilon},
$$
donde $\sigma^2_f$ es la varianza asociada a los efectos fijos, $\sigma^2_\alpha$ es la varianza de los efectos aleatorios y $\sigma^2_\epsilon$ es la varianza asociada a los valores de la variable respuesta.
@nakagawa2017 en su artículo proponen el cálculo de $R^2$ para otros GLMM (generalized linear mixed models). Los autores afirman que el reto al proponer un $R^2$ para los GLMM radica en la definición de $\sigma^2_\epsilon$. En la implementación de $R^2$ para GLMM hecha por @R-MuMIn se pueden obtener tres formas para cuantificar $\sigma^2_\epsilon$, esas formas son: “delta”, “lognormal” y “trigamma”.
En los siguientes ejemplos vamos a ilustrar el uso de la función `r.squaredGLMM()`.
## 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 y en la cual la varianza de los efectos aleatorio $b_0$ ($\sigma^2_{b0}=625$) es más grande que la varianza de las observaciones ($\sigma^2_y=16$).
El objetivo del ejemplo es calcular el $R^2$ para el modelo ajustado correctamente.
El modelo para simular los datos es el siguiente.
\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.
```{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)
x <- runif(n=nobs, min=-5, max=6)
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
y <- rnorm(n=nobs, mean=media, sd=sqrt(16))
datos <- data.frame(obs, grupo, b0, x, media, y)
```
Ahora vamos a ajustar el modelo correcto así:
```{r, message=FALSE}
library(nlme)
fit1 <- lme(y ~ x, random = ~ 1 | grupo, data=datos)
```
A continuación vamos a calcular los $R^2$.
```{r}
library(MuMIn)
r.squaredGLMM(fit1)
```
En la salida anterior se reporta la medida marginal y la medida condicional del $R^2$. De los resultados se observa que el $R^2_{m}$ es mayor que el $R^2_{c}$, eso significa que el modelo completo (con efectos fijos y aleatorio) logra explicar un 98\% de la varianza. En otras palabras, significa que si valió la pena ajustar el modelo con esos efectos aleatorios. Era de esperarse que esto sucediera porque los datos se simularon con $\sigma^2_{b0}=625$.
## Ejemplo: modelo normal sin intercepto aleatorio {-}
En este ejemplo vamos a repetir a simular datos como el ejemplo anterior pero con $\sigma^2_{b0}=0$. La idea es ver si el $R^2$ logra detectar que no se necesita incluir un intercepto aleatorio al ajustar el modelo.
<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.
```{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)
x <- runif(n=nobs, min=-5, max=6)
b0 <- rnorm(n=G, mean=0, sd=0) # Intercepto aleatorio
b0 <- rep(x=b0, each=ni) # El mismo intercepto aleatorio pero repetido
media <- 4 - 6 * x + b0
y <- rnorm(n=nobs, mean=media, sd=sqrt(16))
datos <- data.frame(obs, grupo, b0, x, media, y)
```
Ahora vamos a ajustar el modelo así:
```{r, message=FALSE}
library(nlme)
fit2 <- lme(y ~ x, random = ~ 1 | grupo, data=datos)
```
A continuación vamos a calcular los $R^2$.
```{r}
library(MuMIn)
r.squaredGLMM(fit2)
```
De la salida anterior vemos que ambos $R^2$ son muy parecidos. Esto significa que el porcentaje de la varianza explicada por un modelo con ambos tipos de efectos (fijos y aleatorios) es igual a la varianza explicada por un modelo sólo con efectos fijos, en otras palabras, esto indica que NO vale la pena incluir efectos aleatorios en el modelo. Este resultado era el esperado porque los datos fueron simulados con $\sigma^2_{b0}=0$.
## Ejemplo: modelo clase lme {-}
En este ejemplo vamos a ajustar un modelo para explicar la desde la hipófisis hasta la fisura pterigomaxilar (mm), en función del sexo y la edad de los pacientes (Subject). La base de datos a usar es `Orthodont`. Luego se desea calcular el $R^2$ para el modelo ajustado.
<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 que vamos a usar es el siguiente.
```{r}
library(nlme)
library(MuMIn)
data(Orthodont, package = "nlme")
fm1 <- lme(distance ~ Sex * age, ~ 1 | Subject, data = Orthodont)
r.squaredGLMM(fm1)
```
De la salida anterior se observa que 78.27\% de la variabilidad observada se logra explicar con el modelo ajustado con efectos fijos y un intercepto aleatorio.
## Ejemplo: modelo clase lme4 {-}
En este ejemplo vamos a calcular $R^2$ para el modelo gamma presentado en el ejemplo 1 sobre conductores que está en el capítulo \@ref(glmm-gamma).
<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>
A continuación se presenta el mismo código usado en el ejemplo 1 del capítulo \@ref(glmm-gamma). En la última línea de este código se calcula el $R^2$.
```{r message=FALSE}
library(hglm)
data(semiconductor)
library(lme4)
mod1 <- glmer(y ~ x1 + x2 + x3 + x4 + x5 + x6 + (1 | Device),
data = semiconductor,
family = Gamma(link = log))
library(MuMIn)
r.squaredGLMM(mod1)
```
De la salida anterior se observa que el $R^2_c$ es ligeramente mayor al $R^2_m$. Los valores de $R^2_c$ son alrededor de 77\%, esto indica que el modelo `mod1` con efectos fijos y un intercepto aletorio, logra explicar un 77\% de la variabilidad observada. Este resultado es una evidencia en favor del modelo mixto frente al modelo sólo con efectos fijos.