El uso de dagitty nos ayuda a identificar por qué variables tenemos que condicionar/estratificar. Pero veamos algunos ejemplos si wmulados para que nos quede claro
Intro
Eliminar los caminos de asociación no causales.
Counfounder (Fork)
DAG
Show the code
#| echo: truedag1<-dagify(x~z, y~z+x, exposure ="x", outcome ="y", coords =list( x =c( x =1, z =2, y =3), y =c( x =1, z =2, y =1)))#TODO poner las coordenadas de los nodosggdag(dag1)+theme_dag()
El modelo correcto dado este dag es ajustar por Z. Veamos si es cierto, simulando la estructura del dag.
Mediator o chain
DAG
Show the code
#| echo: truedag_med<-dagify(z~x, y~z , coords =list( x =c( x =1, z =2, y =3), y =c( x =1, z =1, y =1)))#TODO poner las coordenadas de los nodosggdag(dag_med)+theme_dag()
usamos dagitty o ggdag_adjustment_set.
Cuando tenemos un “mediator” podemos estar interesado en el efecto total o en el efecto directo.
Ponemos que el efecto total de x sobre y es 6, pero que está mediado sólo por z, y que el efecto directo de x sobre y es 0.
Efecto total es sin ajustar por z y el efecto directo es ajustando por z
Collider
Normalmente “Don’t touch the collider !!”
DAG
Show the code
collider_dag<-dagify(z~x+y , coords =list( x =c( x =1, z =2, y =3), y =c( x =1, z =2, y =1)))#TODO poner las coordenadas de los nodosggdag(collider_dag)+theme_dag()
En el siguiente ejemplo, donde hay un collider, el verdadero efecto de x sobre y es 0 por construcción, pero si ajustamos por el collider se abre el camino no causal y tenemos sesgo.
Parásito de la precisión
Hay veces en que ajustar por pre-treatments es perjudicial. Si no hay variables de confusión “ocultas” se produce una parasitación de la precisión
DAG
Show the code
dag2<-dagify(x~z, y~x, coords =list( x =c( x =1, z =1, y =3), y =c( x =1, z =2, y =1)))#TODO poner las coordenadas de los nodosggdag(dag2)+theme_dag()
A veces, si ajustamos por un pretreatment pero hay una variable de confusión no observada, podemos incurrir en amplificación del sesgo
DAG
Show the code
dag3<-dagify(x~z+u, y~x+u, coords =list( x =c( x =1, z =1, u =2, y =3), y =c( x =1, z =2, u =1.5, y =1)))#TODO poner las coordenadas de los nodosggdag(dag3)+theme_dag()
---title: "Buenos y malos controles"format: live-html: fig-height: 5 fig-dpi: 300 fig-width: 8 fig-align: center code-fold: true code-link: true code-summary: "Show the code" code-tools: true toc: true toc-depth: 2 resources: - data - docs/web_userengine: knitrexecute: warning: false message: falsewebr: packages: - tidyverse repos: - https://r-lib.r-universe.devknitr: opts_chunk: out.width: 80% fig.showtext: TRUE comment: "#>"---{{< include ./_extensions/r-wasm/live/_knitr.qmd >}}```{r setup, include=FALSE}library(ggdag)library(tidyverse)source("./R/setup.R")source("./R/ggdag-mask.R")```El uso de dagitty nos ayuda a identificar por qué variables tenemos que condicionar/estratificar. Pero veamos algunos ejemplos si wmulados para que nos quede claro ## IntroEliminar los caminos de asociación no causales.## Counfounder (Fork)__DAG__```{r}#| echo: truedag1 <-dagify( x ~ z, y ~ z + x,exposure ="x",outcome ="y",coords =list(x =c(x =1,z =2,y =3 ),y =c(x =1, z =2,y =1 )))#TODO poner las coordenadas de los nodosggdag(dag1) +theme_dag()```usamos dagitty o ggdag_adjustment_set. ```{r}#| echo: truedagitty::adjustmentSets(dag1, exposure ="x", outcome ="y")ggdag_adjustment_set(dag1, exposure ="x", outcome ="y")```El modelo correcto dado este dag es ajustar por Z. Veamos si es cierto, simulando la estructura del dag.```{webr}z <- rnorm(300)x <- 2 + 3*z + rnorm(300)y <- 2 + x + z # Modelo incorrectolm(y ~ x)# Modelo correcto debería dar coeficiente de x = 1lm(y ~ x + z )```## Mediator o chain__DAG__```{r}#| echo: truedag_med <-dagify( z ~ x, y ~ z ,coords =list(x =c(x =1,z =2,y =3 ),y =c(x =1, z =1,y =1 )))#TODO poner las coordenadas de los nodosggdag(dag_med) +theme_dag()```usamos dagitty o ggdag_adjustment_set. Cuando tenemos un "mediator" podemos estar interesado en el efecto total o en el efecto directo.__efecto total__```{r}#| echo: truedagitty::adjustmentSets(dag_med, exposure ="x", outcome ="y", effect ="total")ggdag_adjustment_set(dag_med, exposure ="x", outcome ="y", effect ="total") +theme_dag()```__efecto_directo__```{r}#| echo: truedagitty::adjustmentSets(dag_med, exposure ="x", outcome ="y", effect ="direct")ggdag_adjustment_set(dag_med, exposure ="x", outcome ="y", effect ="direct") +theme_dag()```Simulemos Ponemos que el efecto total de `x` sobre `y` es 6, pero que está mediado sólo por `z`, y que el efecto directo de `x` sobre `y` es 0.```{webr}x <- rnorm(800)z <- 3*x + rnorm(800)y <- 2*z + rnorm(800)d <- tibble(x, y, z)# Dentro de valores similares de Z el efecto de X sobre Y es 0, # por eso hacemos varios cortes en z para verlod |> mutate(z_cut = cut_interval(z, 9)) |> ggplot(aes(x = x, y = y)) + geom_point(aes(color = z_cut)) + geom_smooth(method = "lm", se = FALSE, color = "black", linewidth= 2) + geom_smooth(aes(color= z_cut), method = "lm", se = FALSE)```Efecto total es sin ajustar por `z` y el efecto directo es ajustando por `z````{webr}# efecto totallm(y ~ x, data = d)# efecto directolm(y ~ x + z, data = d)```## ColliderNormalmente "Don't touch the collider !! "DAG```{r}#| echo: truecollider_dag <-dagify( z ~ x + y ,coords =list(x =c(x =1,z =2,y =3 ),y =c(x =1, z =2,y =1 )))#TODO poner las coordenadas de los nodosggdag(collider_dag) +theme_dag()``````{r}#| echo: truedagitty::adjustmentSets(collider_dag, exposure ="x", outcome ="y")ggdag_adjustment_set(collider_dag, exposure ="x", outcome ="y")```En el siguiente ejemplo, donde hay un collider, el verdadero efecto de x sobre y es 0 por construcción, pero si ajustamos por el collider se abre el camino no causal y tenemos sesgo.```{webr, colllider_sim}x <- rnorm(800)y <- rnorm(800)z <- as.factor(rbinom(800,1, plogis(2*x - 2*y + rnorm(800))))d <- tibble(x, y, z)d |> ggplot(aes(x = x, y = y)) + geom_point(aes(color = z)) + geom_smooth(method = "lm", se = FALSE, color = "black", linewidth= 2) + geom_smooth(aes(color= z), method = "lm", se = FALSE)# Prueba a ajustar lo siguiente, puedes descomentar o escribir tu código#lm(y ~ x, data = d)#lm(y ~ x + z, data = d)```## Parásito de la precisiónHay veces en que ajustar por pre-treatments es perjudicial. Si no hay variables de confusión "ocultas"se produce una parasitación de la precisiónDAG```{r}#| echo: truedag2 <-dagify( x ~ z, y ~ x,coords =list(x =c(x =1,z =1,y =3 ),y =c(x =1, z =2,y =1 )))#TODO poner las coordenadas de los nodosggdag(dag2) +theme_dag()``````{r}#| echo: truedagitty::adjustmentSets(dag2, exposure ="x", outcome ="y")ggdag_adjustment_set(dag2, exposure ="x", outcome ="y")``````{webr}z <- rnorm(800)x <- 2 + 2*z + rnorm(800)y <- 2 * x + rnorm(800)d <- tibble(x, y, z)# no hay sesgosummary(lm(y ~ x, data = d))summary(lm(y ~ x + z, data = d))```Simulemos para ver mejor la diferencia```{webr}f <- function(bzx = 2, bxy =2, n = 100){ z <- rnorm(n, 4, 2) x <- 2 + bzx*z + rnorm(n) y <- 2 + bxy*x + rnorm(n) d <- tibble(x, y, z) bx_preciso = coef(lm(y ~ x, data = d))[2] bx_parasito = coef(lm(y ~ x + z, data = d))[2] return(c( unlist(bx_preciso), unlist(bx_parasito)))}sim <- replicate(400, f(n = 200))``````{webr}d <- data.frame(bx_preciso = sim[1,], bx_parasito = sim[2,])d <- d |> pivot_longer(cols =everything())d |> ggplot(aes(x = value, fill = name)) + geom_density()```## Amplificación del sesgoA veces, si ajustamos por un pretreatment pero hay una variable de confusiónno observada, podemos incurrir en amplificación del sesgoDAG```{r}#| echo: truedag3 <-dagify( x ~ z + u, y ~ x + u,coords =list(x =c(x =1,z =1,u =2,y =3 ),y =c(x =1, z =2,u =1.5,y =1 )))#TODO poner las coordenadas de los nodosggdag(dag3) +theme_dag()``````{r}dagitty::adjustmentSets(dag3, exposure ="x", outcome ="y")ggdag_adjustment_set(dag3, exposure ="x", outcome ="y")```Tenemos que ajustar por `u`, pero `u ````{webr}z <- rnorm(200)u <- rnorm(200)x <- 1 + 2*z + u + rnorm(200)y <- 1 + x + u + rnorm(200)d <- tibble(x, y, z)# amgos modelos son erróneos puesto que no podemos ajustar por el confounder no observado# pero ajustar por la variable pre-treatment amplifica el sesgosummary(lm(y ~ x, data = d))summary(lm(y ~ x + z, data = d))```Simulemos para ver mejor la diferencia```{webr}f <- function(bzx = 2, bxy =1, n = 100){ z <- rnorm(n) u <- rnorm(n) x <- 1 + bzx*z + u + rnorm(n) y <- 1 + bxy*x + u + rnorm(n) d <- tibble(x, y, z) bx_sesgado = coef(lm(y ~ x, data = d))[2] bx_sesgo_ampli = coef(lm(y ~ x + z, data = d))[2] return(c( unlist(bx_sesgado), unlist(bx_sesgo_ampli)))}sim <- replicate(400, f(n = 200))``````{webr}d <- data.frame(bx_sesgado = sim[1,], bx_sesgo_ampli = sim[2,])d <- d |> pivot_longer(cols =everything())d |> ggplot(aes(x = value, fill = name)) + geom_density()```Y vemos que se produce la amplificación del sesgo