-
Notifications
You must be signed in to change notification settings - Fork 6
/
501-glmm-Poisson.Rmd
145 lines (106 loc) · 5.46 KB
/
501-glmm-Poisson.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
# GLMM Poisson {#glmm_poisson}
En este capítulo se presenta un ejemplo de glmm con variable respuesta Poisson y está basado en [esta publicación](https://stats.oarc.ucla.edu/r/faq/random-coefficient-poisson-models/).
```{r echo=FALSE, out.width="80%"}
knitr::include_graphics("images/premios_estudiantes.png")
```
A continuación la base de datos a utilizar.
```{r, message = FALSE}
require(foreign)
datos <- read.dta("https://stats.idre.ucla.edu/stat/data/hsbdemo.dta")
# Cambios menores a algunas variables
datos$cid <- factor(datos$cid)
datos$female <- ifelse(datos$female == 'female', 1, 0)
```
Los datos corresponden a una muestra de 200 estudiantes de educación secundaria, los cuales están agrupados en 20 escuelas diferentes. El propósito es modelar el número de premios que recibe un estudiante dada la escuela a la cual pertenece, teniendo en cuenta la variable género como la variable que mejor podría ayudar en la predicción.
Las variables de la base de datos son:
- `awards`: Variable respuesta que hace referencia al número de premios que recibe un estudiante.
- `id`: Identificación única de los estudiantes de educación secundaria que participaron en el estudio.
- `female`: Variable correspondiente al género de los estudiantes.
- `ses`: Variable categórica que representa el estatus socioeconómico del estudiante, teniendo tres distinciones las cuales son: bajo, medio y alto.
- `schtyp`: Tipo de escuela, ya sea privado o público.
- `prog`: Programa de formación del estudiante, ya sea general, vocación o académico.
- `read`: Puntuación del estudiante en comprensión lectora.
- `write`: Puntuación del estudiante en escritura.
- `math`: Puntuación del estudiante en matemáticas.
- `science`: Puntuación del estudiante en ciencias.
- `honors`: Variable que hace referencia si el estudiante presenta matrículas de honor o no.
- `cid`: Variable que indica la escuela a la cual pertenece el estudiante.
Vamos a explorar las primeras líneas de la base de datos.
```{r}
head(datos, n=5)
```
El siguiente código sirve para construir un histograma que muestra el número de premios por cada escuela.
```{r message=FALSE}
library(ggplot2)
ggplot(datos, aes(awards)) +
geom_histogram(binwidth = 0.5) +
facet_wrap(~cid)
```
De la figura anterior vemos que el comportamiento de la variable respuesta `awards` es muy diferente dada la escuela.
En la siguiente figura se relaciona el número de premios `awards` con la variable género.
```{r}
ggplot(datos, aes(factor(awards))) +
geom_bar(aes(fill = factor(female)),
position = "fill") +
geom_hline(yintercept = 0.5) +
ylab("Proporción") + xlab("Número de premios")
```
De la figura anterior se observa una posible relación entre el número de premios está relacionado y el género.
Los modelos lineales generalizados mixtos nos permiten modelar la variable respuesta con la distribución Poisson o la binomial negativa, en este caso, vamos a usar la distribución Poisson.
El primer modelo que vamos a considerar aquí es el siguiente:
\begin{align*}
y_{ij} | b_0 &\sim Poisson(\mu_{ij}) \\
\log(\mu_{ij}) &= \beta_0 + b_{0i} \\
b_0 &\sim N(0, \sigma_{b0}^2)
\end{align*}
con $i=1, 2, \ldots, 20$ y $j=1, 2, \ldots, n_i$.
El modelo anterior se va a ajustar con la función `glmer` del paquete **lme4** de @R-lme4 y utilizando 15 puntos en la aproximación de Gauss para la log-verosimilitud.
```{r message=FALSE}
require(lme4)
m1 <- glmer(awards ~ 1 + (1 | cid), data = datos,
family = poisson(link = "log"), nAGQ = 15)
```
Los resultados del modelo se muestran a continuación.
```{r}
summary(m1)
```
El segundo modelo que se propone es una modificación del modelo 1, agregando como variable explicativa el género femenino. El modelo propuesto es el siguiente:
\begin{align*}
y_{ij} | b_0 &\sim Poisson(\mu_{ij}) \\
\log(\mu_{ij}) &= \beta_0 + \beta_1 female_{ij} + b_{0i} \\
b_0 &\sim N(0, \sigma_{b0}^2)
\end{align*}
con $i=1, 2, \ldots, 20$ y $j=1, 2, \ldots, n_i$.
```{r}
require(lme4)
m2 <- glmer(awards ~ 1 + female + (1 | cid), data = datos,
family = poisson(link = "log"), nAGQ = 15)
```
Los resultados del modelo se muestran a continuación.
```{r}
summary(m2)
```
Para comparar los dos modelos ajustados propuestos colocamos los resultados en una única tabla como se muestra a continuación.
```{r message=FALSE}
library("texreg")
screenreg(list(m1, m2))
```
```{r message=FALSE, results='asis'}
library(stargazer)
stargazer(m1, m2, type = "html")
```
- En los resultados que obtenemos de esta tabla podemos ver que la log-verosimilitud aument para el modelo 2 en relación ál modelo 1.
- Además podemos observar que tanto el AIC como el BIC del modelo 2 son menores que los del modelo 1.
- Por lo tanto, basados en esto, concluimos que el mejor modelo propuesto es el modelo 2, es decir, vale la pena agregarle como variable explicativa el género femenino, va a generar un modelo más confiable.
```{r}
## QQ plot
plot(ranef(m2))
## Caterpillar plot
lattice::dotplot(ranef(m2, condVar = TRUE))
```
Los modelos anteriores se pueden ajustar usando la función `glmmadmb` del paquete *glmmADMB*, a continuación se muestra el código. Para conocer más sobre este paquete se recomienda visitar [este enlace](https://glmmadmb.r-forge.r-project.org/)
```{r eval=FALSE}
library(glmmADMB)
m1_alt <- glmmadmb(awards ~ 1 + (1 | cid), data = datos, family = "poisson", link = "log")
m2_alt <- glmmadmb(awards ~ 1 + female + (1 | cid), data = datos, family = "poisson", link = "log")
```