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

La base de datos epilepsy del paquete HSAUR2 contiene los resultados de un ensayo clínico para comparar dos tratamientos contra la epilepsia en pacientes que sufren de la enfermedad.

  1. Consulte la ayuda de la base de datos para conocer los detalles de las variables recolectadas.

Use el siguiente código para responder las preguntas de abajo y ajustar los modelos.

data("epilepsy", package = "HSAUR2")
epilepsy$x <- as.numeric(epilepsy$period)
epilepsy$y <- epilepsy$seizure.rate
  1. ¿Cuántos pacientes tiene la base de datos?
  2. ¿Cuál fue el mayor número de seguimientos a un paciente?
  3. ¿Cuántos pacientes recibieron Progabide? ¿Cuántos pacientes recibieron el placebo?
  4. Haga un dibujo igualito al mostrado abajo para monitorear la evolución de la variable \(Y\) en función del periodo \(X\) para los pacientes 8, 14, 35, 49. Use colores para diferenciar los pacientes que recibieron Progabide y el placebo.

¿Qué conclusión puede sacar de la figura anterior?

  1. Ajuste el siguiente modelo en el objeto mod1.

\[\begin{align*} y_{ij} | b_0 &\sim Poisson(\lambda_{ij}) \\ \log(\lambda_{ij}) &= \beta_0 + \beta_1 \, x_{ij} + \beta_2 \, treatment_{Progabidei} + b_{0i} \\ b_0 &\sim N(0, \sigma^2_{b0}), \end{align*}\]

donde \(i\) corresponde al paciente y \(j\) a la visita. Note que el tratamiento de referencia es placebo. Use nAGQ=1 para la aproximación por cuadratura de la intergral interna.

  1. Ajuste el siguiente modelo en el objeto mod2.

\[\begin{align*} y_{ij} | b_0 &\sim Poisson(\lambda_{ij}) \\ \log(\lambda_{ij}) &= \beta_0 + \beta_1 \, x_{ij} + \beta_2 \, treatment_{Progabidei} + \beta_3 \, x_{ij} \, treatment_{Progabidei} + b_{0i} \\ b_0 &\sim N(0, \sigma^2_{b0}), \end{align*}\]

Use nAGQ=1 para la aproximación por cuadratura de la intergral interna.

  1. Ajuste el siguiente modelo en el objeto mod3.

\[\begin{align*} y_{ij} | {b_0, b_1} &\sim Poisson(\lambda_{ij}) \\ \log(\lambda_{ij}) &= \beta_0 + \beta_1 \, x_{ij} + \beta_2 \, treatment_{Progabidei} + \beta_3 \, x_{ij} \, treatment_{Progabidei} + b_{0i} + b_{1i} \,x_{ij} \\ \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} & \sigma_{b01} \\ \sigma_{b01} & \sigma^2_{b1} \end{matrix} \right ) \right ) \\ \end{align*}\]

Use nAGQ=1 porque cuando hay más de un efecto aleatorio se debe usar solo un punto de cuadratura.

  1. Use una prueba de razón de verosimilitud sencilla para determinar cuál de los modelos es mejor. Use el siguiente código.
anova(mod1, mod2, mod3)

¿Cuál de los dos modelos es mejor?

  1. Explore la ayuda de la función predict.merMod (o predict). ¿Para qué sirve el argumento type?

  2. Para cada uno de los modelos ajustados calcule el coeficiente de correlación de Pearson entre \(y_{ij}\) y \(\hat{\lambda}_{ij}\), es decir, \(\rho(y_{ij}, \hat{\lambda}_{ij})\).

¿Cuál modelo tiene el mayor \(\rho\)?

  1. Escriba el modelo ajustado que resultó ser el mejor modelo.

  2. Usando el modelo anterior, escriba el modelo ajustado para el paciente #8.

  3. Vuelva a construir la figura anterior con los pacientes 8, 14, 35, 49 pero incluyendo la curva ajustada de \(\hat{\lambda}\) con cada modelo. Abajo se muestra una de las tres figuras.

¿A cuál modelo pertence la figura? ¿mod1, mod2 o mod3?