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.
\[\begin{align*} y_{ij} &\sim N(\mu_{ij}, \sigma^2_y) \\ \mu_{ij} &= 4 - 5 x_{ij} + b_{0i} + b_{1i} x_{ij} \\ \sigma^2_y &= 4 \\ \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}=40 & \sigma_{b01}=3 \\ \sigma_{b01}=3 & \sigma^2_{b1}=50 \end{matrix} \right ) \right ) \\ x_{ij} &\sim U(0, 5) \end{align*}\]
Escriba los 6 elementos del vector de parámetros \(\boldsymbol{\Theta}\) del modelo.
&&& por el código correcto para simular los datos solicitados.ni <- &&&
G <- &&&
nobs <- &&& * &&& # Numero total de observaciones
grupo <- factor(rep(x=1:G, &&&=ni)) # Para crear la variable grupal
obs <- rep(x=1:ni, &&&=G) # Para identificar las obs por grupo
x <- &&& # La covariable
library(MASS) # Libreria para simular obs de Normal bivariada
Sigma <- matrix(c(40, 3, # Matriz de var-cov
3, 50), ncol=2, nrow=2)
b <- mvrnorm(n=G, mu=rep(0, 2), Sigma=Sigma) # Simulando b0 y b1
b <- apply(b, MARGIN=2, function(c) rep(c, each=ni)) # Replicando
b0 <- as.vector(b[, 1]) # Separando los b0
b1 <- as.vector(b[, 2]) # Separando los b1
media <- &&& + &&& * x + b0 + b1 * x # La media
y <- rnorm(n=&&&, mean=&&&, sd=&&&) # La variable respuesta
datos <- &&&(grupo, obs, &&&, &&&, &&&, &&&) # 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.
¿Logra ver en la figura el efecto de \(b_0\) y \(b_1\) para inducir interceptos y pendientes aleatorias? Si no logra verlo con claridad, vuelva a simular un nuevo conjunto de datos.
En este ejemplo con datos artificiales (simulados) se usaron valores grandes para las varianzas para obtener nubes con el patrón esperado. En los casos prácticos esos valores de varianzas no son tan grandes.
Encuentre las estimaciones \(\hat{\boldsymbol{\Theta}}\) del vector \(\boldsymbol{\Theta}\). ¿Qué tan cerca está su estimación del vector real?
Encuentre las predicciones de los \(b_0\) y \(b_1\), luego compárelas con los valores verdaderos.
Pregunte en stackoverflow sobre cómo generar datos para el mismo modelo considerado aquí.