diff --git a/tareas/tarea-08/backdoor-logistico.stan b/tareas/tarea-08/backdoor-logistico.stan new file mode 100644 index 0000000..e9a864b --- /dev/null +++ b/tareas/tarea-08/backdoor-logistico.stan @@ -0,0 +1,70 @@ +data { + int N; + vector[N] trata; + array[N] int res; + vector[N] peso; + +} + +transformed data { + real media_peso; + + // centrar + media_peso = mean(peso); +} + +parameters { + real gamma_0; + real gamma_1; + real gamma_2; +} + +transformed parameters { + vector[N] p_logit_res; + + p_logit_res = gamma_0 + gamma_1 * trata + gamma_2 * (peso - media_peso); + +} + +model { + // modelo de resultado + // bernoulli_logit transforma p_logit_res a una probabilidad: + res ~ bernoulli_logit(p_logit_res); + gamma_0 ~ normal(0, 2); + gamma_1 ~ normal(0, 0.5); + gamma_2 ~ normal(0, 0.5); + + +} +generated quantities { + real dif_trata; + real p_trata; + real p_no_trata; + vector[N] probs; + + // Ahora hacemos simulaciones para comparar una + // población de tratados vs una de no tratados + // marginalizamos sobre peso tomando en cada simulación + // un peso de la población. + for(i in 1:N){ + probs[i] = 1.0 / N; + } + + { + array[2000] int res_trata; + array[2000] int res_no_trata; + for(k in 1:2000){ + real peso_sim = peso[categorical_rng(probs)]; + res_trata[k] = bernoulli_rng( + inv_logit(gamma_0 + gamma_1 * 1 + + gamma_2 * (peso_sim - media_peso))); + res_no_trata[k] = bernoulli_rng( + inv_logit(gamma_0 + gamma_1 * 0 + + gamma_2 * (peso_sim - media_peso))); + } + p_trata = mean(res_trata); + p_no_trata = mean(res_no_trata); + } + dif_trata = p_trata - p_no_trata; + +} \ No newline at end of file diff --git a/tareas/tarea-08/proceso-completo.stan b/tareas/tarea-08/proceso-completo.stan new file mode 100644 index 0000000..cb9deac --- /dev/null +++ b/tareas/tarea-08/proceso-completo.stan @@ -0,0 +1,98 @@ +data { + int N; + array[N] int trata; + array[N] int res; + array[N] real peso; + array[N] int nse; + +} + +transformed data { + real media_peso; + + // centrar + media_peso = mean(peso); +} + +parameters { + real gamma_0; + real gamma_1; + real gamma_2; + real delta_0; + real delta_1; + real alpha_0; + real alpha_1; + real sigma_peso; + real p_nse; +} + +transformed parameters { + vector[N] p_logit_res; + vector[N] p_logit_trata; + vector[N] mu_peso; + vector[N] peso_c; + + for(i in 1:N){ + peso_c[i] = peso[i] - mean(peso); + } + + + p_logit_res = gamma_0 + gamma_1 * to_vector(trata) + gamma_2 * peso_c; + p_logit_trata = delta_0 + delta_1 * to_vector(nse); + mu_peso = alpha_0 + alpha_1 * to_vector(nse); + +} + +model { + // modelo de resultado + // bernoulli_logit transforma p_logit_res a una probabilidad: + nse ~ bernoulli(p_nse); + peso ~ normal(mu_peso, sigma_peso); + trata ~ bernoulli_logit(p_logit_trata); + res ~ bernoulli_logit(p_logit_res); + gamma_0 ~ normal(0, 2); + gamma_1 ~ normal(0, 0.5); + delta_0 ~ normal(0, 2); + delta_1 ~ normal(0, 0.5); + gamma_2 ~ normal(0, 0.5); + alpha_0 ~ normal(0, 20); + alpha_1 ~ normal(0, 10); + sigma_peso ~ normal(0, 20); + + +} +generated quantities { + real dif_trata; + real p_trata; + real p_no_trata; + vector[N] probs; + + // Ahora hacemos simulaciones para comparar una + // población de tratados vs una de no tratados + // marginalizamos sobre peso tomando en cada simulación + // un peso de la población. + for(i in 1:N){ + probs[i] = 1.0 / N; + } + + { + array[2000] int res_trata; + array[2000] int res_no_trata; + for(k in 1:2000){ + // simulamos según el modelo gráfico + real nse_sim = bernoulli_rng(p_nse); + real peso_c_sim = normal_rng(alpha_0 + alpha_1 * nse_sim, sigma_peso); + // pero no incluimos la ecuación de tratamiento, pues lo queremos + // manipular + res_trata[k] = bernoulli_rng( + inv_logit(gamma_0 + gamma_1 * 1 + + gamma_2 * peso_c_sim)); + res_no_trata[k] = bernoulli_rng( + inv_logit(gamma_0 + gamma_1 * 0 + + gamma_2 * peso_c_sim)); + } + p_trata = mean(res_trata); + p_no_trata = mean(res_no_trata); + } + dif_trata = p_trata - p_no_trata; +} diff --git a/tareas/tarea-08/tarea-08.qmd b/tareas/tarea-08/tarea-08.qmd new file mode 100644 index 0000000..9a1b9d9 --- /dev/null +++ b/tareas/tarea-08/tarea-08.qmd @@ -0,0 +1,250 @@ +--- +title: "Tarea 8: puertas traseras" +format: html +--- + + + +En este ejercicio resolveremos el problema que vimos en clase +acerca del criterio de puertas traseras, identificaremos el efecto causal de interés +y usaremos Stan para hacer una estimación de este efecto casual. + +```{r} +library(tidyverse) +library(DiagrammeR) +library(cmdstanr) +``` + + + +El modelo gráfico que consideramos es el siguiente, donde queremos +evaluar la efectividad de un tratamiento en cierta enfermedad. +Los datos que tenemos disponibles (son datos observados de una encuesta, por +ejemplo), son si una persona +utilizó o no un tratamiento, y si se recuperó o no. +No se registró el nivel socioeconómico, pero sabemos que el tratamiento +es caro, de forma que fue accedido más por gente de NSE más alto. + +También que sabemos que para este tipo de tratamiento, el peso de la +persona es un factor importante para la recuperación, independiente +de si tomó el tratamiento o no. + +Nuestros supuestos están en la siguiente gráfica: + +```{r} +#| code-fold: true +grViz(" +digraph { + graph [ranksep = 0.2, rankdir = LR] + node [shape=plaintext] + Trata + Res + node [shape = circle] + NSE + Peso + edge [minlen = 3] + NSE -> Peso + NSE -> Trata + Trata -> Res + Peso -> Res + +} +") + +``` + + + +## Parte 1: Cantidades a estimar teóricas, simulación y aleatorización + +Primero supondremos que conocemos el proceso generador, y veremos cómo +podemos calcular el efecto causal de interés, que es: + +- Si asignáramos el tratamiento al azar en la población, ¿que efecto +tendría? Otra forma de decir esto es decir: ¿cómo se compararía una +población grande de este proceso generador que toma el tratamiento +vs una que no toma el tratamiento? + + +En este caso usamos directamente +las relaciones funcionales entre los nodos que por el momento suponemos conocidas: + +```{r} +inv_logit <- function(x) 1 / (1 + exp(-x)) +simular_datos <- function(n = 10){ + nse <- sample(c(0, 1), n, replace = TRUE) + peso <- rnorm(n, 70 - 7 * nse, 12 + 2 * nse) |> round() + trata <- rbinom(n, 1, 0.8 * nse + 0.2 * (1 - nse)) + ## aquí es donde el tratamiento acúta sobre el resultado + p_trata <- inv_logit(0.5 * trata - 0.2 * (peso - 70)) + res <- rbinom(n, 1, p_trata) + # regresar los datos + tibble(nse, peso, trata, res) +} +``` + +Podríamos con estos datos calcular cómo se relaciona el tratamiento +con el resultado, de forma observacional: + + +```{r} +datos_sim <- simular_datos(5000) +head(datos_sim) +``` + +Nótese, que por lo que sabemos del proceso generador, es incorrecto hacer +simplemente (no contesta la pregunta de interés): + +```{r} + datos_sim |> + group_by(trata) |> + summarise(p_res = mean(res)) |> + pivot_wider(names_from = trata, values_from = p_res) |> + mutate(dif_p_trata = `1` - `0`) +``` + +Sin embargo, como conocemos todos los detalles el proceso generador, podemos +simular el mismo proceso, pero suponiendo que asignamos al azar +el tratamiento. Lo único que tenemos que hacer que es que tratamiento +**no** dependa de nse en nuestra gráfica. Como conocemos todos los detalles, +nuestro nuevo proceso generador modificado por esta intervención es: + + +```{r} +simular_intervencion <- function(n = 10){ + nse <- sample(c(0, 1), n, replace = TRUE) + peso <- rnorm(n, 70 - 7 * nse, 12 + 2 * nse) |> round() + # el tratamiento lo asignamos al azar, no responde a nse ahora + trata <- rbinom(n, 1, 0.5) + #trata <- rbinom(n, 1, 0.8 * nse + 0.2 * (1 - nse)) + p_trata <- inv_logit(0.5 * trata - 0.2 * (peso - 70)) + res <- rbinom(n, 1, p_trata) + tibble(nse, peso, trata, res) +} +``` + +```{r} +datos_sim_intervencion <- simular_intervencion(5000) +datos_sim_intervencion |> + group_by(trata) |> + summarise(p_res = mean(res)) |> + pivot_wider(names_from = trata, values_from = p_res) |> + mutate(dif_p_trata = `1` - `0`) +``` +**Esta es la cantidad correcta que queremos estimar.** + + +**Pregunta 1**: Explica por qué estos dos resultados son diferentes. Justifica +que si planeamos hacer una intervención general sobre la población, +esta útlima cantidad es nuestro objetivo a estimar. + + + +## Parte 2: Usando datos y el critero de puerta trasera + +Ahora supongamos que tenemos los datos pero no conocemos todos los detalles +del proceso generador. Solamente consideramos nuestro modelo gráfico como +supuestos. + +En este caso, sabemos que hay un camino no causal que pasa por NSE, que +llega a tratamiento y resultado. Ese camino está abierto originalmente. + +**Pregunta 2 **: explica por qué (recuerda que no medimos nse) no podemos usar +la fórmula de ajuste que estratifica por padres del tratamiento. Explica sin +embargo por qué podemos usar el criterio de puerta trasera para identificar +el efecto causal de interés. + +De modo que el camino no causal en donde +NSE es una bifuración puede ser bloqueado si condicionamos a peso, que es una +variable que sí observamos. Por el +**criterio de puerta trasera** que vimos en clase, entonces es posible modelar +el efecto del tratamiento sobre la respuesta una vez que estratificamos por +peso. + +Usaremos un modelo logístico, y hacemos las simulaciones para el contraste +dentro del modelo + +```{r} +mod_trata <- cmdstan_model("backdoor-logistico.stan") +print(mod_trata) +``` +```{r} +datos_lista <- list(N = nrow(datos_sim), + trata = datos_sim$trata, res = datos_sim$res, + peso = datos_sim$peso) +ajuste <- mod_trata$sample(data = datos_lista, refresh = 1000) +sims <- ajuste$draws(c("dif_trata"), format = "df") +``` + + +Obtenemos el resultado correcto. + +```{r} +resumen <- ajuste$summary(c("dif_trata")) +resumen +``` + + +**Pregunta 3**: Nuestro modelo logístico incluye tratamiento y peso. ¿Dónde +ocurre y como se marginaliza el peso en el código de arriba? + + +**Pregunta 4**: Imagínate que sólo nos interesara saber cuál es el efecto +del tratamiento para la población con un peso de 70 kg. Modifica el código de Stan para calcularlo. +Verifica tu resultado usando el proceso generador modificado que consideramos arriba. + + + + +## Parte 3: Estimando y usando el proceso generador completo + +Otra manera de estimar nuestra cantidad de interés es modelando +el sistemo de manera completa. En este caso, no es necesario usar cálculo-do +(criterio de la puerta trasera), sino: + +- Modelamos todos los nodos como indica nuestro diagrama causal y con modelos +estadísticos apropiados. +- Simulamos de la gráfica "mutilada" que no tiene flechas al tratamiento, que +vamos a fijar. Para el resto de los nodos, los simulamos según el modelo ajustado. + +Observa que esto es lo mismo que hicimos en la Parte 1, teórica, donde conocíamos +todas las relaciones funcionales. Este es un camino más complicado pero +también da los resultados correctos. + + + +```{r} +mod_completo <- cmdstan_model("proceso-completo.stan") +print(mod_completo) +``` + + +```{r} +datos_lista <- list(N = nrow(datos_sim), + trata = datos_sim$trata, + nse = datos_sim$nse, + res = datos_sim$res, + peso = datos_sim$peso) +ajuste_comp <- mod_trata$sample(data = datos_lista, refresh = 1000) +sims_comp <- ajuste_comp$draws(c("dif_trata"), format = "df") +``` + + +Obtenemos otra vez el resultado correcto. + +```{r} +resumen_comp <- ajuste_comp$summary(c( "dif_trata")) +resumen_comp +``` + +**Pregunta 5**: explica en general por qué si ya hemos modelado el sistema completo +que propone el modelo gráfico, no es necesario usar ningún "truco" para estimar +el efecto causal de interés. + + + + + + + +