Este taller está basado en los ejemplos del libro Data Analysis Using Regression and Multilevel/Hierarchical Models. Los autores del libro dispusieron los datos y los scripts de R en el siguiente enlace http://www.stat.columbia.edu/~gelman/arm/. Visite el enlace para descargar el material necesario.

Preguntas

  1. Simular \(n=5\) observaciones del siguiente modelo normal:

\[\begin{align*} y_i &\sim N(\mu_i, \sigma^2) \\ \mu_i &= 4 - 6 x_i \\ x_i &\sim U(-5, 6) \\ \sigma^2 &= 16 \end{align*}\]

  1. Considere el siguiente modelo normal mixto con intercepto aleatorio \(b_{0i}\).

\[\begin{align*} y_{ij} &\sim N(\mu_{ij}, \sigma^2) \\ \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*}\]

Complete el siguiente código de R para simular \(n_i=50\) observaciones para \(G=10\) grupos, es decir, en total 500 observaciones. Cambie los símbolos &&& por el código correcto para simular los datos solicitados.

ni <- 50
G <- &&&
nobs <- ni * &&&                      # Numero total de observaciones
grupo <- factor(rep(x=1:G, each=&&&)) # Para crear la variable grupal
obs <- rep(x=1:ni, times=G)           # Para identificar las obs por grupo
x <- runif(n=nobs, min=&&&, max=&&&)  # La covariable
b0 <- rnorm(n=G, mean=&&&, sd=&&&)    # El Intercepto aleatorio
b0 <- rep(x=b0, each=ni)              # El intercepto aleatorio pero repetido
media <- &&& - 6 * x + &&&            # La media
&&& <- rnorm(n=nobs, mean=media, sd=&&&) # La variable respuesta
datos <- data.frame(grupo, obs, b0, x, &&&) # Organizando el dataframe

El objeto datos contendrá los datos simulados. Imprima en pantalla los datos y mire explore las variables, encuéntrele sentido a los números.

  1. Use el código de abajo para dibujar los datos simulados, el resultado que usted va a obtener es una figura similar (no igual porque estamos simulando sin fijar la semilla) a la figura que se muestra luego.
library(ggplot2)
ggplot(datos, aes(x, y, color=grupo) ) + 
  geom_point() + 
  labs(colour="Grupo/Cluster")

  1. ¿Qué observa en la figura? ¿Hay algún patrón evidente?

  2. ¿Qué efecto tiene la inclusión del intercept \(b_0\) en la simulación de los datos?

  3. Escriba los 4 elementos del vector de parámetros del modelo con intercepto aleatorio.

  4. Complete el siguiente código para estimar el vector de parámetros \(\boldsymbol{\Theta}\).

library(lme4)
fit <- lmer(&&& ~ x + (1 | &&&), data=&&&)
summary(fit)

El resultado de usar summary(fit) será algo similar a lo que se muestra en la siguiente figura.

Según la figura anterior,

\[\hat{\Theta}=(\hat{\beta}_0=3.53363, \hat{\beta}_1=-5.97805, \hat{\sigma}_y=4.012, \hat{\sigma}_{bo}=21.507)^\top\]

  1. ¿Cómo extraer los efectos fijos?
library(lme4)
## Loading required package: Matrix
fit <- lmer(y ~ x + (1 | grupo), data=datos)

Muy sencillo, solo debe usar la función fixed sobre el modelo ajustado fit así:

fixef(fit)
## (Intercept)           x 
##    3.533631   -5.978053
  1. ¿Cómo extraer las predicciones de interceptos aleatorios?

Muy sencillo, para obtener las predicciones \(\tilde{b}_0\) debe usar la función ranef sobre el modelo ajustado fit así:

ranef(fit)
## $grupo
##    (Intercept)
## 1   -26.277747
## 2   -30.861272
## 3    26.285385
## 4    13.720291
## 5   -20.037461
## 6     8.666852
## 7    20.551113
## 8   -17.733851
## 9    19.722390
## 10    5.964300
## 
## with conditional variances for "grupo"
  1. ¿Qué tan buenas son esas predicciones \(\tilde{b}_0\) de los verdaderos \(b_0\)?

Sencillo, podemos comparar los resultados construyendo una tabla, en la primera columna el verdadero \(b_0\) usado y en la segunda columna el valor predicho \(\tilde{b}_0\), el código podría ser:

cbind(True=unique(b0), ranef(fit)$grupo)
##          True (Intercept)
## 1  -27.216339  -26.277747
## 2  -31.005875  -30.861272
## 3   26.194197   26.285385
## 4   13.279302   13.720291
## 5  -20.716457  -20.037461
## 6    7.552105    8.666852
## 7   20.107027   20.551113
## 8  -18.690492  -17.733851
## 9   19.813066   19.722390
## 10   5.624403    5.964300

¿Qué opina usted de los \(\tilde{b}_0\)?

Nota: esta comparación solo es posible hacerla cuando tenemos datos simulados y disponemos de los b0, en un caso de aplicación con datos reales los b0 son desconocidos y no se podrán comparar como se hizo aquí.

  1. Creando el modelo ajustado para cada grupo.

Recordemos que los efectos fijos estimados fueron \(\hat{\beta}_0 \approx 3.53\) y \(\hat{\beta}_1\approx -5.98\). Para el grupo #1 se obtuvo \(\tilde{b}_{0, i=1}=-26.27\), así que la media del grupo #1 se calcula así:

\[\begin{align*} \hat{\mu}_{1j} &= \hat{\beta}_0 + \hat{\beta}_0 x_{1j} + \tilde{b}_{01} \\ \hat{\mu}_{1j} &= 3.53 -5.98 x_{1j} - 26.27 \\ \hat{\mu}_{1j} &= -22.74 -5.98 x_{1j} \end{align*}\]

Lo anterior se puede resumir en el siguiente modelo.

\[\begin{align*} y_{1j} &\sim N(\hat{\mu}_{1j}, \hat{\sigma}^2_y) \\ \hat{\mu}_{1j} &= -22.74 -5.98 x_{1j} \\ \hat{\sigma}_y &= 4.012 \end{align*}\]

Calcule usted ahora el modelo para las observaciones del grupo #3.

  1. ¿Cómo extraer los elementos para construir los modelos individuales?

Puede usar el siguiente código.

coef(fit)
## $grupo
##    (Intercept)         x
## 1   -22.744116 -5.978053
## 2   -27.327641 -5.978053
## 3    29.819016 -5.978053
## 4    17.253922 -5.978053
## 5   -16.503830 -5.978053
## 6    12.200482 -5.978053
## 7    24.084744 -5.978053
## 8   -14.200221 -5.978053
## 9    23.256020 -5.978053
## 10    9.497931 -5.978053
## 
## attr(,"class")
## [1] "coef.mer"
  1. Predicción.

Considere el modelo ajustado para el grupo #1. ¿Cuál será el valor ajustado de la media \(\hat{\mu}_{1j}\) cuando \(x=3\).

  1. ¿Cuál es la clase del objeto fit?
class(fit)
## [1] "lmerMod"
## attr(,"package")
## [1] "lme4"
  1. Predicción con R

¿Hay alguna forma de obtener predicciones para objetos de la clase lmerMod? Consúltelo!!!