Agregar tarea 10

tareas/tarea-10/datos/hearttransplants.csv

@@ -0,0 +1,95 @@
+e,y
+532,0
+584,0
+672,2
+722,1
+904,1
+1236,0
+950,0
+1405,1
+776,3
+1013,0
+739,0
+1770,1
+821,0
+1115,2
+1164,3
+1164,0
+1303,0
+1774,3
+3585,1
+1193,1
+1213,1
+1232,1
+1517,4
+1520,3
+1862,3
+1888,1
+1247,0
+1381,2
+1643,2
+1660,4
+1827,4
+1486,3
+1593,2
+2265,4
+1524,1
+1759,3
+1309,0
+1529,4
+1677,1
+1654,2
+1785,3
+1979,4
+1767,4
+2465,2
+1750,2
+2458,4
+2383,2
+2717,3
+2282,0
+2115,0
+2852,2
+2856,5
+3174,5
+2369,1
+2557,1
+3859,3
+2641,1
+2741,1
+3055,3
+3513,1
+2728,2
+3354,6
+3814,0
+4014,2
+2612,2
+2815,1
+4294,2
+3450,8
+3628,6
+4219,1
+3932,6
+4082,4
+4203,1
+4022,3
+4636,5
+5571,2
+6436,4
+5344,3
+4445,4
+4705,5
+5039,2
+6043,6
+5121,8
+11260,5
+5789,0
+6044,6
+5569,8
+6130,7
+6249,3
+7002,3
+7851,9
+9573,7
+12050,18
+12131,17

tareas/tarea-10/heart-agregado.stan

@@ -0,0 +1,30 @@
+data {
+  int<lower=0> N;
+  array[N]  int e;
+  array[N]  int y;
+}
+
+parameters {
+  real <lower=0> lambda;
+}
+
+transformed parameters {
+  vector[N] media_hospital;
+  // lambda es por cada 1000 expuestos:
+  for (i in 1:N){
+    media_hospital[i] = lambda * e[i] / 1000;
+  }
+}
+
+model {
+  // partes no determinísticas
+  y ~ poisson(media_hospital);
+  lambda ~ exponential(1);
+}
+
+generated quantities {
+  array[N] int y_sim;
+  for (i in 1:N){
+    y_sim[i] = poisson_rng(media_hospital[i]);
+  }
+}

tareas/tarea-10/heart-individual.stan

@@ -0,0 +1,30 @@
+data {
+  int<lower=0> N;
+  array[N]  int e;
+  array[N]  int y;
+}
+
+parameters {
+  vector<lower=0>[N] lambda;
+}
+
+transformed parameters {
+  vector[N] media_hospital;
+  // lambda es por cada 1000 expuestos:
+  for (i in 1:N){
+    media_hospital[i] = lambda[i] * e[i] / 1000;
+  }
+}
+
+model {
+  // partes no determinísticas
+  y ~ poisson(media_hospital);
+  lambda ~ exponential(1);
+}
+
+generated quantities {
+  array[N] int y_sim;
+  for (i in 1:N){
+    y_sim[i] = poisson_rng(media_hospital[i]);
+  }
+}

tareas/tarea-10/tarea-10.qmd

@@ -0,0 +1,154 @@
+---
+title: "Tarea 10: intro a modelos jerárquicos"
+format: html
+---
+
+```{r}
+library(tidyverse)
+library(kableExtra)
+library(DiagrammeR)
+ggplot2::theme_set(ggplot2::theme_light())
+library(cmdstanr)
+set_cmdstan_path(".cmdstan/cmdstan-2.34.1/")
+```
+
+
+
+Consideramos un ejemplo simple, donde queremos estimar el efecto del 
+hospital en la tasa de mortalidad de pacientes de cirugía de corazón. Este ejemplo
+se puede encontrar en @albert2009bayesian. Plantearemos  3 alternativas de modelación para resolver el problema: modelo de unidades iguales, modelo de unidades independientes y finalmente modelo jerárquico.
+
+En esta tarea sólo consideramos los primeros dos como preparación:
+
+- Modelo donde estimamos una sola tasa común para todos los hospitales. ¿Ajusta este modelo?
+- Modelo donde estimamos una tasa individual para cada hospital, sin tomar en cuenta cómo se comporta el conjunto de los hospitales.
+
+
+
+Tenemos datos todas las cirugías de transplante de corazón llevadas a cabo en Estados Unidos en un periodo de 24 meses, entre octubre de  1987 y diciembre de 1989. Para cada uno de los  131 hospitales, se registró el número de cirugías de transplante de corazón, y el número de muertes durante los 30 días posteriores a la cirugía $y$.
+Además, se cuenta con una una cantidad $e$ de pacientes expuestos por la cirugía
+(nota: ver detalles en las notas).
+
+El diagrama simple que consideraremos es uno donde hospital es causa tanto de 
+su exposición $e$ (por su tamaño), como del número
+de personas fallecidas. A su vez, la exposición $e$ es causa del número de muertes $y$.
+Nos interesa estimar el efecto directo de hospital en el número de muertes. 
+
+
+
+```{r}
+#| message: false
+datos_hosp <- read_csv("./datos/hearttransplants.csv") |> 
+  mutate(hospital = row_number())
+head(datos_hosp)
+```
+
+Consideramos la cantidad $y/e$ como una estimación cruda de la tasa de mortalidad.
+En la siguiente gráfica, observamos que parece ser la variabilidad es alta
+cuando el número de expuestos es relativamente baja. Nótese que
+la tasa de mortalidad no es muy alta en general, y que el número de muertes
+es relativamente bajo en muchos hospitales (puede tomar valores 0, 1, 2, etc.) Esto
+produce variabilidad alta para exposiciones bajas.
+
+
+```{r}
+ggplot(datos_hosp, aes(x = e, y = 1000 * y / e, color = log(1 + y))) +
+  geom_point() + scale_x_log10() + xlab("Número de expuestos e")
+```
+
+## Modelo agregado
+
+Consideramos primero un modelo donde consideramos que todos los hospitales
+tienen una misma tasa de mortalidad. Si $e_j$ es la exposición del hospital $j$ y $y_j$ el número de muertes, entonces podemos considerar un modelo de la forma
+
+$$y_j \sim \text{Poisson}(e_j \lambda),$$
+
+Es decir, el número de muertes es Poisson con valor esperado igual al número
+de expuestos multiplicado por la tasa común de mortalidad. 
+
+```{r}
+#| message: false
+mod_agregado <- cmdstan_model("heart-agregado.stan")
+print(mod_agregado)
+datos_agregado <- list(N = nrow(datos_hosp), y = datos_hosp$y, e = datos_hosp$e)
+ajuste_agregado <- mod_agregado$sample(data = datos_agregado, chains = 4, refresh = 1000)
+```
+
+**Pregunta 1**: checa el código de Stan y describe en tus palabras qué cantidades
+estima.
+
+```{r}
+ajuste_agregado$summary("lambda")
+```
+
+Los diagnósticos básicos de convergencia parecen ser apropiados. Procedemos a hacer un
+chequeo predictivo posterior:
+
+```{r}
+#| warning: false
+#| message: false
+set.seed(912)
+ajuste_agregado$draws("y_sim", format = "df") |> 
+  as_tibble() |> 
+  pivot_longer(cols = starts_with("y_sim"), names_to = "variable") |> 
+  separate(variable, into = c("variable", "hospital"), sep = "[\\[\\]]") |>
+  mutate(hospital = as.integer(hospital)) |>
+  left_join(datos_hosp, by = "hospital") |>
+  filter(hospital %in% sample(1:94, 20)) |>
+  ggplot(aes(x = value)) + geom_histogram(binwidth = 1) +
+  facet_wrap(~ hospital) + 
+  geom_vline(aes(xintercept = y), color = "red")
+```
+
+Y vemos fallas en el ajuste del modelo, con varias observaciones
+en los extremos de las colas.
+
+
+## Modelo individual
+
+Podemos considerar un modelo donde cada hospital tiene su propia tasa de mortalidad.
+
+
+```{r}
+library(cmdstanr)
+mod_ind <- cmdstan_model("heart-individual.stan")
+print(mod_ind)
+datos_ind <- list(N = nrow(datos_hosp), y = datos_hosp$y, e = datos_hosp$e)
+ajuste_ind <- mod_ind$sample(data = datos_ind, chains = 4, refresh = 1000)
+resumen <- ajuste_ind$summary("lambda") |> 
+  select(variable, mean, sd, rhat, ess_bulk)
+resumen |> kable()
+```
+
+
+**Pregunta 2**: checa el código de Stan y describe en tus palabras qué cantidades
+estima.
+
+El problema en este caso es que tenemos intervalos que simplemente no
+son creíbles, en particular con aquellos hospitales que tienen poca exposición. 
+
+```{r}
+#| message: false
+#| warning: false
+set.seed(912)
+ajuste_ind$draws("lambda", format = "df") |> 
+  as_tibble() |> 
+  pivot_longer(cols = starts_with("lambda"), names_to = "variable") |> 
+  separate(variable, into = c("variable", "hospital"), sep = "[\\[\\]]") |>
+  mutate(hospital = as.integer(hospital)) |>
+  left_join(datos_hosp, by = "hospital") |>
+  mutate(hospital = factor(hospital)) |>
+  group_by(hospital, e, y) |> 
+  summarise(inf = quantile(value, 0.1), sup = quantile(value, 0.9)) |>
+  ggplot(aes(x = e)) + geom_linerange(aes(ymin = inf, ymax = sup)) +
+  geom_point(aes(y = 1000 * y / e), color = "red") +
+  scale_x_log10() + xlab("Número de expuestos e") + ylab("Muertes por mil expuestos")
+```
+
+En este caso, la variabilidad es muy alta para hospitales con poca exposición, tanto
+en los datos observados como en los intervalos. Los intervalos no aportan
+mucha información. En este punto utilizar iniciales fuertes
+para las $\lambda_j$ si tenemos la información disponible. Sin embargo, los
+resultados serán altamente sensible a esta información inicial.
+
+