diff --git a/.nojekyll b/.nojekyll
index 71ac78b..7843991 100644
--- a/.nojekyll
+++ b/.nojekyll
@@ -1 +1 @@
-14bc3736
\ No newline at end of file
+c5f0913a
\ No newline at end of file
diff --git a/01-introduccion.html b/01-introduccion.html
index cb3f92f..351ad2c 100644
--- a/01-introduccion.html
+++ b/01-introduccion.html
@@ -230,6 +230,12 @@
5Modelos gráficos y causalidad
+
+
Es decir, el tamaño de los cálculos es una causa común de tratamiento (T) y resultado (M). Veremos más adelante que la decisión de condicionar a el tipo de cálculos proviene de un análisis relativamente simple de este diagrama causal, independientemente de los métodos que usemos para estimar las proporciones de interés (en este ejemplo, examinar las tablas cruzadas es equivalente a hacer estimaciones de máxima verosimlitud).
@@ -744,8 +750,8 @@
Eje
", width =200, height =50)
-
-
+
+
Nótese que el análisis más apropiado no está en los datos: en ambos casos la tabla de datos es exactamente la misma. Los supuestos acerca del proceso que genera los datos sin embargo nos lleva a respuestas opuestas.
Donde vemos ahora que el estado real de cada persona de la prueba es desconocido, aunque el resultado de la prueba depende de ese estado, y la cantidad de positivos que observamos es ahora \(N_{obs}\), que depende también de la sensibilidad y especificidad de la prueba.
@@ -539,8 +545,8 @@
")#, width = 200, height = 50)
-
-
+
+
Usando argumentos como los del modelo original, las distribuciones de esp y sens son beta y podemos incorporarlas en la simulación de la posterior. Nuestra nueva función para simular el proceso generativo es:
Que también podríamos simplificar (suponiendo la \(N\) fija y conocida, pues \(N_+\) y \(M\) dan \(N_{-}\)) como:
@@ -349,8 +355,8 @@
", width =300, height =100)
-
-
+
+
Y ahora construimos el modelo generativo. Supondremos que la muestra de \(N\) personas se toma de manera aleatoria de la población (una población grande, así que podemos ignorar el efecto de muestreo). Supondremos provisionalmente, además, que la prueba es perfecta, es decir, no hay falsos positivos o negativos.
Nótese que no consideramos \(W\to H\), porque podemos pensar en varias intervenciones que podrían cambiar el peso por no cambian la estatura. Por otro lado, es difícil pensar en alguna intervención que cambie la estatura pero no cambie el peso de una persona. Adicionalmente, hay otros factores desconocidos no observados \(U\) que afectan el peso de cada persona adicionalmente a su estatura.
@@ -916,8 +922,8 @@
", width =200, height =50)
-
-
+
+
Omitiendo del diagrama las variables no observadas que también son causas únicamente de \(S\) y \(W, H\):
@@ -942,8 +948,8 @@
", width =200, height =50)
-
-
+
+
Si queremos saber cómo influye el sexo en el peso, este diagrama indica que hay dos tipos de preguntas que podemos hacer:
En este caso, el modelo causal es como sigue: conocemos la distancia \(D\) al hoyo en cada tiro. El éxito (\(Y=1\)) o fracaso (\(Y=0\)) depende de la distancia, junto con la velocidad a la que sale la pelota (muy alto o muy bajo puede dar un tiro fallido), y el ángulo \(\theta\) de salida. Adicionalmente, hay otros factors \(U\) que pueden afectar la probabilidad de éxito. Nótese que no escribiríamos, por ejemplo \(Y \leftarrow D\), porque la distancia no cambia causalmente con el resultado del tiro, aunque es cierto que si intervenimos en la distancia, esperaríamos obtener tasas de éxito diferentes. Igualmente, es necesario poner una flecha de \(V\) a \(D\) y \(V\) a \(Y\).
@@ -2181,7 +2187,7 @@
Warning: 236 of 4000 (6.0%) transitions hit the maximum treedepth limit of 10.
diff --git a/05-dags.html b/05-dags.html
index b2fe615..225bfe0 100644
--- a/05-dags.html
+++ b/05-dags.html
@@ -84,6 +84,7 @@
+
@@ -227,6 +228,12 @@
5Modelos gráficos y causalidad
+
+
Nótese que no describimos exactamente cómo son las funciones que relacionan las variables, sino más bien qué variables son causas directas de qué otras. Por ejemplo, aunque en nuestro ejemplo de arriba \(Y\) puede estar correlacionado con \(Z\), no hay una causa directa a \(Y\), porque cambios en \(Z\) afectan a \(X\), y es el cambio en \(X\) que es causa directa de \(Y\).
@@ -386,8 +393,8 @@
")
-
-
+
+
En este ejemplos no podemos saber \(U1\) y \(U2\), y no nos interesa modelar la física de monedas, manera de lanzarlas, etc. En este ejemplo también no consideraremos qué hace que un día sea soleado o lluvioso (no nos interesa modelar el clima). En este momento, en teoría tenemos ecuaciones determinísticas para todas las variables, y si conocemos todas las variables exógenas \(X,U1,U2,U3\) podríamos determinar exactamente lo que va a suceder con la ganancia, por ejemplo, o cualquier otra variable del sistema.
Sin embargo, si condicionamos a \(Z\), que puede tomar los valores 0 o 1, vemos que \(X\) y \(Y\) son independientes, o dicho de otra manera, la condicional de \(Y\) dada \(Z\) y \(X\) sólo depende de \(Z\):
Un ejemplo con variables continuas podría ser como sigue:
@@ -746,8 +753,8 @@
E
", width =200, height =50)
-
-
+
+
Por la discusión de arriba, es claro que es necesario considerar la edad al casarse si queremos estimar el efecto de tasa de matrimonio en la tasa de divorcio. Es posible que la correlación entre estas dos tasas puede ser explicada solamente por la edad al casarse, y que en realidad al flecha \(M\to D\) sea muy débil o inexistente.
@@ -898,8 +905,8 @@
", width =200, height =50)
-
-
+
+
Es decir, borramos todas las flechas que caen en \(M\) (pues la estamos interveniendo al valor que queramos), y luego simulando \(D\).
Como el NSE es del hogar (una medida general de estatus social), se consideró en principio como una variable pre-tratamiento a la inteligencia de los niños por la que tradicionalmente se controlaba. Burks notó que hacer esto tenía no era apropiado, pues tiene como consecuencia cortar parte del efecto total de la inteligencia sobre el la inteligencia de los hijos. En otras palabras: la inteligencia de los padres hace más probable mejor NSE, y mejor NSE presenta mejores condiciones de desarrollo para sus hijos. Estatificar por esta variable bloquea este efecto.
@@ -1164,8 +1171,8 @@
", width =200, height =50)
-
-
+
+
@@ -1205,9 +1212,9 @@
cor(sims_colisionador |>select(x,y))
-
x y
-x 1.0000000000 0.0004289729
-y 0.0004289729 1.0000000000
+
x y
+x 1.00000e+00 3.52394e-05
+y 3.52394e-05 1.00000e+00
Sin embargo, si condicionamos a \(Z\), que puede tomar los valores 0 o 1:
Consideremos la relación entre Z y Y. Primero vemos que hay dos caminos entre \(Z\) y \(Y\), que son \(p_1:X\gets V \to S\) y \(p_2: Z\to W \gets X \to Y\)
@@ -1672,8 +1679,8 @@
Ejercicio
")
-
-
+
+
@@ -1710,8 +1717,8 @@
E
")
-
-
+
+
@@ -1769,8 +1776,8 @@
")#, width = 200, height = 50)
-
-
+
+
Vimos que para calcular el efecto directo de \(F\) sobre \(W\), por ejemplo, es necesario bloquear el camino que pasa por \(G\) (estratificar por este nodo). Para el efecto total no es necesario condicionar a ningún otro nodo.
En esta sección veremos cómo utilizar la lógica de los diagramas causales que vimos en la sección anterior para entender la posibilidad de identificar efectos causales, es decir, entender si es posible desarrollar estrategias para estimar esos efectos causales. Enfatizamos que este proceso es uno lógico que se deriva de nuestro análisis de las estructuras básicas en DAGs que vimos anteriormente, más que de una colección de “trucos” o “recetas”.
+
+
6.1 Cambiando el proceso generador de datos
+
Comenzamos con el ejemplo más simple de una variable confusora:
+
+
grViz("
+ digraph {
+ node [shape = plaintext];
+ X [label = 'X'];
+ Y [label = 'Y'];
+ U [label = 'U'];
+ X -> Y;
+ U-> X ;
+ U -> Y;
+ {rank = same; X; Y;}
+ }
+ ")
+
+
+
+
+
+
Nos interesa estimar el efecto causal de \(X\) sobre \(Y\). Sucede que en muchas ocasiones existen variables como \(U\) que son causas comunes de \(X\) y \(Y\). Como vimos, esto implica que no podemos simplemente ver la correlación entre \(X\) y \(Y\) para entender el efecto de \(X\) sobre \(Y\), pues una causa común de variación conjunta entre estas dos variables. Esta variable \(U\) puede ser observada o no.
+
Este tipo de confusores ocurren muchas veces en datos observacionales (es decir, de un proceso o sistema que funcione sin intervención de los investigadores). Por ejemplo, si un estudio observa que aquellos que se aplicaron (voluntariamente) un tratamiento \(X\), tienen menor riesgo de hospitalización \(Y\) por cierta enfermadad. Sin embargo, se observa también que aquellos que se aplicaron el tratamiento tienen menos riesgo de tener accidentes viales. Esto indica que la observación de la reducción de riesgo de hospitalización entre los que escogieron el tratamiento probablemente se debe al menos en parte a una variable confusora (por ejemplo, qué tipo de actividades hacen, qué tan cautelosos son, etc.)
+
+
6.1.1 Experimentación
+
Cuando es posible, podemos proponer generar nuevos datos donde alteramos el proceso generador. Una forma muy efectiva y útil, que es muy conveniente cuando es posible, es controlar la asignación del tratamiento. Si en el diagrama anterior, diseñamos un estudio donde observamos a un grupo de personas para las cuales el tratamiento se asignó de acuerdo a un proceso aleatorio, entonces el nuevo diagrama para este nuevo proceso generador es:
+
+
grViz("
+ digraph {
+ node [shape = plaintext];
+ X [label = 'X'];
+ Y [label = 'Y'];
+ R
+ U [label = 'U'];
+ R -> X
+ X -> Y;
+ U -> Y;
+ {rank = same; X; Y;}
+ }
+ ")
+
+
+
+
+
+
Nótese que:
+
+
La variable \(R\) no puede ser endógena (es decir, ninguna flecha del sistema puede incidir en ella), pues se utiliza un dado o algo totalmente no relacionado para asignar el tratamiento. Por ejemplo, también podríamos asignar el tratamiento utilizando la segunda letra del apellido de las personas.
+
No puede existir una flecha de \(U\) a \(X\), pues nada en \(X\) responde a cambios en \(X\), qué solo depende del proceso de aleatorización \(R\).
+
+
En este caso, no es necesario estratificar por ninguna variable y podemos proponer directamente un modelo estadístico para \(Y\) en función de \(X\) que nos permita estimar el efecto causal de \(X\) sobre \(Y\).
+
+
+
+
+
+
+Experimentos
+
+
+
+
Esto describe la idea básica de un experimento simple: es una herramienta para modificar el proceso generador de datos que nos permite identificar efectos causales de manera relativamente simple.
+
Cuando es posible hacer experimentos de calidad, esta puede ser la mejor forma de estimar efectos causales.
+
+
+
En muchos casos, sin embargo, no es posible hacer experimentos de calidad. Hay varias diversas razones, por ejemplo cuando se trata de experimentos que involucran personas:
+
+
No es ético aleatorizar: es totalmente inaceptable asignar aleatoriamente a personas a un tratamientos como fumar 20 cigarros al día, o aleatorizar a niños a recibir educación o no.
+
Aleatorización imposible o imperfecta: no es posible lograr un control total sobre la asignación del tratamiento, y la adherencia al tratamiento asignado de las personas puede variar (por ejemplo, uso de tapabocas en escuelas). A lo más podemos considerar los efectos de una política que intenta tratar a una selección aleatoria de individuos (IIT, o intent-to-treat).
+
+
Así que muchas preguntas causales no están sujetas a modificaciones del proceso generador de datos mediante aleatorización, y es necesario recurrir a otras estrategias.
+
+
+
+
6.2 El operador do
+
Regresamos al diagrama original donde \(U\) es una causa común de \(X\) y \(Y\), y que no tenemos recursos o no es posible hacer un experimento. ¿Existe algún procedimiento estadístico que nos permita estimar el efecto causal de \(X\) sobre \(Y\)?
+
Escribiremos la distribución condicional de la respuesta \(Y\) dada una manipulación de \(X\) como sigue (es decir, en la situación experimental):
+
\[p(Y| do(X=x))\]
+
Esto significa: ¿cómo se distribuye la \(Y\) dado que intervenimos en la población completa (aunque podemos también considerar subpoblaciones más adelante) para poner en \(X=x\)? En primer lugar, notemos que esto no es lo mismo que la distribución condicional usual
+
\[p(Y|X=x),\] que siempre podemos estimar directamente de los datos, y no es la que nos interesa. En el siguiente ejemplo vemos la distinción entre las dos distribuciones:
+
+
Ejemplo (Pearl)
+
Supongamos que tenemos el siguiente modelo del diagrama causal:
+
+
+Código
+
grViz("
+digraph {
+ graph [ranksep = 0.2]
+ node [shape=plaintext]
+
+ edge [minlen = 3]
+ T -> A
+ T -> Z
+
+
+}
+")
+
+
+
+
+
+
+
donde \(T\) es la temperatura, \(A\) son las unidades de agua embotellada vendidas y \(Z\) es la actividad de los mosquitos (medido con muestreo, por ejemplo).
+
No interesa contestar la pregunta: ¿qué tanto influyen las ventas de agua embotellada en la actividad de los mosquitos? Del diagrama, sabemos que no hay ningún camino causal de \(Z\) a \(A\), por lo que nuestra respuesta debería ser igual a 0.
+
Sin embargo, sabemos que estas dos variables están asociadas (por el análisis de DAGs), de manera que describir cómo cambia \(p(Z|A)\) cuando condicionamos a distintos valores de \(A\) no responde nuestra pregunta. La distribución \(p(Z|do(A = a))\) nos dice cómo se distribuye \(Z\) cuando manipulamos \(a\) artificialmente. Por ejemplo, si cerramos todas las tiendas un día haciendo \(do(A=0)\), veríamos que esta variable no tiene efecto sobre la actividad de mosquitos, por ejemplo comparado con \(do(A = 10000)\).
+
Ilustramos la diferencia entre \(p(Y|X)\) y \(p(Y|do(X))\) simulando del ejemplo anterior. Supondremos que sólo consideramos un día del año a lo largo de varios años, para no modelar el comportamiento cíclo de la temperatura:
+
+
simular_t <-function(n =10, dia =150){
+# simular un año, alrededor del día 160 (en junio)
+ t_maxima <-rnorm(n, 28, 2)
+ mosquitos <-rpois(n, 250+10* (t_maxima -28))
+ a_unidades <-rnorm(n, 20000+2000* (t_maxima -28), 2000)
+tibble(t_maxima, a_unidades, mosquitos)
+}
+set.seed(128)
+simular_dias <-simular_t(50)
+
+
Si simulamos, vemos que \(mosquitos\) y \(unidades\) son dependientes, pues tenemos un camino abierto dado por la bifurcación en temperatura:
+
+
ggplot(simular_dias, aes(x = a_unidades, y = mosquitos)) +geom_point() +
+geom_smooth(method ="loess", method.args =list(degree =1)) +
+xlab("Ventas de agua embotellada")
+
+
`geom_smooth()` using formula = 'y ~ x'
+
+
+
+
+
+
+
+
+
+
Sabemos que esta asociación no es causal, pues no hay caminos causales entre estas variables dos variables, pero que hay una dependencia debido a la bifurcación en \(T\). La gráfica muestra que la media condicional \(E[M|A=a]\) depende fuertemente de \(a\), lo que quiere decir que \(p(m|a)\) depende de \(a\) fuertemente.
+
+
+
Una intervención simple
+
En este caso, nos interesaría saber qué sucede si alteramos artificalmente el número de botellas de agua vendidas (puedes imaginar distintas maneras de hacer esto).
+
Cuando hacemos esto, quitamos las aristas que van hacia \(A\), pues \(A\) ya no está determinado por el proceso generador de datos. Tenemos entonces la nueva gráfica:
+
+
+Código
+
grViz("
+digraph {
+ graph [ranksep = 0.2]
+ node [shape=plaintext]
+ A
+ edge [minlen = 3]
+ U_t -> T
+ T -> Z
+ U_m -> Z
+{ rank = same; A; Z }
+}
+")
+
+
+
+
+
+
+
En esta nueva gráfica, \(A\) y \(Z\) son independientes, que es la respuesta correcta. Como cambiamos la gráfica, su proceso generador es diferente al original de los datos observados. Sin embargo, en este ejemplo puedes ver por qué es claro que el cambio que hicimos (manipular \(A\) en lugar de que esté determinado por su proceso generador original) no cambia el modelo de \(Z\), de manera que podemos simular de nuestro nuevo proceso generador donde manipulamos \(A\):
+
+
simular_cirugia <-function(n =10, a_unidades = a_unidades){
+# simular un año, alrededor del día 160 (en junio)
+ t_maxima <-rnorm(n, 28, 2)
+#### cirugía #########
+# ahora a_unidades es fijado por nosotros:
+# a_unidades <- rnorm(n, 20000 + 2000 * (t_maxima - 28), 2000)
+ a_unidades <- a_unidades
+######################
+ mosquitos <-rpois(n, 250+10* (t_maxima -28))
+tibble(t_maxima, a_unidades, mosquitos)
+}
+
+
Y ahora simulamos y graficamos \(p(Z|do(A=a))\) para distintos valores de \(a\):
ggplot(simular_dias_2, aes(x = a_unidades, y = mosquitos)) +
+geom_point() +geom_smooth()
+
+
`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
+
+
+
+
+
+
+
+
+
+
y vemos, como esperaríamos, que no hay relación entre unidades de agua embotellada y mosquitos.
+
+
+
+
6.3 Cálculo-do de Pearl
+
El cálculo do nos da reglas para operar con probabilidades que incluyen nuestro operador do de intervención. En este ejemplo particular, veremos cómo es el argumento:
+
Nótese que al intervenir \(A\) hemos modificado el proceso generador. Si la conjunta original tiene distribución \(p\), escribimos \(p_m\) para la conjunta de la gráfica modificada, de manera que \(p(Z|do(A)) = p_m(Z|A)\): con esto podemos pasar de una pregunta causal (lado izquierdo con operador do) a una estadpística (lado derecho).
+
Aunque intuitivamente vimos cómo simular de esta distribución arriba, especificamos abajo qué reglas son las que nos permiten hacer esto: ¿cómo calculamos \(p_m\)?
+
En primer lugar, consideremos la marginal \(p_m(T)\). Esta marginal es invariante a nuestra cirugía, pues la arista \(T\to A\) que eliminamos \(T\) no afecta el proceso que determina \(T\). De modo que la marginal del proceso modificado es igual a la marginal observada:
+
\[p_m(T) = p(T)\] En segundo lugar, tenemos que
+
\[p_m(Z|T=t,A=a) = p(Z|T=t,A=a),\] Pues el proceso por el cual \(Z\) responde a \(T\) y \(A\) es el mismo, no importa si \(A\) fue modificada artificalmente o no.
+
Juntamos estos argumentos. Primero, por definición,
+
\[p(Z|do(A=a)) = p_m(Z|A=a)\]
+
Por la regla de probabilidad total, podemos condicionar todo a \(T\) y marginalizar. La segunda igualdad la obtenemos por la independencia entre \(T\) y \(Z\) en nuestra gráfica modificada (están \(d\) separadas):
+
\[p_m(z|a) = \int p_m(z|a,t)p_m(t|a)dt = \int p_m(z|a,t)p_m(t)dt\] En segunda igualdad, nótese que cambiamos \(p_m(t|a) = p_m(t)\), lo cual podemos verificar pues en la gráfica modificada \(A\) y \(T\) están \(d\)-separados, lo que implica que son condicionalmente independientes.
+
Finalmente, las últimas dos distribuciones podemos extraerlas de los datos, como explicamos arriba \(p_m(z|t,a) = p(z|t,a)\) y \(p_m(t) = p(t),\) y terminamos con la fórmula:
+
\[p(z|do(a))=p_m(z|a) = \int p(z|a,t)p(t)dt \]
+
Las dos distribuciones de la derecha están en el contexto de \(p\), el proceso generador de datos original. Así que podemos estimarlas de los datos observados.
+
+
Este argumento justifica el proceso que hicimos arriba: simulamos primero \(T\) con su proceso generador, y después simulamos \(Z\) condicional a \(A\) y \(T\)según el proceso generador original, el cual no depende de \(A\) en este ejemplo.
+
+
En el caso de arriba, simulamos de la distribución para entender cómo se distribuía \(Z\) dependiendo de modificaciones a \(A\). Muchas veces nos interesa calcular solamente la esperanza condicional, es decir, cuál es el valor esperado de la variable de interés dado el nivel intervenido, es decir:
+
\(E(Z|do(A=a)) = E_m(Z|A =a),\)
+
que mostramos arriba con la línea ajustada. También quisiéramos calcular contrastes particulares, como qué pasaría si las ventas de agua las aumentamos en 10 mil unidades:
+
\[E(Z|do(A=30000)) - E(Z|do(A=20000)),\] que podemos calcular de manera simple con simulación:
Y vemos que no hay diferencia entre las dos medias.
+
+
Ejemplo
+
Ahora hagamos otro ejemplo donde hay una relación causal que queremos estimar. Imaginemos una ciudad en donde temperaturas altas producen desabasto de agua en algunos hogares, debido a un aumento del riego y uso de agua en general. Nos interesa estimar el efecto del desabasto en las compras de agua embotellada. Nuestro diagrama ahora es:
+
+
+Código
+
grViz("
+digraph {
+ graph [ranksep = 0.2]
+ node [shape=plaintext]
+
+ edge [minlen = 3]
+ U_t -> T
+ T -> A
+ T -> D
+ D -> A
+ U_a -> A
+ U_d -> D
+
+{ rank = same; A; D }
+
+}
+")
+
+
+
+
+
+
+
+
simular_t <-function(n =10, dia =150){
+# simular un año, alrededor del día 160 (en junio)
+ t_maxima <-rnorm(n, 28, 2)
+ u <-rnorm(n, 0, 1)
+ desabasto_agua <-1/(1+exp(-(t_maxima -28) + u))
+ unidades <-rnorm(n, 20000+2000* (t_maxima -28) +8000*desabasto_agua, 2000)
+tibble(t_maxima, unidades, desabasto_agua)
+}
+set.seed(128)
+simular_dias <-simular_t(150)
+
+
+
ggplot(simular_dias, aes(x = desabasto_agua, y = unidades)) +
+geom_point() +geom_smooth()
+
+
`geom_smooth()` using method = 'loess' and formula = 'y ~ x'
+
+
+
+
+
+
+
+
+
+
La correlación parece muy fuerte, sin embargo, sabemos que hay un camino no causal de asociación entre estas dos variables.
+
Igual que en ejemplo anterior, vamos a intervenir teóricamente en el desabasto de agua. Después de la cirugía, nuestro diagrama modificado es:
+
+
+Código
+
grViz("
+digraph {
+ graph [ranksep = 0.2]
+ node [shape=plaintext]
+
+ edge [minlen = 3]
+ U_t -> T
+ T -> A
+ D -> A
+ U_a -> A
+{ rank = same; A; D }
+
+}
+")
+
+
+
+
+
+
+
Ahora queremos calcular \(p(a|do(d)) = p_m(a|d)\) en función de los datos. Siguiendo el mismo argumento que en el ejemplo anterior, sabemos que tenemos que estratificar o condicionar a \(T\) para poder usar nuestro proceso generador de observaciones, y obtenemos:
+
\[p(a|do(d))=p_m(a|d) = \int p(a|d,t)p(t)dt \] Aunque a veces es posible calcular analíticamente el lado derecho analíticamente, podemos simular como hicimos en los ejemplos anteriores:
+
+
simular_cirugia <-function(n =10, da =0){
+# simular un año, alrededor del día 160 (en junio)
+ t_maxima <-rnorm(n, 28, 2)
+### cirugía ####
+#u <- rnorm(n, 0, 1)
+ desabasto_agua <- da
+######
+ unidades <-rnorm(n, 20000+2000* (t_maxima -28) +8000*desabasto_agua, 2000)
+tibble(t_maxima, unidades, desabasto_agua)
+}
+set.seed(128)
+simular_dias_c <-map_df(seq(0, 1, 0.1), \(da) simular_cirugia(1000, da = da))
+
+
+
ggplot(simular_dias_c, aes(x = desabasto_agua, y = unidades)) +
+geom_point() +geom_smooth()
+
+
`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
`geom_smooth()` using method = 'loess' and formula = 'y ~ x'
+
+
+
+
+
+
+
+
+
+
Y este es el efecto causal del desabasto de agua. No tenemos medidas de incertidumbre pues conocemos todos los parámetros de los modelos. La media condicional parece ser lineal, así que podríamos resumir con un modelo lineal:
+
+
# Modelo 1 (con datos de intervención)
+lm(unidades ~ desabasto_agua, simular_dias_c)
Aproximadamente, cada incremento en puntos porcentuales de 10% en desabasto incrementa las ventas en unas 800 unidades. Compara con el análisis donde no estratificamos o controlamos por la temperatura:
En resumen, tenemos la primera regla de Pearl de inferencia causal:
+
+
+
+
+
+
+Fórmula de ajuste (Pearl)
+
+
+
+
Sea \(G\) donde los padres de \(X\) son \(Z_1,Z_2\). El efecto causal total de \(X\) en \(Y\) se puede calcular como
+
\[p(y|do(x)) = \int p(y|x, z_1,z_2) p(z_1,z_2)\, dz_1dz_2\] Es decir, condicionamos al valor de \(x\) y todos los padres de \(X\) para calcular \(p(y|x,z_1,z_2)\), y después marginalizamos sobre los padres.
+
+
+
Esta fórmula se extiende a más de dos padres \(Z_1,Z_2,Z_3,\ldots, Z_k\).
+
+
+
+
+
+
+Tip
+
+
+
+
A este proceso se llama de diferentes maneras en distintos contextos:
+
+
Estamos calculando el efecto causal estratificando por las variables \(z\).
+
Controlamos por las variables \(z\) para calcular el efecto causal.
+
+
+
+
Podemos pensar en esta fórmula de dos maneras: en primer lugar, si estamos modelando toda nuestra gráfica causal, podemos simular de la conjunta de la gráfica mutilada:
+
+
Fijando el nivel del tratamiento \(T\)
+
Simulando \(p(z_1,z_2,\ldots, z_k)\) de nuestro modelo completo (y tomar sólo los valores de las \(z\)’s).
+
Usar \(t\) y las \(z\) simuladas para simular \(y\).
+
Al final, nótese que nos quedan simulaciones de \(p_m(y|t)\) (marginalizamos sobre las \(z\)).
+
+
El otro enfoque busca sólo construir modelos para la parte que nos interesa:
+
+
Construir un modelo separado para \(p(z_1, z_2,\ldots, z_k) = p(z)\) (que puede ser difícil si tenemos muchas variables) a partir los datos. Podemos también simular tomando al azar esta variables de nuestros datos.
+
Construir un modelo \(p(y|t, z)\) para simular la \(y\) a partir de los datos.
+
Marginalizar sobre las \(z\)’s para quedarnos con \(p_m(y|t)\)
+
+
Finalmente, si tenemos un modelo \(p(y| t, z)\) podemos también investigar cómo se comporta \(E[y|t_2,z] - E[y|t_1,z]\) para distintos combinaciones de valores de \(Z\).
+
Nota 1: Con este principio podemos resolver algunos problemas, pero no todos. Veremos que en algunos casos existen padres que no son observados, por ejemplo, no es posible condicionar para usar la fórmula de ajuste y es necesario desarrollar otras estrategias.
+
Nota 2: En regresión lineal, cuando incluímos una variable en el modelo (que consideramos una variable control), estamos estratificando por ella: por ejemplo, en el modelo lineal \(U\sim N(m_u(d,t), \sigma_u)\), donde
+
\[m_u = \beta_0 +\beta_1 d + \beta_2 t\] Estamos calculando un estimador para cada valor de \(T=t\), que es:
+
\[m_u = (\beta_0 + \beta_2 t) + \beta_1 d = \gamma_0 + \gamma_1 d\] Esta es una de las maneras más simples de obtener el efecto de \(d\) estratificando por, o controlando por \(t\), siempre y cuando los modelos lineales sean apropiados.
+
Nótese que en este último caso, tenemos que el efecto de \(d\) no depende de las covariables, de forma que no es necesario hacer el promedio sobre la conjunta, es decir, suponemos que el efecto causal es el mismo independientemente de los valores de las variables de control. Sin embargo, este no siempre es el caso.
+
Nota 3 Sin nuestro modelo \(p(y|t,z)\) es lineal, y nos interesa calcular el efecto causal promedio de la variable \(t\), no es necesario promediar por la conjunta de \(p(z)\). Bajo estas condiciones, el efecto causal promedio está simplemente por el coeficiente de \(t\) en el modelo lineal. Sin embargo, si este no es el caso, entonces para estimar el efecto causal promedio es necesario promediar apropiadamente según la fórmula de ajuste.
+
+
+
6.5 Bloqueando puertas traseras
+
En las partes anteriores vimos que estratificando por los padres de la variable de tratamiento \(X\) podemos construir un estimador del efecto de \(X\) sobre otra variable \(Y\), pasando de una distribución observacional a una conceptualmente experimental (dado que los supuestos causales sean aproximadamente correctos).
+
Sin embargo, esta aplicación de la fórmula de ajuste no funciona si existen padres que no fueron observados, y por tanto no podemos estratificar por ellos. El siguiente método (ajuste por “puerta trasera”) nos da una técnica adicional que podemos usar dado ciertos tipos de estructura en nuestro modelo causal, y presenta una mejoría sobre la fórmula de ajuste simple (veremos también por ejemplo, que a veces podemos usar menos variables que padres de la variable de interés). Nótese que una vez más, este criterio sólo depende de la gráfica causal \(G\) asociada a nuestro modelo, y no los modelos locales que utilizemos para modelar la condicional de cada nodo.
+
+
+
+
+
+
+Ajuste de puerta trasera (Pearl)
+
+
+
+
Si tenemos dos variables \(T\) y \(Y\) en una gráfica \(G\), un conjunto \(Z\) de variables satisface el criterio de puerta trasera relativo a \(T\) y \(Y\) cuando \(Z\) bloquea cualquier camino entre \(T\) y \(Y\) que tenga una arista que incida en \(T\), y ninguna variable de \(Z\) es descendiente de \(T\).
+
En tal caso, podemos utilizar la fórmula de ajuste, pero en lugar de estratificar por los padres de \(T\), estratificamos por las variables en \(Z\)
+
+
+
La idea es:
+
+
Queremos bloquear todos los caminos no causales entre \(T\) y \(Y\).
+
Queremos no perturbar todos los caminos dirigidos de \(T\) a \(Y\) (caminos causales).
+
No queremos activar caminos no causales entre \(T\) y \(Y\) al condicionar.
+
+
Cumplimos 1 al estratificar por variables que bloquean los caminos que son causas de \(T\), pues estos caminos no son causales y distorsionan la relación entre \(T\) y \(Y\). Al mismo tiempo, no bloqueamos caminos causales porque ningúna variable de \(Z\) es descendiente de \(T\), de modo que se satisface el criterio 2 (todos los caminos causales comienzan con \(T\to\)). Finalmente, al excluir descendientes de \(T\) también implica que no condicionamos a colisionadores del tipo \(T\to \cdots \to Z_1\gets Y\), pues esto activa un camino no causal entre \(T\) y \(Y\) (se cumple 3).
+
+
Ejemplo (Pearl)
+
Consideramos primero este ejemplo simple, donde queremos evaluar la efectividad de un tratamiento en cierta enfermedad. Los datos que tenemos disponibles son si una persona recibió 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. Nuestros supuestos están en la siguiente gráfica:
+
+
+Código
+
grViz("
+digraph {
+ graph [ranksep = 0.2, rankdir = LR]
+ node [shape=plaintext]
+ Trata
+ Res
+ node [shape = circle]
+ NSE
+ Peso
+ U
+ edge [minlen = 3]
+ NSE -> Peso
+ NSE -> Trata
+ Trata -> Res
+ Peso -> Res
+ U -> NSE
+ U -> Peso
+}
+")
+
+
+
+
+
+
+
Observamos que no podemos directamente usar la fórmula de ajuste pues NSE no es una variable observada.
+
En esta circunstancia no podríamos identificar el efecto causal, pues existen un caminos abiertos no causales. Quizá el tratamiento no es muy efectivo, y parece ser bueno pues fue aplicado a personas con menor peso que las que no recibieron el tratamiento, a través del efecto de NSE. Sin embargo, supón que tuviéramos disponible la variable Peso:
+
+
+Código
+
grViz("
+digraph {
+ graph [ranksep = 0.2, rankdir = LR]
+ node [shape=plaintext]
+ Trata
+ Res
+ Peso
+ node [shape = circle]
+ NSE
+ U
+ edge [minlen = 3]
+ NSE -> Peso
+ NSE -> Trata
+ Trata -> Res
+ Peso -> Res
+ U -> NSE
+ U -> Peso
+}
+")
+
+
+
+
+
+
+
En este caso, todavía no podemos aplicar la fórmula original de ajuste pues no conocemos \(NSE\). Sin embargo, podemos bloquear los caminos no causales estratificando por Peso, y entonces podemos usar el criterio de puerta trasera para identificar el efecto del tratamiento, aún cuando no tengamos NSE.
# A tibble: 2 × 5
+ trata res n p dif
+ <int> <int> <int> <dbl> <dbl>
+1 0 1 2678 0.533 NA
+2 1 1 3686 0.741 0.208
+
+
+
Sabemos que esta diferencia en respuesta puede estar confundida por un camino no causal. El verdadero efecto casual podemos calcularlo en nuestras simulaciones como sigue a partir de nuestro modelo (igualmente, usamos una muestra muy grande):
# A tibble: 2 × 5
+ trata res n p dif
+ <int> <int> <int> <dbl> <dbl>
+1 0 1 5002 0.500 NA
+2 1 1 7344 0.735 0.235
+
+
+
Suponiendo nuestro diagrama, queremos estimar estratificando por peso. Podríamos usar un sólo modelo logístico, pero pueden ser más simples los cálculos si construimos nuestro modelo en stan. En este caso, podríamos calcular las diferencias para un peso particular, por ejemplo 70 kg (en lugar de modelar estaturas para producir una estimación de diferencia promedio).
Warning: Dropping 'draws_df' class as required metadata was removed.
+
+
+
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
+
+
+
+
+
+
+
+
+
+
Y obtenemos una estimación correcta del efecto en 70 kg. Podríamos también calcular el efecto en distintos pesos (nuestro estimador es una curva), promediar estimando una distribución de pesos modelada, o tomar una distribución fija de pesos para modelar (cada una de estas estrategias tiene propósitos diferentes).
+
Si queremos tener un efecto promedio, podemos modelar los pesos. Otra estrategia es promediar sobre los valores observados de la muestra. Nótese que esto ignora una parte de la incertidumbre proveniente de la muestra particular usada.
Warning: Dropping 'draws_df' class as required metadata was removed.
+
+
+
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
+
+
+
+
+
+
+
+
+
+
Y recuperamos nuevamente el efecto verdadero que mostramos arriba.
+
+
+
+
6.6 Reglas del cálculo-do (opcional)
+
Existen tres axiomas básicos del cálculo-do de las que se derivan los demás resultados, como veremos en el siguiente ejemplo del criterio de la puerta delantera.
+
Antes de verlas, un resumen rápido de las reglas es el siguiente:
+
+
La regla 1 nos dice que las distribuciones asociadas a intervenciones satisfacen también la equivalencia de \(d\)-separación e independencia condicional: si \(Y\) y \(Z\) están \(d\)-separadas dado en la gráfica manipulada, entonces \(p(y | do(x), z) = p(y|do(x))\).
+
La regla 2 es el criterio de la puerta trasera: si condicionamos a variables \(W\) que bloquean toda puerta trasera de \(X\) a \(Y\), podemos cambiar \(do(x)\) por \(x\): \(p(y | do(x), w) = p(y | x, w)\).
+
La regla 3 expresa que si no hay caminos causales de \(X\) a \(Y\), entonces \(p(y|do(x)) = p(y)\).
+
+
+
+
+
+
+
+Completitud (Shpitser, Pearl)
+
+
+
+
Si un efecto causal es identificable (puede expresarse en términos de cantidades observacionales), entonces puede derivarse una estrategia de identificación a partir de las tres reglas del cálculo-do.
+
+
+
Nota: esto no excluye que bajo ciertas hipótesis adicionales a las de nuestra gráfica causal (por ejemplo cómo se comportan las distribuciones particulares qeu componen el modelo), sea posible identificar efectos causales con otros medios que van más allá del cálculo-do.
+
Con más generalidad, abajo están estas reglas (donde condicionamos a más variables o hacemos más intervenciones, y afinamos las condiciones):
+
Denotamos por \(G_m\) la gráfica mutilada por \(do(x)\), donde quitamos todas las aristas que entran en \(X\). Los tres axiomas son:
+
Regla 1 Ignorar observaciones: Si \(Y\) y \(Z\) están \(d\)-separados por \(X\) y \(W\) en \(G_m\),
+
\[ p(y|do(x), z, w) = p(y|do(x), w)\] O en otras palabras, si \(p_m\) es la conjunta para \(G_m\),
+
\[p_m(y|x,z,w) = p_m(y|x, w)\] es cierto si \(Y\) y \(Z\) están \(d\)-separados por \(X\) y \(W\) en \(G_m\) (condicionalmente independientes). Así que esta regla es independencia condicional dado \(d\)-separación, pero para la gráfica intervenida.
+
Regla 2 Usando observaciones como intervenciones:
+
Si \(Y\) y \(Z\) están \(d\)-separados por \(X\) y \(W\) en \(G_m\) quitándole todas las aristas que salen de \(Z\), entonces
Si \(Z\) y \(Y\) están \(d\)-separadas por \(X\) y \(W\) en la gráfica \(G_m\) donde además quitamos cualquier arista a \(Z\) si \(Z\) no es antecesor de \(W\) en \(G_m\), entonces:
+
\[ p(y|do(x), do(z), w) = p(y|do(x), w)\]
+
+
+
6.7 El criterio de puerta delantera
+
En algunos casos, puede ser que no sea posible bloquear algún camino no causal con variables observadas. Un ejemplo clásico es el de la discusión acerca de la relación de fumar con cáncer de pulmón. Algunos estadísticos plantearon que los estudios de asociación entre fumar y cáncer de pulmón podrían tener efectos gravemente confundidos, por ejemplo, por aspectos genéticos que hacen a una persona propensa a fumar al mismo tiempo que aumenta su probabilidad de fumar:
+
+
+Código
+
grViz("
+digraph {
+ graph [ranksep = 0.2]
+ node [shape=plaintext]
+ F
+ C
+ node [shape = circle]
+ U
+ edge [minlen = 3]
+ U -> F
+ U -> C
+ F -> C
+{rank= same; C; F}
+}
+")
+
+
+
+
+
+
+
En este caso, el efecto de fumar (\(F\)) sobre cáncer (\(C\)) no es identificable pues no podemos condicionar a la variable de Genotipo (\(U\)). Supongamos que tenemos una medida adicional, que es la cantidad de depósitos de alquitrán den los pulmones de los pacientes. Este es es afectado por \(F\), y a su vez, el alquitrán incrementa la probabilidad de cáncer:
+
+
+Código
+
grViz("
+digraph {
+ graph [ranksep = 0.2]
+ node [shape=plaintext]
+ F
+ C
+ A
+ node [shape = circle]
+ U
+ edge [minlen = 3]
+ U -> F
+ U -> C
+ F -> A
+ A -> C
+{rank= same; C; F; A}
+}
+")
+
+
+
+
+
+
+
La idea es primero estimar el efecto de \(F\) sobre \(A\), y después estimar el efecto de \(A\) sobre \(C\). La “composición” de estos dos efectos, dado el diagrama, debe darnos el estimador correcto. Primero consideramos el efecto de \(F\) sobre \(A\), y tenemos que (regla 2)
+
\[p(a|do(f)) = p(a|f),\] La igualdad se debe a que una vez que condicionamos a \(F\) no hay puertas traseras entre \(F\) y \(A\) (pues no condicionamos a \(C\)). Esta dependencia causal la podemos entonces estimar de los datos.
+
El efecto de \(A\) sobre \(C\) también es identificable, pues el camino no causal se bloquea cuando condicionamos a \(A\), de forma que por la fórmula de ajuste:
+
\[p(c|do(a)) = \int p(c|a, f') p(f')\, df'\]
+
Ahora encadenamos estas dos ecuaciones:
+
\[p(c|do(f)) = \int p(c|do(a))p(a|f)\,da\]
+
que equivale en simulación a: dado un valor de \(F\), simulamos \(A=a\) con nuestro modelo ajustado con datos naturales. Ahora intervenimos \(A\) con el valor a que obtuvimos y simulamos \(C\). Sin embargo, para hacer este último paso con datos naturales, necesitamos usar el criterio de puerta trasera como explicamos arriba: simulamos entonces \(f´\) de \(p(f)\), y después simulamos \(C\) en función de \(a\) y \(f´\) (con una distribución construida a partir de datos).
+
Requerimos en este caso construir y estimar la condicional \(p(c|a, f)\) basado en los datos.
+
En fórmula, en general, se escribe como:
+
+
+
+
+
+
+Criterio de fuerta delantera (Pearl)
+
+
+
+
Decimos que un conjunto de variables \(A\) satisface el criterio de puerta delantera en relación a las variables \(F\) y \(C\) cuando:
+
+
\(A\) intercepta todos las cadenas dirigidos de \(F\) a \(C\)
+
No hay ningún camino activo de puerta trasera de \(F\) a \(A\)
+
Todos los caminos de puerta trasera de \(A\) a \(C\) están bloqueados por \(F\).
+
+
Si \(A\) satisface el criterio de puerta delantera en relación a \(F\) y \(C\), entonces el efecto causal de \(F\) en \(C\) es identificable y está dado por la fórmula:
Warning: `qplot()` was deprecated in ggplot2 3.4.0.
+
+
+
+
+
+
+
+
+
+
¿Cómo se ve la relación de fumador con cáncer? En esta gráfica mostramos también el valor de la variable no observada \(U\). Nótese que parte de la correlación positiva que existe es debido a esta variable \(U\).
En efecto causal promedio de fumar, en cada nivel, sobre la incidencia de cáncer de pulmón, suponiendo nuestro proceso generador. Nótese que la relación no es tan fuerte como observamos en los datos naturales (en rojo). Esto se debe a que en los datos naturales, las personas existe una causa común entre no fumar y prevenir cáncer de pulmón.
+
+
+
Ejemplo: estimación con puerta delantera
+
Veamos cómo sería la estimación si tuviéramos datos disponible, y si es que podemos recuperar el efecto correcto dados los datos observados y la técnica de puerta delantera.
+
Nótese que sólo necesitamos \(p(c|a, f), p(a|f)\) y \(p(f)\). Estos son modelos estadísticos con el que podemos identificar el efecto que nos interesa. Una vez que los estimemos, podemos usar simulación:
+
+
Fijamos una \(f\).
+
Simulamos una \(a\) del modelo \(p(a|f)\)
+
Para calcular \(\int p(c|a,f')p(f')\), tenemos que simular un valor \(f'\) de la marginal de \(p(f)\), y luego, sustituir junto la \(a\) de 1 para simular una \(c\) de \(p(c|a, f')\).
+
Consideramos solamente \(c\) y \(f\) para resumir el efecto.
Y parece que hemos obtenido una estimación razonable del efecto causal de fumar sobre cáncer. Recordemos también que debemos ser cuidadosos al comparar intervalos que salen del mismo modelo por su nivel de traslape.
+
Por ejemplo, si quisiéramos calcular contrastes con el nivel 2 de tabaquismo:
Nota: nótese como en este ejemplo hemos evitado incluir en nuestro modelo la variable no observada \(U\), gracias al procedimiento de puerta delantera descrito arriba.
+
Es posible sin embargo intentar un modelo completo bayesiano, sin necesidad de recordar la fórmula. El procedimiento, que es más difícil de ajustar: considera una variable latente \(U\) no observada, y es necesario definir cómo puede ser su relación con sus descendientes. Es necesario más cuidado en definir formas funcionales e iniciales apropiadas para que los muestreadores funcionen apropiadamente.
+
+
+
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/06-calculo-do_files/figure-html/unnamed-chunk-10-1.png b/06-calculo-do_files/figure-html/unnamed-chunk-10-1.png
new file mode 100644
index 0000000..108b869
Binary files /dev/null and b/06-calculo-do_files/figure-html/unnamed-chunk-10-1.png differ
diff --git a/06-calculo-do_files/figure-html/unnamed-chunk-14-1.png b/06-calculo-do_files/figure-html/unnamed-chunk-14-1.png
new file mode 100644
index 0000000..5db4c1e
Binary files /dev/null and b/06-calculo-do_files/figure-html/unnamed-chunk-14-1.png differ
diff --git a/06-calculo-do_files/figure-html/unnamed-chunk-17-1.png b/06-calculo-do_files/figure-html/unnamed-chunk-17-1.png
new file mode 100644
index 0000000..75de885
Binary files /dev/null and b/06-calculo-do_files/figure-html/unnamed-chunk-17-1.png differ
diff --git a/06-calculo-do_files/figure-html/unnamed-chunk-18-1.png b/06-calculo-do_files/figure-html/unnamed-chunk-18-1.png
new file mode 100644
index 0000000..9e46610
Binary files /dev/null and b/06-calculo-do_files/figure-html/unnamed-chunk-18-1.png differ
diff --git a/06-calculo-do_files/figure-html/unnamed-chunk-30-1.png b/06-calculo-do_files/figure-html/unnamed-chunk-30-1.png
new file mode 100644
index 0000000..045c8fc
Binary files /dev/null and b/06-calculo-do_files/figure-html/unnamed-chunk-30-1.png differ
diff --git a/06-calculo-do_files/figure-html/unnamed-chunk-33-1.png b/06-calculo-do_files/figure-html/unnamed-chunk-33-1.png
new file mode 100644
index 0000000..2d486bd
Binary files /dev/null and b/06-calculo-do_files/figure-html/unnamed-chunk-33-1.png differ
diff --git a/06-calculo-do_files/figure-html/unnamed-chunk-37-1.png b/06-calculo-do_files/figure-html/unnamed-chunk-37-1.png
new file mode 100644
index 0000000..664bc72
Binary files /dev/null and b/06-calculo-do_files/figure-html/unnamed-chunk-37-1.png differ
diff --git a/06-calculo-do_files/figure-html/unnamed-chunk-38-1.png b/06-calculo-do_files/figure-html/unnamed-chunk-38-1.png
new file mode 100644
index 0000000..8e7c203
Binary files /dev/null and b/06-calculo-do_files/figure-html/unnamed-chunk-38-1.png differ
diff --git a/06-calculo-do_files/figure-html/unnamed-chunk-39-1.png b/06-calculo-do_files/figure-html/unnamed-chunk-39-1.png
new file mode 100644
index 0000000..6373ac1
Binary files /dev/null and b/06-calculo-do_files/figure-html/unnamed-chunk-39-1.png differ
diff --git a/06-calculo-do_files/figure-html/unnamed-chunk-43-1.png b/06-calculo-do_files/figure-html/unnamed-chunk-43-1.png
new file mode 100644
index 0000000..190a7e7
Binary files /dev/null and b/06-calculo-do_files/figure-html/unnamed-chunk-43-1.png differ
diff --git a/06-calculo-do_files/figure-html/unnamed-chunk-44-1.png b/06-calculo-do_files/figure-html/unnamed-chunk-44-1.png
new file mode 100644
index 0000000..568f2e7
Binary files /dev/null and b/06-calculo-do_files/figure-html/unnamed-chunk-44-1.png differ
diff --git a/06-calculo-do_files/figure-html/unnamed-chunk-6-1.png b/06-calculo-do_files/figure-html/unnamed-chunk-6-1.png
new file mode 100644
index 0000000..2846970
Binary files /dev/null and b/06-calculo-do_files/figure-html/unnamed-chunk-6-1.png differ
diff --git a/index.html b/index.html
index a74bb2d..346dfd1 100644
--- a/index.html
+++ b/index.html
@@ -161,6 +161,12 @@
5Modelos gráficos y causalidad
+
+
diff --git a/search.json b/search.json
index 3e5604b..7fd82a6 100644
--- a/search.json
+++ b/search.json
@@ -24,7 +24,7 @@
"href": "01-introduccion.html#diagramas-causales",
"title": "1 Introducción",
"section": "",
- "text": "Causas y mecanismos\n\n\n\nLas razones de cómo hacemos análisis estadístico (que procedimiento o algoritmo seleccionamos, por ejemplo) en un problema dado no están en los datos observados, las causas de los datos.\n\n\n\n\n\nEjemplo (cálculos renales)\nEste es un estudio real acerca de tratamientos para cálculos renales (Julious y Mullee (1994)). Pacientes se asignaron de una forma no controlada a dos tipos de tratamientos para reducir cálculos renales. Para cada paciente, conocemos el tipo de ćalculos que tenía (grandes o chicos) y si el tratamiento tuvo éxito o no.\nLa tabla original tiene 700 renglones (cada renglón es un paciente)\n\ncalculos <- read_csv(\"../datos/kidney_stone_data.csv\")\nnames(calculos) <- c(\"tratamiento\", \"tamaño\", \"éxito\")\ncalculos <- calculos |> \n mutate(tamaño = ifelse(tamaño == \"large\", \"grandes\", \"chicos\")) |> \n mutate(resultado = ifelse(éxito == 1, \"mejora\", \"sin_mejora\")) |> \n select(tratamiento, tamaño, resultado)\nnrow(calculos)\n\n[1] 700\n\n\ny se ve como sigue (muestreamos algunos renglones):\n\ncalculos |> \n sample_n(10) |> kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\ntamaño\nresultado\n\n\n\n\nB\ngrandes\nmejora\n\n\nB\nchicos\nmejora\n\n\nA\ngrandes\nmejora\n\n\nA\ngrandes\nmejora\n\n\nA\nchicos\nsin_mejora\n\n\nA\ngrandes\nmejora\n\n\nB\nchicos\nmejora\n\n\nA\ngrandes\nsin_mejora\n\n\nB\nchicos\nmejora\n\n\nB\nchicos\nsin_mejora\n\n\n\n\n\n\n\nAunque estos datos contienen información de 700 pacientes, los datos pueden resumirse sin pérdida de información contando como sigue:\n\ncalculos_agregada <- calculos |> \n group_by(tratamiento, tamaño, resultado) |> \n count()\ncalculos_agregada |> kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\ntamaño\nresultado\nn\n\n\n\n\nA\nchicos\nmejora\n81\n\n\nA\nchicos\nsin_mejora\n6\n\n\nA\ngrandes\nmejora\n192\n\n\nA\ngrandes\nsin_mejora\n71\n\n\nB\nchicos\nmejora\n234\n\n\nB\nchicos\nsin_mejora\n36\n\n\nB\ngrandes\nmejora\n55\n\n\nB\ngrandes\nsin_mejora\n25\n\n\n\n\n\n\n\nComo en este caso nos interesa principalmente la tasa de éxito de cada tratamiento, podemos mejorar mostrando como sigue:\n\ncalculos_agregada |> pivot_wider(names_from = resultado, values_from = n) |> \n mutate(total = mejora + sin_mejora) |> \n mutate(prop_mejora = round(mejora / total, 2)) |> \n select(tratamiento, tamaño, total, prop_mejora) |> \n arrange(tamaño) |> \n kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\ntamaño\ntotal\nprop_mejora\n\n\n\n\nA\nchicos\n87\n0.93\n\n\nB\nchicos\n270\n0.87\n\n\nA\ngrandes\n263\n0.73\n\n\nB\ngrandes\n80\n0.69\n\n\n\n\n\n\n\nEsta tabla descriptiva es una reescritura de los datos, y no hemos resumido nada todavía. Pero es apropiada para empezar a contestar la pregunta:\n\n¿Qué indican estos datos acerca de qué tratamiento es mejor? ¿Acerca del tamaño de cálculos grandes o chicos?\n\nSupongamos que otro analista decide comparar los pacientes que recibieron cada tratamiento, ignorando la variable de tamaño:\n\ncalculos |> group_by(tratamiento) |> \n summarise(prop_mejora = mean(resultado == \"mejora\") |> round(2)) |> \n kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\nprop_mejora\n\n\n\n\nA\n0.78\n\n\nB\n0.83\n\n\n\n\n\n\n\ny parece ser que el tratamiento \\(B\\) es mejor que el \\(A\\). Esta es una paradoja (un ejemplo de la paradoja de Simpson) . Si un médico no sabe que tipo de cálculos tiene el paciente, ¿entonces debería recetar \\(B\\)? ¿Si sabe debería recetar \\(A\\)? Esta discusión parece no tener mucho sentido.\nPodemos investigar por qué está pasando esto considerando la siguiente tabla, que solo examina cómo se asignó el tratamiento dependiendo del tipo de cálculos de cada paciente:\n\ncalculos |> group_by(tratamiento, tamaño) |> count() |> \n kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\ntamaño\nn\n\n\n\n\nA\nchicos\n87\n\n\nA\ngrandes\n263\n\n\nB\nchicos\n270\n\n\nB\ngrandes\n80\n\n\n\n\n\n\n\nNuestra hipótesis aquí es que la decisión de qué tratamiento usar depende del tamaño de los cálculos. En este caso, hay una decisión pues A es una cirugía y B es un procedimiento menos invasivo, y se prefiere utilizar el tratamiento \\(A\\) para cálculos grandes, y \\(B\\) para cálculos chicos. Esto quiere decir que en la tabla total el tratamiento \\(A\\) está en desventaja porque se usa en casos más difíciles, pero el tratamiento \\(A\\) parece ser en general mejor. La razón es probablemente un proceso de optimización de recursos y riesgo que hacen los doctores.\n\nEn este caso, una mejor respuesta a la pregunta de qué tratamiento es mejor es la que presenta los datos desagregados.\nLa tabla desagregada de asignación del tratamiento nos informa acerca de cómo se está distribuyendo el tratamiento en los pacientes.\n\n\n\n\n\n\n\nNota\n\n\n\nLos resúmenes descriptivos acompañados de hipótesis causales acerca del proceso generador de datos, nos guía hacia descripciones interpretables de los datos.\n\n\nLas explicaciones no son tan simples y, otra vez, interviene el comportamiento de doctores, tratamientos, y distintos tipos de padecimientos.\nPodemos codificar la información causal con un diagrama:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n T \n M \n C\n edge [minlen = 3]\n T -> M\n C -> T\n C -> M\n{ rank = same; M; T }\n}\n\", width = 200, height = 50)\n\n\n\n\n\n\nEs decir, el tamaño de los cálculos es una causa común de tratamiento (T) y resultado (M). Veremos más adelante que la decisión de condicionar a el tipo de cálculos proviene de un análisis relativamente simple de este diagrama causal, independientemente de los métodos que usemos para estimar las proporciones de interés (en este ejemplo, examinar las tablas cruzadas es equivalente a hacer estimaciones de máxima verosimlitud).\n\n\nEjemplo (cálculos renales 2)\nContrastemos el ejemplo anterior usando exactamente la misma tabla de datos, pero con el supuesto de un proceso generador diferente. En este caso, los tratamientos son para mejorar alguna enfermedad del corazón. Sabemos que parte del efecto de este tratamiento ocurre gracias a una baja en presión arterial de los pacientes, así que después de administrar el tratamiento, se toma la presión arterial de los pacientes. Ahora tenemos la tabla agregada y desagregada como sigue:\n\ncorazon <- calculos |> \n select(tratamiento, presión = tamaño, resultado) |> \n mutate(presión = ifelse(presión == \"grandes\", \"alta\", \"baja\"))\ncorazon_agregada <- corazon |> \n group_by(tratamiento, presión, resultado) |> \n count()\ncorazon_agregada |> pivot_wider(names_from = resultado, values_from = n) |> \n mutate(total = mejora + sin_mejora) |> \n mutate(prop_mejora = round(mejora / total, 2)) |> \n select(tratamiento, presión, total, prop_mejora) |> \n arrange(presión) |> \n kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\npresión\ntotal\nprop_mejora\n\n\n\n\nA\nalta\n263\n0.73\n\n\nB\nalta\n80\n0.69\n\n\nA\nbaja\n87\n0.93\n\n\nB\nbaja\n270\n0.87\n\n\n\n\n\n\n\n\ncorazon |> group_by(tratamiento) |> \n summarise(prop_mejora = mean(resultado == \"mejora\") |> round(2)) |> \n kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\nprop_mejora\n\n\n\n\nA\n0.78\n\n\nB\n0.83\n\n\n\n\n\n\n\n¿Cuál creemos que es el mejor tratamiento en este caso? ¿Deberíamos usar la tabla agregada o la desagregada por presión?\n\nEn este caso, la tabla agregada es más apropiada (B es mejor tratamiento).\nLa razón es que presión en este caso es una consecuencia de tomar el tratamiento, y como las tablas muestran, B es más exitoso en bajar la presión de los pacientes.\nSi sólo comparamos dentro de los grupos de presión baja o de presión alta, ignoramos lo más importante del tratamiento en la probabilidad de mejorar.\n\nNuestros supuestos causales podemos mostrarlos con el siguiente diagrama:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n P\n T \n M \n edge [minlen = 3]\n T -> P\n P -> M\n T -> M\n{ rank = same; M; T}\n}\n\", width = 200, height = 50)\n\n\n\n\n\n\nNótese que el análisis más apropiado no está en los datos: en ambos casos la tabla de datos es exactamente la misma. Los supuestos acerca del proceso que genera los datos sin embargo nos lleva a respuestas opuestas.",
+ "text": "Causas y mecanismos\n\n\n\nLas razones de cómo hacemos análisis estadístico (que procedimiento o algoritmo seleccionamos, por ejemplo) en un problema dado no están en los datos observados, las causas de los datos.\n\n\n\n\n\nEjemplo (cálculos renales)\nEste es un estudio real acerca de tratamientos para cálculos renales (Julious y Mullee (1994)). Pacientes se asignaron de una forma no controlada a dos tipos de tratamientos para reducir cálculos renales. Para cada paciente, conocemos el tipo de ćalculos que tenía (grandes o chicos) y si el tratamiento tuvo éxito o no.\nLa tabla original tiene 700 renglones (cada renglón es un paciente)\n\ncalculos <- read_csv(\"../datos/kidney_stone_data.csv\")\nnames(calculos) <- c(\"tratamiento\", \"tamaño\", \"éxito\")\ncalculos <- calculos |> \n mutate(tamaño = ifelse(tamaño == \"large\", \"grandes\", \"chicos\")) |> \n mutate(resultado = ifelse(éxito == 1, \"mejora\", \"sin_mejora\")) |> \n select(tratamiento, tamaño, resultado)\nnrow(calculos)\n\n[1] 700\n\n\ny se ve como sigue (muestreamos algunos renglones):\n\ncalculos |> \n sample_n(10) |> kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\ntamaño\nresultado\n\n\n\n\nA\ngrandes\nsin_mejora\n\n\nB\ngrandes\nmejora\n\n\nB\nchicos\nmejora\n\n\nA\ngrandes\nsin_mejora\n\n\nB\nchicos\nmejora\n\n\nB\ngrandes\nmejora\n\n\nA\ngrandes\nmejora\n\n\nB\ngrandes\nsin_mejora\n\n\nB\ngrandes\nmejora\n\n\nB\nchicos\nmejora\n\n\n\n\n\n\n\nAunque estos datos contienen información de 700 pacientes, los datos pueden resumirse sin pérdida de información contando como sigue:\n\ncalculos_agregada <- calculos |> \n group_by(tratamiento, tamaño, resultado) |> \n count()\ncalculos_agregada |> kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\ntamaño\nresultado\nn\n\n\n\n\nA\nchicos\nmejora\n81\n\n\nA\nchicos\nsin_mejora\n6\n\n\nA\ngrandes\nmejora\n192\n\n\nA\ngrandes\nsin_mejora\n71\n\n\nB\nchicos\nmejora\n234\n\n\nB\nchicos\nsin_mejora\n36\n\n\nB\ngrandes\nmejora\n55\n\n\nB\ngrandes\nsin_mejora\n25\n\n\n\n\n\n\n\nComo en este caso nos interesa principalmente la tasa de éxito de cada tratamiento, podemos mejorar mostrando como sigue:\n\ncalculos_agregada |> pivot_wider(names_from = resultado, values_from = n) |> \n mutate(total = mejora + sin_mejora) |> \n mutate(prop_mejora = round(mejora / total, 2)) |> \n select(tratamiento, tamaño, total, prop_mejora) |> \n arrange(tamaño) |> \n kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\ntamaño\ntotal\nprop_mejora\n\n\n\n\nA\nchicos\n87\n0.93\n\n\nB\nchicos\n270\n0.87\n\n\nA\ngrandes\n263\n0.73\n\n\nB\ngrandes\n80\n0.69\n\n\n\n\n\n\n\nEsta tabla descriptiva es una reescritura de los datos, y no hemos resumido nada todavía. Pero es apropiada para empezar a contestar la pregunta:\n\n¿Qué indican estos datos acerca de qué tratamiento es mejor? ¿Acerca del tamaño de cálculos grandes o chicos?\n\nSupongamos que otro analista decide comparar los pacientes que recibieron cada tratamiento, ignorando la variable de tamaño:\n\ncalculos |> group_by(tratamiento) |> \n summarise(prop_mejora = mean(resultado == \"mejora\") |> round(2)) |> \n kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\nprop_mejora\n\n\n\n\nA\n0.78\n\n\nB\n0.83\n\n\n\n\n\n\n\ny parece ser que el tratamiento \\(B\\) es mejor que el \\(A\\). Esta es una paradoja (un ejemplo de la paradoja de Simpson) . Si un médico no sabe que tipo de cálculos tiene el paciente, ¿entonces debería recetar \\(B\\)? ¿Si sabe debería recetar \\(A\\)? Esta discusión parece no tener mucho sentido.\nPodemos investigar por qué está pasando esto considerando la siguiente tabla, que solo examina cómo se asignó el tratamiento dependiendo del tipo de cálculos de cada paciente:\n\ncalculos |> group_by(tratamiento, tamaño) |> count() |> \n kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\ntamaño\nn\n\n\n\n\nA\nchicos\n87\n\n\nA\ngrandes\n263\n\n\nB\nchicos\n270\n\n\nB\ngrandes\n80\n\n\n\n\n\n\n\nNuestra hipótesis aquí es que la decisión de qué tratamiento usar depende del tamaño de los cálculos. En este caso, hay una decisión pues A es una cirugía y B es un procedimiento menos invasivo, y se prefiere utilizar el tratamiento \\(A\\) para cálculos grandes, y \\(B\\) para cálculos chicos. Esto quiere decir que en la tabla total el tratamiento \\(A\\) está en desventaja porque se usa en casos más difíciles, pero el tratamiento \\(A\\) parece ser en general mejor. La razón es probablemente un proceso de optimización de recursos y riesgo que hacen los doctores.\n\nEn este caso, una mejor respuesta a la pregunta de qué tratamiento es mejor es la que presenta los datos desagregados.\nLa tabla desagregada de asignación del tratamiento nos informa acerca de cómo se está distribuyendo el tratamiento en los pacientes.\n\n\n\n\n\n\n\nNota\n\n\n\nLos resúmenes descriptivos acompañados de hipótesis causales acerca del proceso generador de datos, nos guía hacia descripciones interpretables de los datos.\n\n\nLas explicaciones no son tan simples y, otra vez, interviene el comportamiento de doctores, tratamientos, y distintos tipos de padecimientos.\nPodemos codificar la información causal con un diagrama:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n T \n M \n C\n edge [minlen = 3]\n T -> M\n C -> T\n C -> M\n{ rank = same; M; T }\n}\n\", width = 200, height = 50)\n\n\n\n\n\n\nEs decir, el tamaño de los cálculos es una causa común de tratamiento (T) y resultado (M). Veremos más adelante que la decisión de condicionar a el tipo de cálculos proviene de un análisis relativamente simple de este diagrama causal, independientemente de los métodos que usemos para estimar las proporciones de interés (en este ejemplo, examinar las tablas cruzadas es equivalente a hacer estimaciones de máxima verosimlitud).\n\n\nEjemplo (cálculos renales 2)\nContrastemos el ejemplo anterior usando exactamente la misma tabla de datos, pero con el supuesto de un proceso generador diferente. En este caso, los tratamientos son para mejorar alguna enfermedad del corazón. Sabemos que parte del efecto de este tratamiento ocurre gracias a una baja en presión arterial de los pacientes, así que después de administrar el tratamiento, se toma la presión arterial de los pacientes. Ahora tenemos la tabla agregada y desagregada como sigue:\n\ncorazon <- calculos |> \n select(tratamiento, presión = tamaño, resultado) |> \n mutate(presión = ifelse(presión == \"grandes\", \"alta\", \"baja\"))\ncorazon_agregada <- corazon |> \n group_by(tratamiento, presión, resultado) |> \n count()\ncorazon_agregada |> pivot_wider(names_from = resultado, values_from = n) |> \n mutate(total = mejora + sin_mejora) |> \n mutate(prop_mejora = round(mejora / total, 2)) |> \n select(tratamiento, presión, total, prop_mejora) |> \n arrange(presión) |> \n kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\npresión\ntotal\nprop_mejora\n\n\n\n\nA\nalta\n263\n0.73\n\n\nB\nalta\n80\n0.69\n\n\nA\nbaja\n87\n0.93\n\n\nB\nbaja\n270\n0.87\n\n\n\n\n\n\n\n\ncorazon |> group_by(tratamiento) |> \n summarise(prop_mejora = mean(resultado == \"mejora\") |> round(2)) |> \n kable() |> \n kable_paper(full_width = FALSE)\n\n\n\n\ntratamiento\nprop_mejora\n\n\n\n\nA\n0.78\n\n\nB\n0.83\n\n\n\n\n\n\n\n¿Cuál creemos que es el mejor tratamiento en este caso? ¿Deberíamos usar la tabla agregada o la desagregada por presión?\n\nEn este caso, la tabla agregada es más apropiada (B es mejor tratamiento).\nLa razón es que presión en este caso es una consecuencia de tomar el tratamiento, y como las tablas muestran, B es más exitoso en bajar la presión de los pacientes.\nSi sólo comparamos dentro de los grupos de presión baja o de presión alta, ignoramos lo más importante del tratamiento en la probabilidad de mejorar.\n\nNuestros supuestos causales podemos mostrarlos con el siguiente diagrama:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n P\n T \n M \n edge [minlen = 3]\n T -> P\n P -> M\n T -> M\n{ rank = same; M; T}\n}\n\", width = 200, height = 50)\n\n\n\n\n\n\nNótese que el análisis más apropiado no está en los datos: en ambos casos la tabla de datos es exactamente la misma. Los supuestos acerca del proceso que genera los datos sin embargo nos lleva a respuestas opuestas.",
"crumbs": [
"1Introducción"
]
@@ -224,7 +224,7 @@
"href": "03-modelos-genericos.html#ampliando-el-modelo",
"title": "4 Componentes de modelación 1",
"section": "4.4 Ampliando el modelo",
- "text": "4.4 Ampliando el modelo\nEntre los adultos humanos, hombres y mujeres tienen distintas distribuciones de peso y estatura. La variable \\(S\\) (sexo) influye tanto en estatura como en peso. La relación la consideramos causalmente partiendo en \\(S\\):\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.3, rankdir = LR]\n node [shape=circle]\n U\n V\n Z\n node [shape=plaintext]\n H\n W\n S\n edge [minlen = 3]\n H -> W\n U -> W\n S -> H\n S -> W\n V -> H\n Z -> S\n}\n\", width = 200, height = 50)\n\n\n\n\n\n\nOmitiendo del diagrama las variables no observadas que también son causas únicamente de \\(S\\) y \\(W, H\\):\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.3, rankdir = LR]\n node [shape=circle]\n \n node [shape=plaintext]\n H\n W\n S\n edge [minlen = 3]\n H -> W\n S -> H\n S -> W\n\n}\n\", width = 200, height = 50)\n\n\n\n\n\n\nSi queremos saber cómo influye el sexo en el peso, este diagrama indica que hay dos tipos de preguntas que podemos hacer:\n\n¿Cuál es el efecto causal de \\(S\\) sobre \\(W\\) (efecto total) ?\n¿Cuál es el efecto causal directo de \\(S\\) sobre \\(W\\)? Es decir, que no actúa a través de \\(H\\).\n\nAunque tenemos un solo modelo causal, pueden construirse distintos modelos estadísticos para contestar cada pregunta. El modelo causal nos dice que si no tenemos causas comunes de \\(S\\) y \\(H\\) y \\(W\\), entonces podemos estimar el efecto total de \\(S\\) sobre \\(W\\) (esto lo formalizaremos más adelante).\nEmpezamos con el efecto total. Para esto, podemos usar el modelo lineal e ignorar la estatura, donde \\(S_i=2\\) si el individuo \\(i\\) es hombre y \\(S_i=1\\) si el individuo \\(i\\) es mujer.\n\\[\n\\begin{align}\nW_i &\\sim N(\\alpha_{S_i}, \\sigma)\\\\\n\\alpha_1,\\alpha_2 &\\sim N(60, 10) \\\\\n\\sigma &\\sim N^+(0, 20) \\\\\n\\end{align}\n\\] Nótese que tenemos dos posibles medias para el peso, una para hombres y otra para mujeres. La estatura no nos importa porque la pregunta es acerca del efecto total de sexo sobre estatura. Para las iniciales podemos seguir un argumento similar al de arriba.\nNota: esta parametrización es más conveniente que utilizar un indicador (o dummy) de sexo en términos de interpetación y en términos de poner iniciales acordes con el conocimiento del área, aunque estadísticamente son equivalentes.\nEl modelo generador simplificado para este caso puede ser:\n\nsim_peso_mod_s <- function(S, alpha, sigma){\n n <- length(S)\n W <- rnorm(n, alpha[S], sigma)\n tibble(alpha_1 = alpha[1], alpha_2 = alpha[2], \n sigma, S = S, W = W)\n}\n\nDado este modelo generador, ¿cuál es el efecto causal de sexo? Tenemos que definir esta cantidad en términos del modelo. En nuestro caso, definiremos el efecto causal promedio sobre la población, que definimos como la diferencia promedio de estaturas de dos poblaciones: una compuesta enteramente por hombres y otra por mujeres.\n\nset.seed(2021)\n# Fjamos mismos valores de los parámetros para simular dos\n# poblaciones\nsim_hombres <- sim_peso_mod_s(rep(2, 1000), c(55, 70), 10)\nsim_mujeres <- sim_peso_mod_s(rep(1, 1000), c(55, 70), 10)\nmean(sim_hombres$W - sim_mujeres$W)\n\n[1] 14.75203\n\n\n\nVerificación a priori\nAhora generamos una población con estos parámetros y vemos si podemos recuperar el efecto causal promedio sobre la población. Nuestro modelo es como definimos arriba:\n\nlibrary(cmdstanr)\nmod_peso <- cmdstan_model(\"./src/peso-estatura-2.stan\")\nprint(mod_peso)\n\ndata {\n int<lower=0> N;\n vector[N] w;\n array[N] int s;\n}\n\nparameters {\n array[2] real alpha;\n real <lower=0> sigma;\n}\n\ntransformed parameters {\n\n}\n\nmodel {\n // modelo para peso\n w ~ normal(alpha[s], sigma);\n // también se puede escribir como\n // for (i in 1:N) {\n // w[i] ~ normal(alpha[s[i]], sigma);\n // }\n // iniciales\n alpha ~ normal(60, 10);\n sigma ~ normal(0, 20);\n}\n\n\nSimulamos datos y ajustamos el modelo, usando los mismos parámetros fijos:\n\nS_sim <- sample(c(1,2), 1000, replace = TRUE)\ndatos_sim_tbl <- sim_peso_mod_s(S_sim, c(55, 70), 10)\n\n\nmod_2_fit <- mod_peso$sample(\n data = list(N = nrow(datos_sim_tbl), \n s = datos_sim_tbl$S, \n w = datos_sim_tbl$W),\n refresh = 0, seed = 221\n)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.1 seconds.\nChain 2 finished in 0.1 seconds.\nChain 3 finished in 0.1 seconds.\nChain 4 finished in 0.1 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.1 seconds.\nTotal execution time: 0.6 seconds.\n\n\n\nmod_2_fit$summary(c(\"alpha\", \"sigma\"))\n\n# A tibble: 3 × 10\n variable mean median sd mad q5 q95 rhat ess_bulk ess_tail\n <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>\n1 alpha[1] 55.0 55.0 0.435 0.429 54.2 55.7 1.00 4611. 2866.\n2 alpha[2] 70.9 70.9 0.435 0.444 70.2 71.6 1.00 4350. 2891.\n3 sigma 9.77 9.76 0.223 0.227 9.41 10.1 1.00 4210. 3013.\n\n\nNótese que la diferencia de medias poblacionales es de alrededor de 15 cm, que es lo que esperábamos según el cálculo de arriba. Podemos replicar el cálculo que hicimos arriba directamente usando simulación:\n\nPara cada simulación de la posterior calculamos una población hipotética de hombres y otras de mujeres (mismos parámetros)\nCalculamos la diferencia de medias poblacionales\nResumimos con la posterior.\n\nEsto es fácil hacerlo directamente en Stan, pero en este ejemplo lo calcularemos manualmente:\n\nsims_post_tbl <- mod_2_fit$draws() |> as_draws_df() |> \n as_tibble()\nsimular_diferencia_post <- function(sims_post_tbl){\n # Simulamos parámetros de la posterior\n pars <- sample_n(sims_post_tbl, 1) |> \n select(starts_with(\"alpha\"), sigma)\n # Simulamos datos\n sims_hombres <- sim_peso_mod_s(rep(2, 1000), \n alpha = c(pars$`alpha[1]`, pars$`alpha[2]`), pars$sigma)\n sims_mujeres <- sim_peso_mod_s(rep(1, 1000), \n c(pars$`alpha[1]`, pars$`alpha[2]`), pars$sigma)\n diferencia <- mean(sims_hombres$W - sims_mujeres$W)\n # Calculamos la diferencia de medias\n tibble(diferencia = diferencia) |> bind_cols(pars)\n}\n\n\nsimular_diferencia_post(sims_post_tbl)\n\n# A tibble: 1 × 4\n diferencia `alpha[1]` `alpha[2]` sigma\n <dbl> <dbl> <dbl> <dbl>\n1 15.3 55.5 70.4 9.60\n\n\nY ahora calculamos el resumen de interés, que es la posterior del contraste o diferencia entre las dos poblaciones simuladas. Comparamos con la línea en rojo que es la cantidad que establecimos a estimar:\n\nmap_df(1:4000, ~ simular_diferencia_post(sims_post_tbl) |> \n mutate(rep = .x)) |> \nggplot(aes(x = diferencia)) +\n geom_histogram(bins = 50) +\n labs(x = \"Efecto de sexo en estatura hombres vs mujeres (cm)\") +\n geom_vline(xintercept = mean(sim_hombres$W - sim_mujeres$W), \n color = \"red\", linewidth = 1.5)\n\n\n\n\n\n\n\n\nPuedes repetir este ejercicio para distintos valores de los parámetros, como hicimos en los ejemplos de arriba.",
+ "text": "4.4 Ampliando el modelo\nEntre los adultos humanos, hombres y mujeres tienen distintas distribuciones de peso y estatura. La variable \\(S\\) (sexo) influye tanto en estatura como en peso. La relación la consideramos causalmente partiendo en \\(S\\):\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.3, rankdir = LR]\n node [shape=circle]\n U\n V\n Z\n node [shape=plaintext]\n H\n W\n S\n edge [minlen = 3]\n H -> W\n U -> W\n S -> H\n S -> W\n V -> H\n Z -> S\n}\n\", width = 200, height = 50)\n\n\n\n\n\n\nOmitiendo del diagrama las variables no observadas que también son causas únicamente de \\(S\\) y \\(W, H\\):\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.3, rankdir = LR]\n node [shape=circle]\n \n node [shape=plaintext]\n H\n W\n S\n edge [minlen = 3]\n H -> W\n S -> H\n S -> W\n\n}\n\", width = 200, height = 50)\n\n\n\n\n\n\nSi queremos saber cómo influye el sexo en el peso, este diagrama indica que hay dos tipos de preguntas que podemos hacer:\n\n¿Cuál es el efecto causal de \\(S\\) sobre \\(W\\) (efecto total) ?\n¿Cuál es el efecto causal directo de \\(S\\) sobre \\(W\\)? Es decir, que no actúa a través de \\(H\\).\n\nAunque tenemos un solo modelo causal, pueden construirse distintos modelos estadísticos para contestar cada pregunta. El modelo causal nos dice que si no tenemos causas comunes de \\(S\\) y \\(H\\) y \\(W\\), entonces podemos estimar el efecto total de \\(S\\) sobre \\(W\\) (esto lo formalizaremos más adelante).\nEmpezamos con el efecto total. Para esto, podemos usar el modelo lineal e ignorar la estatura, donde \\(S_i=2\\) si el individuo \\(i\\) es hombre y \\(S_i=1\\) si el individuo \\(i\\) es mujer.\n\\[\n\\begin{align}\nW_i &\\sim N(\\alpha_{S_i}, \\sigma)\\\\\n\\alpha_1,\\alpha_2 &\\sim N(60, 10) \\\\\n\\sigma &\\sim N^+(0, 20) \\\\\n\\end{align}\n\\] Nótese que tenemos dos posibles medias para el peso, una para hombres y otra para mujeres. La estatura no nos importa porque la pregunta es acerca del efecto total de sexo sobre estatura. Para las iniciales podemos seguir un argumento similar al de arriba.\nNota: esta parametrización es más conveniente que utilizar un indicador (o dummy) de sexo en términos de interpetación y en términos de poner iniciales acordes con el conocimiento del área, aunque estadísticamente son equivalentes.\nEl modelo generador simplificado para este caso puede ser:\n\nsim_peso_mod_s <- function(S, alpha, sigma){\n n <- length(S)\n W <- rnorm(n, alpha[S], sigma)\n tibble(alpha_1 = alpha[1], alpha_2 = alpha[2], \n sigma, S = S, W = W)\n}\n\nDado este modelo generador, ¿cuál es el efecto causal de sexo? Tenemos que definir esta cantidad en términos del modelo. En nuestro caso, definiremos el efecto causal promedio sobre la población, que definimos como la diferencia promedio de estaturas de dos poblaciones: una compuesta enteramente por hombres y otra por mujeres.\n\nset.seed(2021)\n# Fjamos mismos valores de los parámetros para simular dos\n# poblaciones\nsim_hombres <- sim_peso_mod_s(rep(2, 1000), c(55, 70), 10)\nsim_mujeres <- sim_peso_mod_s(rep(1, 1000), c(55, 70), 10)\nmean(sim_hombres$W - sim_mujeres$W)\n\n[1] 14.75203\n\n\n\nVerificación a priori\nAhora generamos una población con estos parámetros y vemos si podemos recuperar el efecto causal promedio sobre la población. Nuestro modelo es como definimos arriba:\n\nlibrary(cmdstanr)\nmod_peso <- cmdstan_model(\"./src/peso-estatura-2.stan\")\nprint(mod_peso)\n\ndata {\n int<lower=0> N;\n vector[N] w;\n array[N] int s;\n}\n\nparameters {\n array[2] real alpha;\n real <lower=0> sigma;\n}\n\ntransformed parameters {\n\n}\n\nmodel {\n // modelo para peso\n w ~ normal(alpha[s], sigma);\n // también se puede escribir como\n // for (i in 1:N) {\n // w[i] ~ normal(alpha[s[i]], sigma);\n // }\n // iniciales\n alpha ~ normal(60, 10);\n sigma ~ normal(0, 20);\n}\n\n\nSimulamos datos y ajustamos el modelo, usando los mismos parámetros fijos:\n\nS_sim <- sample(c(1,2), 1000, replace = TRUE)\ndatos_sim_tbl <- sim_peso_mod_s(S_sim, c(55, 70), 10)\n\n\nmod_2_fit <- mod_peso$sample(\n data = list(N = nrow(datos_sim_tbl), \n s = datos_sim_tbl$S, \n w = datos_sim_tbl$W),\n refresh = 0, seed = 221\n)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.1 seconds.\nChain 2 finished in 0.1 seconds.\nChain 3 finished in 0.1 seconds.\nChain 4 finished in 0.1 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.1 seconds.\nTotal execution time: 0.7 seconds.\n\n\n\nmod_2_fit$summary(c(\"alpha\", \"sigma\"))\n\n# A tibble: 3 × 10\n variable mean median sd mad q5 q95 rhat ess_bulk ess_tail\n <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>\n1 alpha[1] 55.0 55.0 0.435 0.429 54.2 55.7 1.00 4611. 2866.\n2 alpha[2] 70.9 70.9 0.435 0.444 70.2 71.6 1.00 4350. 2891.\n3 sigma 9.77 9.76 0.223 0.227 9.41 10.1 1.00 4210. 3013.\n\n\nNótese que la diferencia de medias poblacionales es de alrededor de 15 cm, que es lo que esperábamos según el cálculo de arriba. Podemos replicar el cálculo que hicimos arriba directamente usando simulación:\n\nPara cada simulación de la posterior calculamos una población hipotética de hombres y otras de mujeres (mismos parámetros)\nCalculamos la diferencia de medias poblacionales\nResumimos con la posterior.\n\nEsto es fácil hacerlo directamente en Stan, pero en este ejemplo lo calcularemos manualmente:\n\nsims_post_tbl <- mod_2_fit$draws() |> as_draws_df() |> \n as_tibble()\nsimular_diferencia_post <- function(sims_post_tbl){\n # Simulamos parámetros de la posterior\n pars <- sample_n(sims_post_tbl, 1) |> \n select(starts_with(\"alpha\"), sigma)\n # Simulamos datos\n sims_hombres <- sim_peso_mod_s(rep(2, 1000), \n alpha = c(pars$`alpha[1]`, pars$`alpha[2]`), pars$sigma)\n sims_mujeres <- sim_peso_mod_s(rep(1, 1000), \n c(pars$`alpha[1]`, pars$`alpha[2]`), pars$sigma)\n diferencia <- mean(sims_hombres$W - sims_mujeres$W)\n # Calculamos la diferencia de medias\n tibble(diferencia = diferencia) |> bind_cols(pars)\n}\n\n\nsimular_diferencia_post(sims_post_tbl)\n\n# A tibble: 1 × 4\n diferencia `alpha[1]` `alpha[2]` sigma\n <dbl> <dbl> <dbl> <dbl>\n1 15.3 55.5 70.4 9.60\n\n\nY ahora calculamos el resumen de interés, que es la posterior del contraste o diferencia entre las dos poblaciones simuladas. Comparamos con la línea en rojo que es la cantidad que establecimos a estimar:\n\nmap_df(1:4000, ~ simular_diferencia_post(sims_post_tbl) |> \n mutate(rep = .x)) |> \nggplot(aes(x = diferencia)) +\n geom_histogram(bins = 50) +\n labs(x = \"Efecto de sexo en estatura hombres vs mujeres (cm)\") +\n geom_vline(xintercept = mean(sim_hombres$W - sim_mujeres$W), \n color = \"red\", linewidth = 1.5)\n\n\n\n\n\n\n\n\nPuedes repetir este ejercicio para distintos valores de los parámetros, como hicimos en los ejemplos de arriba.",
"crumbs": [
"4Componentes de modelación 1"
]
@@ -274,7 +274,7 @@
"href": "03-modelos-genericos.html#modelos-genéricos-para-ajustar-curvas",
"title": "4 Componentes de modelación 1",
"section": "4.9 Modelos genéricos para ajustar curvas",
- "text": "4.9 Modelos genéricos para ajustar curvas\nOtra posibilidad es utilizar un modelo más flexible creando variables derivadas de la distancia. En este caso, quizá podemos ajustar una curva que sea aceptable desde el punto de vista predictivo, pero no podremos aprender mucho acerca de cómo funciona la probabilidad de éxitos de los tiros de putts\n\n\n\n\n\n\nSplines y ajuste de curvas\n\n\n\nLos splines nos dan una manera estándar de ajustar curvas más flexibles, de tipo polinomial por tramos. Usualmente son numéricamente más conveniente que polinomios.\n\n\nAunque hay muchos tipos de splines (los más comunes son B-splines), para este problema consideraremos una base de splines cuadráticos que resultan en curvas monótonas (I-splines). Puedes ver más detalles de splines en McElreath (2020)\nEn este caso, haremos expansión de entradas de las siguiente manera. Supongamos que tenemos la variable de distancia \\(d\\) que va de 0 a 750 cm, por ejemplo. Construimos entradas derivadas de la siguiente manera:\n\nlibrary(splines2)\nnudos <- c(25, 50, 100, 200, 400)\ndistancias <- seq(0, 750, 1)\nsplines_tbl <- iSpline(distancias, knots = nudos, \n Boundary.knots = c(0, 750), degree = 2, intercept = FALSE) |> \n as_tibble() |> \n mutate(d = distancias) |> \n pivot_longer(-d, names_to = \"spline\", values_to = \"valor\")\nggplot(splines_tbl) +\n geom_line(aes(x = d, y = valor, color = spline)) +\n geom_vline(xintercept = nudos, color = \"red\", linetype = 2) \n\n\n\n\n\n\n\n\nEsta gráfica muestra cómo para cada distancia \\(x\\) generamos valores \\(x_1,\\ldots, x_p\\) que son variables derivadas de \\(x\\). Podemos entonces obtener más flexibilidad hacer regresión en estas nuevas \\(p\\) variables en lugar de usar solamente \\(x\\). Por la elección de la base, obsérvese que siempre que \\(\\beta_1, \\ldots, \\beta_p\\) sean no negativos, entonces la función \\[\\alpha + \\beta_1 x_1 + \\cdots + \\beta_p x_p\\] será monótona no decreciente, que es lo que necesitamos para este problema.\nNuestra función generadora para este modelo puede ser:\n\nsimular_putts <- function(distancias, nudos) {\n # Simular intercepto\n alpha <- rnorm(1, 4, 2)\n # Simular coeficientes de splines\n beta <- - abs(rnorm(7, 0, 1.5))\n # Calcular splines para distancias dadas\n mat_splines <- splines2::iSpline(distancias, \n Boundary.knots = c(0, 750), knots = nudos, degree = 2, intercept = FALSE) \n # Calcular probabilidad de éxito con regresión logística\n p <- 1 / (1 + exp(- alpha - mat_splines %*% beta))\n tibble(y = rbinom(length(distancias), 1, p), p = p, d = distancias) |> \n select(d, p, y) |> \n mutate(alpha = alpha, beta = list(beta))\n}\n\n\nset.seed(8123)\ndistancias <- seq(1, 600, 5) |> rep(each = 5)\nsimular_putts(distancias, nudos) |> \n ggplot(aes(x = d, y = y)) +\n geom_jitter(height = 0.1) +\n labs(x = \"Distancia (cm)\", y = \"Éxito\") +\n geom_smooth(span = 1, se = FALSE)\n\n`geom_smooth()` using method = 'loess' and formula = 'y ~ x'\n\n\n\n\n\n\n\n\n\nY podemos hacer simulaciones a priori para entender nuestros supuestos:\n\nmap_df(1:100, \\(x) simular_putts(distancias, nudos) |> mutate(id = x)) |> \n ggplot(aes(x = d, y = p, group = id)) +\n geom_line(alpha = 0.2) +\n labs(x = \"Distancia (cm)\", y = \"Probabilidad de Éxito\")\n\n\n\n\n\n\n\n\nAhora construimos nuestro nuevo modelo en Stan, donde \\(x\\) será la matriz de splines (entradas derivadas como se explicó arriba):\n\n#! message: false\nlibrary(cmdstanr)\nmod_logistica_splines <- cmdstan_model(\"./src/golf-logistico-splines.stan\")\nprint(mod_logistica_splines)\n\ndata {\n int<lower=0> N;\n int<lower=0> p;\n array[N] int n;\n vector[N] d;\n matrix[N, p] x;\n array[N] int y;\n}\nparameters {\n real alpha;\n array[p] real<upper=0> beta;\n}\nmodel {\n for(i in 1:N){\n y[i] ~ binomial_logit(n[i], alpha + dot_product(x[i,], to_vector(beta)));\n }\n alpha ~ normal(4, 2);\n beta ~ normal(0, 1.5);\n}\n\n\n\nset.seed(1225)\nmat_splines <- splines2::iSpline(30.48 * datos_golf$x, \n Boundary.knots = c(0, 750), knots = nudos, degree = 2, intercept = FALSE) \najuste <- mod_logistica_splines$sample(\n data = list(N = nrow(datos_golf), p = ncol(mat_splines),\n d = 30.48 * datos_golf$x, \n x = mat_splines,\n y = datos_golf$y, n = datos_golf$n), \n refresh = 1000, init = 0.1, \n step_size = 0.1, adapt_delta = 0.99)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 1 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 1 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 1 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 1 finished in 2.9 seconds.\nChain 2 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 2 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 2 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 2 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 2 finished in 3.2 seconds.\nChain 3 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 3 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 3 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 3 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 3 finished in 3.2 seconds.\nChain 4 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 4 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 4 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 4 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 4 finished in 4.1 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 3.4 seconds.\nTotal execution time: 13.7 seconds.\n\n\nWarning: 236 of 4000 (6.0%) transitions hit the maximum treedepth limit of 10.\nSee https://mc-stan.org/misc/warnings for details.\n\nsims <- ajuste$draws(c(\"alpha\", \"beta\"), format = \"df\")\n\nresumen <- ajuste$summary()\n\n\nresumen\n\n# A tibble: 9 × 10\n variable mean median sd mad q5 q95 rhat ess_bulk\n <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>\n1 lp__ -2911. -2911. 2.22 2.13 -2916. -2908. 1.00 1193.\n2 alpha 4.88 4.80 0.939 0.934 3.49 6.57 1.00 1575.\n3 beta[1] -0.974 -0.803 0.748 0.743 -2.39 -0.0803 1.00 1876.\n4 beta[2] -1.24 -1.12 0.796 0.839 -2.70 -0.156 1.00 1665.\n5 beta[3] -1.91 -1.92 0.269 0.268 -2.35 -1.46 1.00 1548.\n6 beta[4] -1.03 -1.03 0.226 0.229 -1.40 -0.669 1.00 1501.\n7 beta[5] -1.23 -1.24 0.265 0.264 -1.63 -0.758 1.00 1650.\n8 beta[6] -0.403 -0.350 0.289 0.292 -0.949 -0.0319 1.00 1718.\n9 beta[7] -0.645 -0.529 0.521 0.485 -1.69 -0.0490 1.00 2091.\n# ℹ 1 more variable: ess_tail <dbl>\n\n\nAhora simulamos la posterior y la contrastamos con los datos:\n\nd <- 30.48 * seq(0, 20, 0.5)\nmat_splines_pred <- splines2::iSpline(30.48 * seq(0, 20, 0.5), \n Boundary.knots = c(0, 750), knots = nudos, degree = 2,\n intercept = FALSE) \nsims_2 <- sims |> group_by(.draw, .chain, .iteration) |> nest() \ngrafs <- purrr::map(sims_2$data, function(pars){\n pars <- as.numeric(pars)\n alpha <- pars[1]\n beta <- pars[2:8]\n p <- 1/(1 + exp(- alpha - mat_splines_pred %*% beta))\n tibble(p = as.numeric(p), d = d)\n})\nsims_graf_tbl <- sims_2 |> add_column(graf = grafs) |> select(-data) |> \n ungroup() |> \n slice_sample(n = 100) |> \n select(.draw, graf) |> \n unnest(graf) \n\n\nsims_graf_tbl |> \n ggplot(aes(x = d, y = p)) +\n geom_line(aes(group = .draw), alpha = 0.1) +\n labs(x = \"Distancia (cm)\", y = \"Probabilidad de Éxito\") +\n geom_point(data = resumen_golf, color = \"red\") +\n geom_linerange(data = resumen_golf, \n aes(ymin = p - 2 * sqrt(p * (1 - p) / n), \n ymax = p + 2 * sqrt(p * (1 - p) / n)),\n color = \"red\")\n\n\n\n\n\n\n\n\nEste modelo ajusta mejor, y puede ser usado para hacer comparaciones de probabilidad de éxito a diferentes distancias. Su defecto es que no es interpetable como nuestro modelo anterior (aprendemos poco sobre cómo funcionan los putts), y es considerablemente más difícil de ajustar.\nPuedes ver más de splines en McElreath (2020), y en Hastie, Tibshirani, y Friedman (2017). Puedes revisar también este caso de Stan que explica cómo utilizar splines de forma más general en Stan.\n\n\n\n\nGelman, Andrew, y Deborah Nolan. 2002. «A Probability Model for Golf Putting». Teaching Statistics 24 (septiembre): 93-95. https://doi.org/10.1111/1467-9639.00097.\n\n\nHastie, Trevor, Robert Tibshirani, y Jerome Friedman. 2017. The Elements of Statistical Learning. Springer Series en Statistics. Springer New York Inc. http://web.stanford.edu/~hastie/ElemStatLearn/.\n\n\nHolmes, Brian W. 1991. «Putting: How a golf ball and hole interact». American Journal of Physics 59 (2): 129-36. https://doi.org/10.1119/1.16592.\n\n\nMcElreath, R. 2020. Statistical Rethinking: A Bayesian Course with Examples in R and Stan. A Chapman & Hall libro. CRC Press. https://books.google.com.mx/books?id=Ie2vxQEACAAJ.\n\n\nPenner, Albert. 2002. «The physics of putting». Canadian Journal of Physics 80 (febrero): 83-96. https://doi.org/10.1139/p01-137.",
+ "text": "4.9 Modelos genéricos para ajustar curvas\nOtra posibilidad es utilizar un modelo más flexible creando variables derivadas de la distancia. En este caso, quizá podemos ajustar una curva que sea aceptable desde el punto de vista predictivo, pero no podremos aprender mucho acerca de cómo funciona la probabilidad de éxitos de los tiros de putts\n\n\n\n\n\n\nSplines y ajuste de curvas\n\n\n\nLos splines nos dan una manera estándar de ajustar curvas más flexibles, de tipo polinomial por tramos. Usualmente son numéricamente más conveniente que polinomios.\n\n\nAunque hay muchos tipos de splines (los más comunes son B-splines), para este problema consideraremos una base de splines cuadráticos que resultan en curvas monótonas (I-splines). Puedes ver más detalles de splines en McElreath (2020)\nEn este caso, haremos expansión de entradas de las siguiente manera. Supongamos que tenemos la variable de distancia \\(d\\) que va de 0 a 750 cm, por ejemplo. Construimos entradas derivadas de la siguiente manera:\n\nlibrary(splines2)\nnudos <- c(25, 50, 100, 200, 400)\ndistancias <- seq(0, 750, 1)\nsplines_tbl <- iSpline(distancias, knots = nudos, \n Boundary.knots = c(0, 750), degree = 2, intercept = FALSE) |> \n as_tibble() |> \n mutate(d = distancias) |> \n pivot_longer(-d, names_to = \"spline\", values_to = \"valor\")\nggplot(splines_tbl) +\n geom_line(aes(x = d, y = valor, color = spline)) +\n geom_vline(xintercept = nudos, color = \"red\", linetype = 2) \n\n\n\n\n\n\n\n\nEsta gráfica muestra cómo para cada distancia \\(x\\) generamos valores \\(x_1,\\ldots, x_p\\) que son variables derivadas de \\(x\\). Podemos entonces obtener más flexibilidad hacer regresión en estas nuevas \\(p\\) variables en lugar de usar solamente \\(x\\). Por la elección de la base, obsérvese que siempre que \\(\\beta_1, \\ldots, \\beta_p\\) sean no negativos, entonces la función \\[\\alpha + \\beta_1 x_1 + \\cdots + \\beta_p x_p\\] será monótona no decreciente, que es lo que necesitamos para este problema.\nNuestra función generadora para este modelo puede ser:\n\nsimular_putts <- function(distancias, nudos) {\n # Simular intercepto\n alpha <- rnorm(1, 4, 2)\n # Simular coeficientes de splines\n beta <- - abs(rnorm(7, 0, 1.5))\n # Calcular splines para distancias dadas\n mat_splines <- splines2::iSpline(distancias, \n Boundary.knots = c(0, 750), knots = nudos, degree = 2, intercept = FALSE) \n # Calcular probabilidad de éxito con regresión logística\n p <- 1 / (1 + exp(- alpha - mat_splines %*% beta))\n tibble(y = rbinom(length(distancias), 1, p), p = p, d = distancias) |> \n select(d, p, y) |> \n mutate(alpha = alpha, beta = list(beta))\n}\n\n\nset.seed(8123)\ndistancias <- seq(1, 600, 5) |> rep(each = 5)\nsimular_putts(distancias, nudos) |> \n ggplot(aes(x = d, y = y)) +\n geom_jitter(height = 0.1) +\n labs(x = \"Distancia (cm)\", y = \"Éxito\") +\n geom_smooth(span = 1, se = FALSE)\n\n`geom_smooth()` using method = 'loess' and formula = 'y ~ x'\n\n\n\n\n\n\n\n\n\nY podemos hacer simulaciones a priori para entender nuestros supuestos:\n\nmap_df(1:100, \\(x) simular_putts(distancias, nudos) |> mutate(id = x)) |> \n ggplot(aes(x = d, y = p, group = id)) +\n geom_line(alpha = 0.2) +\n labs(x = \"Distancia (cm)\", y = \"Probabilidad de Éxito\")\n\n\n\n\n\n\n\n\nAhora construimos nuestro nuevo modelo en Stan, donde \\(x\\) será la matriz de splines (entradas derivadas como se explicó arriba):\n\n#! message: false\nlibrary(cmdstanr)\nmod_logistica_splines <- cmdstan_model(\"./src/golf-logistico-splines.stan\")\nprint(mod_logistica_splines)\n\ndata {\n int<lower=0> N;\n int<lower=0> p;\n array[N] int n;\n vector[N] d;\n matrix[N, p] x;\n array[N] int y;\n}\nparameters {\n real alpha;\n array[p] real<upper=0> beta;\n}\nmodel {\n for(i in 1:N){\n y[i] ~ binomial_logit(n[i], alpha + dot_product(x[i,], to_vector(beta)));\n }\n alpha ~ normal(4, 2);\n beta ~ normal(0, 1.5);\n}\n\n\n\nset.seed(1225)\nmat_splines <- splines2::iSpline(30.48 * datos_golf$x, \n Boundary.knots = c(0, 750), knots = nudos, degree = 2, intercept = FALSE) \najuste <- mod_logistica_splines$sample(\n data = list(N = nrow(datos_golf), p = ncol(mat_splines),\n d = 30.48 * datos_golf$x, \n x = mat_splines,\n y = datos_golf$y, n = datos_golf$n), \n refresh = 1000, init = 0.1, \n step_size = 0.1, adapt_delta = 0.99)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 1 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 1 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 1 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 1 finished in 3.0 seconds.\nChain 2 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 2 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 2 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 2 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 2 finished in 3.2 seconds.\nChain 3 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 3 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 3 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 3 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 3 finished in 3.2 seconds.\nChain 4 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 4 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 4 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 4 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 4 finished in 4.1 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 3.4 seconds.\nTotal execution time: 13.8 seconds.\n\n\nWarning: 236 of 4000 (6.0%) transitions hit the maximum treedepth limit of 10.\nSee https://mc-stan.org/misc/warnings for details.\n\nsims <- ajuste$draws(c(\"alpha\", \"beta\"), format = \"df\")\n\nresumen <- ajuste$summary()\n\n\nresumen\n\n# A tibble: 9 × 10\n variable mean median sd mad q5 q95 rhat ess_bulk\n <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>\n1 lp__ -2911. -2911. 2.22 2.13 -2916. -2908. 1.00 1193.\n2 alpha 4.88 4.80 0.939 0.934 3.49 6.57 1.00 1575.\n3 beta[1] -0.974 -0.803 0.748 0.743 -2.39 -0.0803 1.00 1876.\n4 beta[2] -1.24 -1.12 0.796 0.839 -2.70 -0.156 1.00 1665.\n5 beta[3] -1.91 -1.92 0.269 0.268 -2.35 -1.46 1.00 1548.\n6 beta[4] -1.03 -1.03 0.226 0.229 -1.40 -0.669 1.00 1501.\n7 beta[5] -1.23 -1.24 0.265 0.264 -1.63 -0.758 1.00 1650.\n8 beta[6] -0.403 -0.350 0.289 0.292 -0.949 -0.0319 1.00 1718.\n9 beta[7] -0.645 -0.529 0.521 0.485 -1.69 -0.0490 1.00 2091.\n# ℹ 1 more variable: ess_tail <dbl>\n\n\nAhora simulamos la posterior y la contrastamos con los datos:\n\nd <- 30.48 * seq(0, 20, 0.5)\nmat_splines_pred <- splines2::iSpline(30.48 * seq(0, 20, 0.5), \n Boundary.knots = c(0, 750), knots = nudos, degree = 2,\n intercept = FALSE) \nsims_2 <- sims |> group_by(.draw, .chain, .iteration) |> nest() \ngrafs <- purrr::map(sims_2$data, function(pars){\n pars <- as.numeric(pars)\n alpha <- pars[1]\n beta <- pars[2:8]\n p <- 1/(1 + exp(- alpha - mat_splines_pred %*% beta))\n tibble(p = as.numeric(p), d = d)\n})\nsims_graf_tbl <- sims_2 |> add_column(graf = grafs) |> select(-data) |> \n ungroup() |> \n slice_sample(n = 100) |> \n select(.draw, graf) |> \n unnest(graf) \n\n\nsims_graf_tbl |> \n ggplot(aes(x = d, y = p)) +\n geom_line(aes(group = .draw), alpha = 0.1) +\n labs(x = \"Distancia (cm)\", y = \"Probabilidad de Éxito\") +\n geom_point(data = resumen_golf, color = \"red\") +\n geom_linerange(data = resumen_golf, \n aes(ymin = p - 2 * sqrt(p * (1 - p) / n), \n ymax = p + 2 * sqrt(p * (1 - p) / n)),\n color = \"red\")\n\n\n\n\n\n\n\n\nEste modelo ajusta mejor, y puede ser usado para hacer comparaciones de probabilidad de éxito a diferentes distancias. Su defecto es que no es interpetable como nuestro modelo anterior (aprendemos poco sobre cómo funcionan los putts), y es considerablemente más difícil de ajustar.\nPuedes ver más de splines en McElreath (2020), y en Hastie, Tibshirani, y Friedman (2017). Puedes revisar también este caso de Stan que explica cómo utilizar splines de forma más general en Stan.\n\n\n\n\nGelman, Andrew, y Deborah Nolan. 2002. «A Probability Model for Golf Putting». Teaching Statistics 24 (septiembre): 93-95. https://doi.org/10.1111/1467-9639.00097.\n\n\nHastie, Trevor, Robert Tibshirani, y Jerome Friedman. 2017. The Elements of Statistical Learning. Springer Series en Statistics. Springer New York Inc. http://web.stanford.edu/~hastie/ElemStatLearn/.\n\n\nHolmes, Brian W. 1991. «Putting: How a golf ball and hole interact». American Journal of Physics 59 (2): 129-36. https://doi.org/10.1119/1.16592.\n\n\nMcElreath, R. 2020. Statistical Rethinking: A Bayesian Course with Examples in R and Stan. A Chapman & Hall libro. CRC Press. https://books.google.com.mx/books?id=Ie2vxQEACAAJ.\n\n\nPenner, Albert. 2002. «The physics of putting». Canadian Journal of Physics 80 (febrero): 83-96. https://doi.org/10.1139/p01-137.",
"crumbs": [
"4Componentes de modelación 1"
]
@@ -314,7 +314,7 @@
"href": "05-dags.html#regla-del-producto-y-simulación",
"title": "5 Modelos gráficos y causalidad",
"section": "5.3 Regla del producto y simulación",
- "text": "5.3 Regla del producto y simulación\nEl orden del modelo gráfico también nos indica cómo simular las variables de la gráfica. Como cada modelo gráfico nos da una factorización de la conjunta, podemos utlizar esta para simular datos una vez que conocemos o estimamos las relaciones de dependencia directa. Empezamos con las variables exógenas (que no tienen padres) y vamos simulando hacia adelante.\n\nEjemplo\nEn nuestro ejemplo simulamos primero \\(X\\) y \\(D\\). A partir de \\(X\\) podemos simular \\(X_1\\) y \\(S_2\\), y a partir de \\(D\\), junto con \\(S_1\\) y \\(S_2\\), podemos simular \\(G\\). En nuestro ejemplo tendríamos\n\nsimular_juego <- function(N){\n x <- runif(N)\n d <- sample(c(\"lluvioso\",\"soleado\"), N, replace = TRUE, prob = c(0.3,0.7))\n s1 <- rbinom(N, 5, x)\n s2 <- rbinom(N, 5, x)\n g <- ifelse(d==\"lluvioso\", s1+s2, s1)\n tibble(x, d, s1, s2, g)\n}\nsimular_juego(5)\n\n# A tibble: 5 × 5\n x d s1 s2 g\n <dbl> <chr> <int> <int> <int>\n1 0.683 soleado 4 5 4\n2 0.0541 soleado 2 0 2\n3 0.810 soleado 4 5 4\n4 0.235 soleado 1 1 1\n5 0.120 soleado 0 2 0",
+ "text": "5.3 Regla del producto y simulación\nEl orden del modelo gráfico también nos indica cómo simular las variables de la gráfica. Como cada modelo gráfico nos da una factorización de la conjunta, podemos utlizar esta para simular datos una vez que conocemos o estimamos las relaciones de dependencia directa. Empezamos con las variables exógenas (que no tienen padres) y vamos simulando hacia adelante.\n\nEjemplo\nEn nuestro ejemplo simulamos primero \\(X\\) y \\(D\\). A partir de \\(X\\) podemos simular \\(X_1\\) y \\(S_2\\), y a partir de \\(D\\), junto con \\(S_1\\) y \\(S_2\\), podemos simular \\(G\\). En nuestro ejemplo tendríamos\n\nsimular_juego <- function(N){\n x <- runif(N)\n d <- sample(c(\"lluvioso\",\"soleado\"), N, replace = TRUE, prob = c(0.3,0.7))\n s1 <- rbinom(N, 5, x)\n s2 <- rbinom(N, 5, x)\n g <- ifelse(d==\"lluvioso\", s1+s2, s1)\n tibble(x, d, s1, s2, g)\n}\nsimular_juego(5)\n\n# A tibble: 5 × 5\n x d s1 s2 g\n <dbl> <chr> <int> <int> <int>\n1 0.715 soleado 4 2 4\n2 0.585 lluvioso 3 4 7\n3 0.285 soleado 2 3 2\n4 0.528 soleado 4 1 4\n5 0.410 soleado 2 3 2",
"crumbs": [
"5Modelos gráficos y causalidad"
]
@@ -334,7 +334,7 @@
"href": "05-dags.html#bifurcaciones-o-causa-común",
"title": "5 Modelos gráficos y causalidad",
"section": "5.5 Bifurcaciones o causa común",
- "text": "5.5 Bifurcaciones o causa común\nEn el siguiente ejemplo, llamamos a \\(Z\\) una causa que es común a \\(X\\) y \\(Y\\).\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n X\n Y\n Z\n edge [minlen = 3]\n Z -> X\n Z -> Y\n}\n\", width = 200, height = 50)\n\n\n\n\n\n\nEn este caso,\n\n\\(X\\) y \\(Y\\) tienen asociación\nSi condicionamos (o estratificamos) con \\(Z\\), entonces \\(X\\) y \\(Y\\) son condicionalmente independientes.\n\nEste tipo de estructura también se llama bifurcación, o decimos más tradicionalmente que \\(Z\\) es un confusor en esta gráfica. Variación en \\(Z\\) produce variación conjunta de \\(X\\) y \\(Y\\).\nPor ejemplo, podríamos encontrar que el uso de aspirina \\(X\\) está asociado a una mortalidad más alta \\(Y\\). Una causa común es enfermedad grave que produce dolor (\\(Z\\)). Sin embargo, si condicionamos a personas sanas, veríamos que no hay relación entre uso de aspirina y mortalidad, igualmente veríamos que entre las personas enfermas el uso de aspirina no les ayuda a vivir más tiempo.\nEn este caso, tenemos:\n\\[p(x, y, z) = p(z)p(x|z)p(y|z)\\] Y como el lado izquierdo es igual (en general) a \\(p(x,y|z)p(z)\\), obtenemos la independiencia condicional de \\(X\\) y \\(Y\\) dado \\(Z\\).\n\nEjemplo (simulación)\n\nrbern <- function(n, prob){\n rbinom(n, 1, prob = prob)\n} \nsimular_confusor <- function(n = 10){\n z <- rbern(n, p = 0.5) |> as.numeric()\n x <- rbern(n, p = z * 0.3 + (1 - z) * 0.8)\n y <- rbinom(n, 4, z * 0.9 + (1 - z) * 0.3)\n tibble(x, z, y)\n}\nsims_confusor <- simular_confusor(50000)\n\n\\(X\\) y \\(Y\\) están asociadas\n\nsims_confusor |> select(x, y) |> \n count(x, y) |> \n group_by(x) |> \n mutate(p_cond = n / sum(n)) |>\n select(x, y, p_cond) |> \nggplot(aes(x = y, y = p_cond, fill = factor(x))) +\n geom_col(position = \"dodge\") +\n labs(subtitle = \"Condicional de Y dada X\")\n\n\n\n\n\n\n\n\nLo cual lo vemos también si calculamos la correlación:\n\ncor(sims_confusor |> select(x,y)) |> round(3)\n\n x y\nx 1.000 -0.421\ny -0.421 1.000\n\n\nSin embargo, si condicionamos a \\(Z\\), que puede tomar los valores 0 o 1, vemos que \\(X\\) y \\(Y\\) son independientes, o dicho de otra manera, la condicional de \\(Y\\) dada \\(Z\\) y \\(X\\) sólo depende de \\(Z\\):\n\nsims_confusor |> \n count(x, y, z) |> \n group_by(x, z) |> \n mutate(p_cond = n / sum(n)) |>\n select(x, y, z, p_cond) |> \nggplot(aes(x = y, y = p_cond, fill = factor(x))) +\n geom_col(position = \"dodge\") + facet_wrap(~ z) +\n labs(subtitle = \"Condicional de Y dada X y Z\")\n\n\n\n\n\n\n\n\nUna consecuencia es por ejemplo que la correlación debe ser cero:\n\ncor(sims_confusor |> filter(z == 1) |> select(x,y)) |> round(3)\n\n x y\nx 1.000 -0.005\ny -0.005 1.000\n\ncor(sims_confusor |> filter(z == 0) |> select(x,y)) |> round(3)\n\n x y\nx 1.000 0.005\ny 0.005 1.000\n\n\nUn ejemplo con variables continuas podría ser como sigue:\n\nsimular_bifurcacion <- function(n = 10){\n z <- rbern(n, p = 0.5)\n x <- rnorm(n, 100 + 20 * z, 15)\n y <- rnorm(n, 100 + 30 * z, 20)\n tibble(x, z, y)\n}\nsims_bifurcacion <- simular_bifurcacion(5000)\n\n\\(X\\) y \\(Y\\) son dependientes (por ejemplo si vemos la media condicional de \\(Y\\) dado \\(X\\):\n\nggplot(sims_bifurcacion, aes(x = x, y = y, colour = z)) + \n geom_point(alpha = 0.2) +\n geom_smooth(span = 1, se = FALSE)\n\n\n\n\n\n\n\n\nSi condicionamos a \\(Z\\), no hay dependencia entre \\(X\\) y \\(Y\\)\n\nggplot(sims_bifurcacion, aes(x = x, y = y, colour = z, group = z)) + \n geom_point(alpha = 0.2) +\n geom_smooth(span = 2)\n\n\n\n\n\n\n\n\n\n\nEjemplo: matrimonio y divorcio\nEn este ejemplo de McElreath (2020), se muestra que regiones de Estados Unidos con tasas más altas de matrimonio también tienen tasas más altas de divorcio.\n\ndata(WaffleDivorce)\nWaffleDivorce |> \n ggplot(aes(x = Marriage, y = Divorce)) +\n geom_point() +\n geom_smooth(method = \"lm\")\n\n`geom_smooth()` using formula = 'y ~ x'\n\n\n\n\n\n\n\n\n\nAunque esta es una correlación clara, lo que nos interesa en este caso el efecto causal \\(M\\to D\\). Es importante notar que hay considerable variabilidad de la edad promedio al casarse a lo largo de los estados:\n\nWaffleDivorce |> \n ggplot(aes(sample = MedianAgeMarriage)) +\n geom_qq() +\n geom_qq_line()\n\n\n\n\n\n\n\n\nPara el modelo causal, tenemos que considerar las siguientes afirmaciones que no son muy difíciles de justificar:\n\nLa edad promedio al casarse de cada estado es un factor que influye en la tasa de divorcio (menor edad a casarse implica mayores tasas de divorcio, pues las parejas tienen más tiempo para divorciarse, porque la gente cambia más cuando es joven).\nAdicionalmente, si la gente tiende a casarse más joven, en cualquier momento hay más gente con probabilidad de casarse, por lo que esperaríamos que la edad al casarse también influye en la tasa de matrimonio.\n\nEsto implica que tenemos que considerar una causa común de la edad al casarse en nuestro diagrama causal:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n M\n D\n Edad\n edge [minlen = 3]\n Edad -> M\n Edad -> D\n M -> D\n{rank=same; M; D;}\n\n}\n\", width = 200, height = 50)\n\n\n\n\n\n\nPor la discusión de arriba, es claro que es necesario considerar la edad al casarse si queremos estimar el efecto de tasa de matrimonio en la tasa de divorcio. Es posible que la correlación entre estas dos tasas puede ser explicada solamente por la edad al casarse, y que en realidad al flecha \\(M\\to D\\) sea muy débil o inexistente.\nYa que tenemos este modelo causal básico, tendríamos que proponer un proceso generador, proponer un modelo estadístico, y probar nuestra estimación. Este paso nos lo saltaremos (ver sección anterior), aunque sigue siendo necesario.\nPor el momento recordemos que si condicionamos (se dice también estratificar) por edad al casarse, y no vemos relación condicional entre las dos tasas, la relación que vimos en los datos es factible que haya aparecido por la causa común que induce correlación. Una manera en que estratificamos o condicionamos a una variable continua en un modelo lineal, como sigue:\n\\[D_i\\sim N(\\mu_i, \\sigma)\\] donde \\[\\mu_i = \\alpha + \\beta_M M_i + \\beta_E Edad_i\\] ¿De qué manera estamos estratificando por edad en este ejemplo? Obsérvese que para cada Edad que fijemos, la relación entre \\(M\\) y \\(D\\) es:\n\\[\\mu_i = (\\alpha + \\beta_E Edad) + \\beta_M M_i \\] Cada valor de \\(E\\) produce una relación diferente entre \\(M\\) y \\(D\\) (en este caso particular, una recta diferente con distinta altura).\nAhora tenemos que poner iniciales para terminar nuestro modelo estadístico. En este punto poner iniciales informadas para estos coeficientes puede ser complicado (depende de cuánta demografía sabemos). Podemos usar un enfoque más simple, considerando las variables estandarizadas. De esta forma podemos poner iniciales más estándar. Utilizaremos\n\nescalar <- function(x){\n (x - mean(x))/sd(x)\n}\nWaffleDivorce <- WaffleDivorce |> \n mutate(Marriage_est = escalar(Marriage), \n Divorce_est = escalar(Divorce), \n MedianAgeMarriage_est = escalar(MedianAgeMarriage))\ndatos_lista <- list(\n N = nrow(WaffleDivorce),\n d_est = WaffleDivorce$Divorce_est, \n m_est = WaffleDivorce$Marriage_est, \n edad_est = WaffleDivorce$MedianAgeMarriage_est)\n\n\nmod_mat_div <- cmdstan_model(\"./src/matrimonio-divorcio-1.stan\")\nprint(mod_mat_div)\n\ndata {\n int<lower=0> N;\n vector[N] d_est;\n vector[N] m_est;\n vector[N] edad_est;\n}\n\nparameters {\n real alpha;\n real beta_M;\n real beta_E;\n real <lower=0> sigma;\n}\n\ntransformed parameters {\n vector[N] w_media;\n // determinístico dado parámetros\n w_media = alpha + beta_M * m_est + beta_E * edad_est;\n}\n\nmodel {\n // partes no determinísticas\n d_est ~ normal(w_media, sigma);\n alpha ~ normal(0, 1);\n beta_M ~ normal(0, 0.5);\n beta_E ~ normal(0, 0.5);\n sigma ~ normal(0, 1);\n}\n\ngenerated quantities {\n real dif;\n {\n //simulamos 50 estados\n int M = 50;\n array[M] real dif_sim;\n for(i in 1:M){\n real edad_sim_est = normal_rng(0, 1);\n // fijamos el valor de M en 0 y 1 para el modelo con do(M)\n real M_sim_0 = normal_rng(alpha * beta_M * 0 + beta_E * edad_sim_est, sigma);\n real M_sim_1 = normal_rng(alpha * beta_M * 1 + beta_E * edad_sim_est, sigma);\n dif_sim[i] = M_sim_1 - M_sim_0;\n }\n dif = mean(dif_sim);\n }\n\n}\n\n\n\nsims_mod <- mod_mat_div$sample(data = datos_lista, \n chains = 4, \n init = 0.1, step_size = 0.1,\n iter_warmup = 1000, \n iter_sampling = 1000,\n refresh = 0)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.1 seconds.\nChain 2 finished in 0.1 seconds.\nChain 3 finished in 0.1 seconds.\nChain 4 finished in 0.1 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.1 seconds.\nTotal execution time: 0.6 seconds.\n\n\n\nresumen <- sims_mod$summary(c(\"alpha\", \"beta_M\", \"beta_E\", \"sigma\"))\n\n\nresumen |> \n ggplot(aes(x = variable, y = mean, ymin = q5, ymax = q95)) +\n geom_hline(yintercept = 0, color = \"red\") +\n geom_point() +\n geom_linerange() +\n coord_flip()\n\n\n\n\n\n\n\n\nY el resultado que obtenemos es que no observamos un efecto considerable de las tasas de matrimonio en las tasas de divorcio, una vez que estratificamos por la causa común de edad de matrimonio. Este ejemplo es simple y podemos ver el efecto causal directo en un sólo coeficiente \\(\\beta_M\\), pero de todas formas haremos contrastes como hicimos en la parte anterior.\n\n\n5.5.1 Simulando intervenciones\nLa manera más directa de definir efecto causal, bajo nuestros supuestos causales, es a través de intervenciones (imaginarias o reales).\n\n\n\n\n\n\nNota\n\n\n\nEntendemos saber una causa como poder predecir correctamente las consecuencias de una intervención en el sistema generador de datos.\n\n\nEn nuestro caso, el diagrama de arriba muestra nuestro modelo causal. Si nosotros alteramos este proceso causal, interviniendo en la tasa de matrimonio, la distribución de matrimonio ya no depende de la Edad (pues está bajo nuestro control). Esto quiere decir que ahora consideramos el siguiente diagrama, en donde la nueva dependendencia del divorcio del matrimonio la escribiremos como \\(p(D|do(M))\\):\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n M\n D\n Edad\n edge [minlen = 3]\n Edad -> D\n M -> D\n{rank=same; M; D;}\n\n}\n\", width = 200, height = 50)\n\n\n\n\n\n\nEs decir, borramos todas las flechas que caen en \\(M\\) (pues la estamos interveniendo al valor que queramos), y luego simulando \\(D\\).\nEn nuestro ejemplo (ve el código de Stan de arriba, la parte de generated quantities) simularemos los 50 estados bajo dos intervenciones: todos tienen la tasa promedio de matrimonio vs. los 50 estados con tasa de matrimonio un error estándar por encima de la tasa promedio. Repetimos esta comparación sobre todas las simulaciones de la posterior:\n\nsims_tbl <- sims_mod$draws(format = \"df\") |> \n select(dif) \nsims_tbl |> summarize(\n q5 = quantile(dif, 0.05),\n q95 = quantile(dif, 0.95)\n)\n\n# A tibble: 1 × 2\n q5 q95\n <dbl> <dbl>\n1 -0.277 0.268\n\n\n\nggplot(sims_tbl, aes(x = dif)) +\n geom_histogram(bins = 50) +\n geom_vline(xintercept = 0, color = \"red\")\n\n\n\n\n\n\n\n\nEn este caso, vemos que el resultado de la intervención no tienen una tendencia clara hacia incrementar o disminuir la tasa de divorcio, aunque existe variabilidad por la incertidumbre que tenemos acerca de las relaciones modeladas.\n\n\n\n\n\n\nTip\n\n\n\nLa relación que vimos entre matrimonio y divorcio en nuestro ejemplo es probablemente producida por la causa común Edad, y no necesariamente es causal.\n\n\nFinalmente, antes de terminar sería apropiado hacer chequeos predictivos posteriores, pero por el momento los omitiremos para avanzar en los otros tipos de estructuras básicas en los DAGs.",
+ "text": "5.5 Bifurcaciones o causa común\nEn el siguiente ejemplo, llamamos a \\(Z\\) una causa que es común a \\(X\\) y \\(Y\\).\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n X\n Y\n Z\n edge [minlen = 3]\n Z -> X\n Z -> Y\n}\n\", width = 200, height = 50)\n\n\n\n\n\n\nEn este caso,\n\n\\(X\\) y \\(Y\\) tienen asociación\nSi condicionamos (o estratificamos) con \\(Z\\), entonces \\(X\\) y \\(Y\\) son condicionalmente independientes.\n\nEste tipo de estructura también se llama bifurcación, o decimos más tradicionalmente que \\(Z\\) es un confusor en esta gráfica. Variación en \\(Z\\) produce variación conjunta de \\(X\\) y \\(Y\\).\nPor ejemplo, podríamos encontrar que el uso de aspirina \\(X\\) está asociado a una mortalidad más alta \\(Y\\). Una causa común es enfermedad grave que produce dolor (\\(Z\\)). Sin embargo, si condicionamos a personas sanas, veríamos que no hay relación entre uso de aspirina y mortalidad, igualmente veríamos que entre las personas enfermas el uso de aspirina no les ayuda a vivir más tiempo.\nEn este caso, tenemos:\n\\[p(x, y, z) = p(z)p(x|z)p(y|z)\\] Y como el lado izquierdo es igual (en general) a \\(p(x,y|z)p(z)\\), obtenemos la independiencia condicional de \\(X\\) y \\(Y\\) dado \\(Z\\).\n\nEjemplo (simulación)\n\nrbern <- function(n, prob){\n rbinom(n, 1, prob = prob)\n} \nsimular_confusor <- function(n = 10){\n z <- rbern(n, p = 0.5) |> as.numeric()\n x <- rbern(n, p = z * 0.3 + (1 - z) * 0.8)\n y <- rbinom(n, 4, z * 0.9 + (1 - z) * 0.3)\n tibble(x, z, y)\n}\nsims_confusor <- simular_confusor(50000)\n\n\\(X\\) y \\(Y\\) están asociadas\n\nsims_confusor |> select(x, y) |> \n count(x, y) |> \n group_by(x) |> \n mutate(p_cond = n / sum(n)) |>\n select(x, y, p_cond) |> \nggplot(aes(x = y, y = p_cond, fill = factor(x))) +\n geom_col(position = \"dodge\") +\n labs(subtitle = \"Condicional de Y dada X\")\n\n\n\n\n\n\n\n\nLo cual lo vemos también si calculamos la correlación:\n\ncor(sims_confusor |> select(x,y)) |> round(3)\n\n x y\nx 1.000 -0.426\ny -0.426 1.000\n\n\nSin embargo, si condicionamos a \\(Z\\), que puede tomar los valores 0 o 1, vemos que \\(X\\) y \\(Y\\) son independientes, o dicho de otra manera, la condicional de \\(Y\\) dada \\(Z\\) y \\(X\\) sólo depende de \\(Z\\):\n\nsims_confusor |> \n count(x, y, z) |> \n group_by(x, z) |> \n mutate(p_cond = n / sum(n)) |>\n select(x, y, z, p_cond) |> \nggplot(aes(x = y, y = p_cond, fill = factor(x))) +\n geom_col(position = \"dodge\") + facet_wrap(~ z) +\n labs(subtitle = \"Condicional de Y dada X y Z\")\n\n\n\n\n\n\n\n\nUna consecuencia es por ejemplo que la correlación debe ser cero:\n\ncor(sims_confusor |> filter(z == 1) |> select(x,y)) |> round(3)\n\n x y\nx 1.000 -0.004\ny -0.004 1.000\n\ncor(sims_confusor |> filter(z == 0) |> select(x,y)) |> round(3)\n\n x y\nx 1.000 -0.014\ny -0.014 1.000\n\n\nUn ejemplo con variables continuas podría ser como sigue:\n\nsimular_bifurcacion <- function(n = 10){\n z <- rbern(n, p = 0.5)\n x <- rnorm(n, 100 + 20 * z, 15)\n y <- rnorm(n, 100 + 30 * z, 20)\n tibble(x, z, y)\n}\nsims_bifurcacion <- simular_bifurcacion(5000)\n\n\\(X\\) y \\(Y\\) son dependientes (por ejemplo si vemos la media condicional de \\(Y\\) dado \\(X\\):\n\nggplot(sims_bifurcacion, aes(x = x, y = y, colour = z)) + \n geom_point(alpha = 0.2) +\n geom_smooth(span = 1, se = FALSE)\n\n\n\n\n\n\n\n\nSi condicionamos a \\(Z\\), no hay dependencia entre \\(X\\) y \\(Y\\)\n\nggplot(sims_bifurcacion, aes(x = x, y = y, colour = z, group = z)) + \n geom_point(alpha = 0.2) +\n geom_smooth(span = 2)\n\n\n\n\n\n\n\n\n\n\nEjemplo: matrimonio y divorcio\nEn este ejemplo de McElreath (2020), se muestra que regiones de Estados Unidos con tasas más altas de matrimonio también tienen tasas más altas de divorcio.\n\ndata(WaffleDivorce)\nWaffleDivorce |> \n ggplot(aes(x = Marriage, y = Divorce)) +\n geom_point() +\n geom_smooth(method = \"lm\")\n\n`geom_smooth()` using formula = 'y ~ x'\n\n\n\n\n\n\n\n\n\nAunque esta es una correlación clara, lo que nos interesa en este caso el efecto causal \\(M\\to D\\). Es importante notar que hay considerable variabilidad de la edad promedio al casarse a lo largo de los estados:\n\nWaffleDivorce |> \n ggplot(aes(sample = MedianAgeMarriage)) +\n geom_qq() +\n geom_qq_line()\n\n\n\n\n\n\n\n\nPara el modelo causal, tenemos que considerar las siguientes afirmaciones que no son muy difíciles de justificar:\n\nLa edad promedio al casarse de cada estado es un factor que influye en la tasa de divorcio (menor edad a casarse implica mayores tasas de divorcio, pues las parejas tienen más tiempo para divorciarse, porque la gente cambia más cuando es joven).\nAdicionalmente, si la gente tiende a casarse más joven, en cualquier momento hay más gente con probabilidad de casarse, por lo que esperaríamos que la edad al casarse también influye en la tasa de matrimonio.\n\nEsto implica que tenemos que considerar una causa común de la edad al casarse en nuestro diagrama causal:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n M\n D\n Edad\n edge [minlen = 3]\n Edad -> M\n Edad -> D\n M -> D\n{rank=same; M; D;}\n\n}\n\", width = 200, height = 50)\n\n\n\n\n\n\nPor la discusión de arriba, es claro que es necesario considerar la edad al casarse si queremos estimar el efecto de tasa de matrimonio en la tasa de divorcio. Es posible que la correlación entre estas dos tasas puede ser explicada solamente por la edad al casarse, y que en realidad al flecha \\(M\\to D\\) sea muy débil o inexistente.\nYa que tenemos este modelo causal básico, tendríamos que proponer un proceso generador, proponer un modelo estadístico, y probar nuestra estimación. Este paso nos lo saltaremos (ver sección anterior), aunque sigue siendo necesario.\nPor el momento recordemos que si condicionamos (se dice también estratificar) por edad al casarse, y no vemos relación condicional entre las dos tasas, la relación que vimos en los datos es factible que haya aparecido por la causa común que induce correlación. Una manera en que estratificamos o condicionamos a una variable continua en un modelo lineal, como sigue:\n\\[D_i\\sim N(\\mu_i, \\sigma)\\] donde \\[\\mu_i = \\alpha + \\beta_M M_i + \\beta_E Edad_i\\] ¿De qué manera estamos estratificando por edad en este ejemplo? Obsérvese que para cada Edad que fijemos, la relación entre \\(M\\) y \\(D\\) es:\n\\[\\mu_i = (\\alpha + \\beta_E Edad) + \\beta_M M_i \\] Cada valor de \\(E\\) produce una relación diferente entre \\(M\\) y \\(D\\) (en este caso particular, una recta diferente con distinta altura).\nAhora tenemos que poner iniciales para terminar nuestro modelo estadístico. En este punto poner iniciales informadas para estos coeficientes puede ser complicado (depende de cuánta demografía sabemos). Podemos usar un enfoque más simple, considerando las variables estandarizadas. De esta forma podemos poner iniciales más estándar. Utilizaremos\n\nescalar <- function(x){\n (x - mean(x))/sd(x)\n}\nWaffleDivorce <- WaffleDivorce |> \n mutate(Marriage_est = escalar(Marriage), \n Divorce_est = escalar(Divorce), \n MedianAgeMarriage_est = escalar(MedianAgeMarriage))\ndatos_lista <- list(\n N = nrow(WaffleDivorce),\n d_est = WaffleDivorce$Divorce_est, \n m_est = WaffleDivorce$Marriage_est, \n edad_est = WaffleDivorce$MedianAgeMarriage_est)\n\n\nmod_mat_div <- cmdstan_model(\"./src/matrimonio-divorcio-1.stan\")\nprint(mod_mat_div)\n\ndata {\n int<lower=0> N;\n vector[N] d_est;\n vector[N] m_est;\n vector[N] edad_est;\n}\n\nparameters {\n real alpha;\n real beta_M;\n real beta_E;\n real <lower=0> sigma;\n}\n\ntransformed parameters {\n vector[N] w_media;\n // determinístico dado parámetros\n w_media = alpha + beta_M * m_est + beta_E * edad_est;\n}\n\nmodel {\n // partes no determinísticas\n d_est ~ normal(w_media, sigma);\n alpha ~ normal(0, 1);\n beta_M ~ normal(0, 0.5);\n beta_E ~ normal(0, 0.5);\n sigma ~ normal(0, 1);\n}\n\ngenerated quantities {\n real dif;\n {\n //simulamos 50 estados\n int M = 50;\n array[M] real dif_sim;\n for(i in 1:M){\n real edad_sim_est = normal_rng(0, 1);\n // fijamos el valor de M en 0 y 1 para el modelo con do(M)\n real M_sim_0 = normal_rng(alpha * beta_M * 0 + beta_E * edad_sim_est, sigma);\n real M_sim_1 = normal_rng(alpha * beta_M * 1 + beta_E * edad_sim_est, sigma);\n dif_sim[i] = M_sim_1 - M_sim_0;\n }\n dif = mean(dif_sim);\n }\n\n}\n\n\n\nsims_mod <- mod_mat_div$sample(data = datos_lista, \n chains = 4, \n init = 0.1, step_size = 0.1,\n iter_warmup = 1000, \n iter_sampling = 1000,\n refresh = 0)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 finished in 0.1 seconds.\nChain 2 finished in 0.1 seconds.\nChain 3 finished in 0.1 seconds.\nChain 4 finished in 0.1 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.1 seconds.\nTotal execution time: 0.6 seconds.\n\n\n\nresumen <- sims_mod$summary(c(\"alpha\", \"beta_M\", \"beta_E\", \"sigma\"))\n\n\nresumen |> \n ggplot(aes(x = variable, y = mean, ymin = q5, ymax = q95)) +\n geom_hline(yintercept = 0, color = \"red\") +\n geom_point() +\n geom_linerange() +\n coord_flip()\n\n\n\n\n\n\n\n\nY el resultado que obtenemos es que no observamos un efecto considerable de las tasas de matrimonio en las tasas de divorcio, una vez que estratificamos por la causa común de edad de matrimonio. Este ejemplo es simple y podemos ver el efecto causal directo en un sólo coeficiente \\(\\beta_M\\), pero de todas formas haremos contrastes como hicimos en la parte anterior.\n\n\n5.5.1 Simulando intervenciones\nLa manera más directa de definir efecto causal, bajo nuestros supuestos causales, es a través de intervenciones (imaginarias o reales).\n\n\n\n\n\n\nNota\n\n\n\nEntendemos saber una causa como poder predecir correctamente las consecuencias de una intervención en el sistema generador de datos.\n\n\nEn nuestro caso, el diagrama de arriba muestra nuestro modelo causal. Si nosotros alteramos este proceso causal, interviniendo en la tasa de matrimonio, la distribución de matrimonio ya no depende de la Edad (pues está bajo nuestro control). Esto quiere decir que ahora consideramos el siguiente diagrama, en donde la nueva dependendencia del divorcio del matrimonio la escribiremos como \\(p(D|do(M))\\):\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n M\n D\n Edad\n edge [minlen = 3]\n Edad -> D\n M -> D\n{rank=same; M; D;}\n\n}\n\", width = 200, height = 50)\n\n\n\n\n\n\nEs decir, borramos todas las flechas que caen en \\(M\\) (pues la estamos interveniendo al valor que queramos), y luego simulando \\(D\\).\nEn nuestro ejemplo (ve el código de Stan de arriba, la parte de generated quantities) simularemos los 50 estados bajo dos intervenciones: todos tienen la tasa promedio de matrimonio vs. los 50 estados con tasa de matrimonio un error estándar por encima de la tasa promedio. Repetimos esta comparación sobre todas las simulaciones de la posterior:\n\nsims_tbl <- sims_mod$draws(format = \"df\") |> \n select(dif) \nsims_tbl |> summarize(\n q5 = quantile(dif, 0.05),\n q95 = quantile(dif, 0.95)\n)\n\n# A tibble: 1 × 2\n q5 q95\n <dbl> <dbl>\n1 -0.269 0.280\n\n\n\nggplot(sims_tbl, aes(x = dif)) +\n geom_histogram(bins = 50) +\n geom_vline(xintercept = 0, color = \"red\")\n\n\n\n\n\n\n\n\nEn este caso, vemos que el resultado de la intervención no tienen una tendencia clara hacia incrementar o disminuir la tasa de divorcio, aunque existe variabilidad por la incertidumbre que tenemos acerca de las relaciones modeladas.\n\n\n\n\n\n\nTip\n\n\n\nLa relación que vimos entre matrimonio y divorcio en nuestro ejemplo es probablemente producida por la causa común Edad, y no necesariamente es causal.\n\n\nFinalmente, antes de terminar sería apropiado hacer chequeos predictivos posteriores, pero por el momento los omitiremos para avanzar en los otros tipos de estructuras básicas en los DAGs.",
"crumbs": [
"5Modelos gráficos y causalidad"
]
@@ -344,7 +344,7 @@
"href": "05-dags.html#cadenas-o-mediación",
"title": "5 Modelos gráficos y causalidad",
"section": "5.6 Cadenas o mediación",
- "text": "5.6 Cadenas o mediación\nEn este caso tenemos:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2, rankdir=LR]\n node [shape=plaintext]\n X\n Y\n Z\n edge [minlen = 3]\n X -> Z\n Z -> Y\n}\n\", width = 150, height = 20)\n\n\n\n\n\n\nEn este caso,\n\nExiste asociación entre \\(X\\) y \\(Y\\), pero no existe relación directa entre ellas. Decimos que \\(Z\\) es un mediador del efecto de \\(X\\) sobre \\(Y\\).\nSi condicionamos a un valor de \\(Z\\), \\(X\\) y \\(Y\\) son condicionalmente independientes.\n\nPodemos pensar en \\(Z\\) como un mediador del efecto de \\(X\\) sobre \\(Y\\). Si no permitimos que \\(Z\\) varíe, entonces la información de \\(X\\) no fluye a \\(Y\\).\nPor ejemplo, si \\(X\\) tomar o no una medicina para el dolor de cabeza, \\(Z\\) es dolor de cabeza y \\(Y\\) es bienestar general, \\(X\\) y \\(Y\\) están relacionadas. Sin embargo, si condicionamos a un valor fijo de dolor de cabeza, no hay relación entre tomar la medicina y bienestar general.\nEn términos de factorización, podemos checar la independencia condicional: como \\(p(x,y,z) = p(x)p(z|x)p(y|z)\\), entonces\n\\[p(x, y | z) = p(x,y,z) / p(z) = (p(x)(z|x)) (p(y|z) / p(z))\\] y vemos que el lado izquierdo se factoriza en una parte que sólo involucra a \\(x\\) y \\(z\\) y otro factor que sólo tiene a \\(y\\) y \\(z\\): no hay términos que incluyan conjuntamente a \\(x\\), \\(y\\) y \\(z\\). Podemos de cualquier forma continuar notando\n\\[p(x)p(z|x)/p(z) = p(x,z)/p(z) = p(x | z)\\] de modo que\n\\[p(x, y | z) = p(x|z) p(y|z) \\]\nY mostramos un ejemplo simulado:\n\nrbern <- function(n, prob){\n rbinom(n, 1, prob = prob)\n} \nsimular_mediador <- function(n = 10){\n x <- rbern(n, p = 0.5) |> as.numeric()\n z <- rbern(n, p = x * 0.8 + (1 - x) * 0.3)\n y <- rbinom(n, 2, z * 0.7 + (1 - z) * 0.5)\n tibble(x, z, y)\n}\nsims_mediador <- simular_mediador(50000)\n\n\\(X\\) y \\(Y\\) son dependientes:\n\nsims_mediador |> select(x, y) |> \n count(x, y) |> \n group_by(x) |> \n mutate(p_cond = n / sum(n)) |>\n select(x, y, p_cond) |> \nggplot(aes(x = y, y = p_cond, fill = factor(x))) +\n geom_col(position = \"dodge\") +\n labs(subtitle = \"Condicional de Y dada X\")\n\n\n\n\n\n\n\n\nSin embargo, si condicionamos a \\(Z\\), que puede tomar los valores 0 o 1:\n\nsims_mediador |> \n count(x, y, z) |> \n group_by(x, z) |> \n mutate(p_cond = n / sum(n)) |>\n select(x, y, z, p_cond) |> \nggplot(aes(x = y, y = p_cond, fill = factor(x))) +\n geom_col(position = \"dodge\") + facet_wrap(~ z) +\n labs(subtitle = \"Condicional de Y dada X y Z\")\n\n\n\n\n\n\n\n\nY vemos que la condicional de \\(Y\\) dada \\(Z\\) y \\(X\\) sólo depende de \\(Z\\). Una consecuencia es por ejemplo que la correlación debe ser cero:\n\ncor(sims_mediador |> filter(z == 1) |> select(x,y)) |> round(3)\n\n x y\nx 1.00 -0.01\ny -0.01 1.00\n\ncor(sims_mediador |> filter(z == 0) |> select(x,y)) |> round(3)\n\n x y\nx 1.000 0.006\ny 0.006 1.000\n\n\nPodemos también hacer un ejemplo continuo:\n\nsimular_mediador <- function(n = 10){\n x <- rnorm(n, 100, 10)\n prob <- 1 / (1 + exp(-(x - 100)/5))\n z <- rbern(n, p = prob)\n y <- rnorm(n, 100 + 30 * z, 15)\n tibble(x, z, y)\n}\nsims_mediador <- simular_mediador(2000)\n\n\\(X\\) y \\(Y\\) son dependientes (por ejemplo si vemos la media condicional de \\(Y\\) dado \\(X\\):\n\nggplot(sims_mediador, aes(x = x, y = y, colour = z)) + geom_point() +\n geom_smooth(span = 1, se = FALSE)\n\n`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = \"cs\")'\n\n\nWarning: The following aesthetics were dropped during statistical transformation: colour\nℹ This can happen when ggplot fails to infer the correct grouping structure in\n the data.\nℹ Did you forget to specify a `group` aesthetic or to convert a numerical\n variable into a factor?\n\n\n\n\n\n\n\n\n\nSi condicionamos a \\(Z\\), no hay dependencia entre \\(X\\) y \\(Y\\)\n\nggplot(sims_mediador, aes(x = x, y = y, colour = z, group = z)) + \n geom_point() +\n geom_smooth(span = 2)\n\n`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = \"cs\")'\n\n\n\n\n\n\n\n\n\nNótese que en este ejemplo sí hay un efecto causal de \\(X\\) sobre \\(Y\\), pero está mediado por otra variable \\(Z\\). Si condicionamos a \\(Z\\), no hay relación entre \\(X\\) y \\(Y\\). El análisis condicionado podría llevarnos a una conclusión errónea de que \\(X\\) no influye sobre \\(Y\\).\n\n\n\n\n\n\nTip\n\n\n\nNota que no existe una diferencia estadística entre una bifurcación y una cadena: en ambos casos, las variables \\(X\\) y \\(Y\\) están correlacionadas, y son independientes una vez que condicionamos o estratificamos por \\(Z\\). Sin embargo, su tratamiento en inferencia causal es muy diferente.\n\n\n\nSesgo post-tratamiento\nEn McElreath (2020) se discute que en algunos estudios experimentales, se estratifica por variables que son consecuencia del tratamiento. Esto induce sesgo post-tratamiento, lo cual puede llevar a equivocaciones en donde parece que el tratamiento no tiene efecto cuando sí lo tiene. Incluso bajo condiciones de experimento (donde el tratamiento es asignado al azar) estratificar por mediadores es una mala idea. Ver más en McElreath (2020), donde por ejemplo cita una fuente que en estudios experimentales de Ciencia Política, casi la mitad de ellos sufre de este tipo de sesgo por estratificación por mediadores.\n\n\nEjemplo: Burks\nEste ejemplo es de Pearl y Mackenzie (2018). En 1926 Burks recolectó datos sobre qué tanto podría esperarse que la inteligencia de padres se hereda a los hijos (medido según una prueba de IQ). Construyó un diagrama parecido al de abajo:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape = circle]\n U\n node [shape=plaintext]\n edge [minlen = 3]\n IntPadres -> NSE\n NSE -> IntHijos\n U -> NSE\n U -> IntHijos\n IntPadres -> IntHijos\n{rank = same; U}\n}\n\")\n\n\n\n\n\n\nComo el NSE es del hogar (una medida general de estatus social), se consideró en principio como una variable pre-tratamiento a la inteligencia de los niños por la que tradicionalmente se controlaba. Burks notó que hacer esto tenía no era apropiado, pues tiene como consecuencia cortar parte del efecto total de la inteligencia sobre el la inteligencia de los hijos. En otras palabras: la inteligencia de los padres hace más probable mejor NSE, y mejor NSE presenta mejores condiciones de desarrollo para sus hijos. Estatificar por esta variable bloquea este efecto.\nAdicionalmente, como veremos, condicionar a NSE abre un camino no causal entre Inteligencia de Padres e Hijos.",
+ "text": "5.6 Cadenas o mediación\nEn este caso tenemos:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2, rankdir=LR]\n node [shape=plaintext]\n X\n Y\n Z\n edge [minlen = 3]\n X -> Z\n Z -> Y\n}\n\", width = 150, height = 20)\n\n\n\n\n\n\nEn este caso,\n\nExiste asociación entre \\(X\\) y \\(Y\\), pero no existe relación directa entre ellas. Decimos que \\(Z\\) es un mediador del efecto de \\(X\\) sobre \\(Y\\).\nSi condicionamos a un valor de \\(Z\\), \\(X\\) y \\(Y\\) son condicionalmente independientes.\n\nPodemos pensar en \\(Z\\) como un mediador del efecto de \\(X\\) sobre \\(Y\\). Si no permitimos que \\(Z\\) varíe, entonces la información de \\(X\\) no fluye a \\(Y\\).\nPor ejemplo, si \\(X\\) tomar o no una medicina para el dolor de cabeza, \\(Z\\) es dolor de cabeza y \\(Y\\) es bienestar general, \\(X\\) y \\(Y\\) están relacionadas. Sin embargo, si condicionamos a un valor fijo de dolor de cabeza, no hay relación entre tomar la medicina y bienestar general.\nEn términos de factorización, podemos checar la independencia condicional: como \\(p(x,y,z) = p(x)p(z|x)p(y|z)\\), entonces\n\\[p(x, y | z) = p(x,y,z) / p(z) = (p(x)(z|x)) (p(y|z) / p(z))\\] y vemos que el lado izquierdo se factoriza en una parte que sólo involucra a \\(x\\) y \\(z\\) y otro factor que sólo tiene a \\(y\\) y \\(z\\): no hay términos que incluyan conjuntamente a \\(x\\), \\(y\\) y \\(z\\). Podemos de cualquier forma continuar notando\n\\[p(x)p(z|x)/p(z) = p(x,z)/p(z) = p(x | z)\\] de modo que\n\\[p(x, y | z) = p(x|z) p(y|z) \\]\nY mostramos un ejemplo simulado:\n\nrbern <- function(n, prob){\n rbinom(n, 1, prob = prob)\n} \nsimular_mediador <- function(n = 10){\n x <- rbern(n, p = 0.5) |> as.numeric()\n z <- rbern(n, p = x * 0.8 + (1 - x) * 0.3)\n y <- rbinom(n, 2, z * 0.7 + (1 - z) * 0.5)\n tibble(x, z, y)\n}\nsims_mediador <- simular_mediador(50000)\n\n\\(X\\) y \\(Y\\) son dependientes:\n\nsims_mediador |> select(x, y) |> \n count(x, y) |> \n group_by(x) |> \n mutate(p_cond = n / sum(n)) |>\n select(x, y, p_cond) |> \nggplot(aes(x = y, y = p_cond, fill = factor(x))) +\n geom_col(position = \"dodge\") +\n labs(subtitle = \"Condicional de Y dada X\")\n\n\n\n\n\n\n\n\nSin embargo, si condicionamos a \\(Z\\), que puede tomar los valores 0 o 1:\n\nsims_mediador |> \n count(x, y, z) |> \n group_by(x, z) |> \n mutate(p_cond = n / sum(n)) |>\n select(x, y, z, p_cond) |> \nggplot(aes(x = y, y = p_cond, fill = factor(x))) +\n geom_col(position = \"dodge\") + facet_wrap(~ z) +\n labs(subtitle = \"Condicional de Y dada X y Z\")\n\n\n\n\n\n\n\n\nY vemos que la condicional de \\(Y\\) dada \\(Z\\) y \\(X\\) sólo depende de \\(Z\\). Una consecuencia es por ejemplo que la correlación debe ser cero:\n\ncor(sims_mediador |> filter(z == 1) |> select(x,y)) |> round(3)\n\n x y\nx 1.000 -0.001\ny -0.001 1.000\n\ncor(sims_mediador |> filter(z == 0) |> select(x,y)) |> round(3)\n\n x y\nx 1.000 -0.004\ny -0.004 1.000\n\n\nPodemos también hacer un ejemplo continuo:\n\nsimular_mediador <- function(n = 10){\n x <- rnorm(n, 100, 10)\n prob <- 1 / (1 + exp(-(x - 100)/5))\n z <- rbern(n, p = prob)\n y <- rnorm(n, 100 + 30 * z, 15)\n tibble(x, z, y)\n}\nsims_mediador <- simular_mediador(2000)\n\n\\(X\\) y \\(Y\\) son dependientes (por ejemplo si vemos la media condicional de \\(Y\\) dado \\(X\\):\n\nggplot(sims_mediador, aes(x = x, y = y, colour = z)) + geom_point() +\n geom_smooth(span = 1, se = FALSE)\n\n`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = \"cs\")'\n\n\nWarning: The following aesthetics were dropped during statistical transformation: colour\nℹ This can happen when ggplot fails to infer the correct grouping structure in\n the data.\nℹ Did you forget to specify a `group` aesthetic or to convert a numerical\n variable into a factor?\n\n\n\n\n\n\n\n\n\nSi condicionamos a \\(Z\\), no hay dependencia entre \\(X\\) y \\(Y\\)\n\nggplot(sims_mediador, aes(x = x, y = y, colour = z, group = z)) + \n geom_point() +\n geom_smooth(span = 2)\n\n`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = \"cs\")'\n\n\n\n\n\n\n\n\n\nNótese que en este ejemplo sí hay un efecto causal de \\(X\\) sobre \\(Y\\), pero está mediado por otra variable \\(Z\\). Si condicionamos a \\(Z\\), no hay relación entre \\(X\\) y \\(Y\\). El análisis condicionado podría llevarnos a una conclusión errónea de que \\(X\\) no influye sobre \\(Y\\).\n\n\n\n\n\n\nTip\n\n\n\nNota que no existe una diferencia estadística entre una bifurcación y una cadena: en ambos casos, las variables \\(X\\) y \\(Y\\) están correlacionadas, y son independientes una vez que condicionamos o estratificamos por \\(Z\\). Sin embargo, su tratamiento en inferencia causal es muy diferente.\n\n\n\nSesgo post-tratamiento\nEn McElreath (2020) se discute que en algunos estudios experimentales, se estratifica por variables que son consecuencia del tratamiento. Esto induce sesgo post-tratamiento, lo cual puede llevar a equivocaciones en donde parece que el tratamiento no tiene efecto cuando sí lo tiene. Incluso bajo condiciones de experimento (donde el tratamiento es asignado al azar) estratificar por mediadores es una mala idea. Ver más en McElreath (2020), donde por ejemplo cita una fuente que en estudios experimentales de Ciencia Política, casi la mitad de ellos sufre de este tipo de sesgo por estratificación por mediadores.\n\n\nEjemplo: Burks\nEste ejemplo es de Pearl y Mackenzie (2018). En 1926 Burks recolectó datos sobre qué tanto podría esperarse que la inteligencia de padres se hereda a los hijos (medido según una prueba de IQ). Construyó un diagrama parecido al de abajo:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape = circle]\n U\n node [shape=plaintext]\n edge [minlen = 3]\n IntPadres -> NSE\n NSE -> IntHijos\n U -> NSE\n U -> IntHijos\n IntPadres -> IntHijos\n{rank = same; U}\n}\n\")\n\n\n\n\n\n\nComo el NSE es del hogar (una medida general de estatus social), se consideró en principio como una variable pre-tratamiento a la inteligencia de los niños por la que tradicionalmente se controlaba. Burks notó que hacer esto tenía no era apropiado, pues tiene como consecuencia cortar parte del efecto total de la inteligencia sobre el la inteligencia de los hijos. En otras palabras: la inteligencia de los padres hace más probable mejor NSE, y mejor NSE presenta mejores condiciones de desarrollo para sus hijos. Estatificar por esta variable bloquea este efecto.\nAdicionalmente, como veremos, condicionar a NSE abre un camino no causal entre Inteligencia de Padres e Hijos.",
"crumbs": [
"5Modelos gráficos y causalidad"
]
@@ -354,7 +354,7 @@
"href": "05-dags.html#colisionador-o-causas-alternativas",
"title": "5 Modelos gráficos y causalidad",
"section": "5.7 Colisionador o causas alternativas",
- "text": "5.7 Colisionador o causas alternativas\nEn este caso, a \\(Z\\) también le llamamos un colisionador. Este es el caso que puede ser más difícil de entender en un principio. Consiste de la siguiente estructura:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n X\n Y\n Z\n edge [minlen = 3]\n X -> Z\n Y -> Z\n}\n\", width = 200, height = 50)\n\n\n\n\n\n\n\nEn este caso \\(X\\) y \\(Y\\) son independientes. Tanto \\(X\\) como \\(Y\\) influyen en \\(Z\\).\nSin embargo, si condicionamos a \\(Z\\) entonces \\(X\\) y \\(Y\\) están asociados.\n\nPor ejemplo, si observamos que el pasto está mojado, entonces saber que no llovió implica que probablemente se encendieron los aspersores.\nComo la conjunta se factoriza como:\n\\[p(x,y,z) = p(x)p(y)p(z|x,y)\\] Entonces integrando sobre \\(Z\\):\n\\[p(x,y) = \\int p(x,y,z)dz = p(x)p(y)\\int p(z|x,y)\\, dz\\] pero \\(p(z|x,y)\\) integra uno porque es una densidad, de forma que \\(x\\) y \\(y\\) son independientes.\nMostramos un ejemplo simulado:\n\nsimular_colisionador <- function(n = 10){\n x <- rbern(n, 0.5) \n y <- rbinom(n, 2, 0.7)\n z <- rbern(n, p = 0.1 + 0.7 * x * (y > 1)) \n tibble(x, z, y)\n}\nsims_colisionador <- simular_colisionador(50000)\n\n\\(X\\) y \\(Y\\) son independientes:\n\nsims_colisionador|> select(x, y) |> \n count(x, y) |> \n group_by(x) |> \n mutate(p_cond = n / sum(n)) |>\n select(x, y, p_cond) |> \nggplot(aes(x = y, y = p_cond, fill = factor(x))) +\n geom_col(position = \"dodge\") +\n labs(subtitle = \"Condicional de Y dada X\")\n\n\n\n\n\n\n\ncor(sims_colisionador |> select(x,y))\n\n x y\nx 1.0000000000 0.0004289729\ny 0.0004289729 1.0000000000\n\n\nSin embargo, si condicionamos a \\(Z\\), que puede tomar los valores 0 o 1:\n\nsims_colisionador |> \n count(x, y, z) |> \n group_by(x, z) |> \n mutate(p_cond = n / sum(n)) |>\n select(x, y, z, p_cond) |> \nggplot(aes(x = y, y = p_cond, fill = factor(x))) +\n geom_col(position = \"dodge\") + facet_wrap(~ z) +\n labs(subtitle = \"Condicional de Y dada X y Z\")\n\n\n\n\n\n\n\n\nY vemos que la condicional de \\(Y\\) dada \\(Z\\) y \\(X\\) depende de \\(X\\) y de \\(Z\\).\nLas correlaciones condicionales, por ejemplo, no son cero:\n\nprint(\"Dado Z = 0\")\n\n[1] \"Dado Z = 0\"\n\ncor(sims_colisionador |> filter(z == 0) |> select(x,y)) |> round(3)\n\n x y\nx 1.000 -0.274\ny -0.274 1.000\n\nprint(\"Dado Z = 1\")\n\n[1] \"Dado Z = 1\"\n\ncor(sims_colisionador |> filter(z == 1) |> select(x,y)) |> round(3)\n\n x y\nx 1.000 0.362\ny 0.362 1.000\n\n\nOtro ejemplo con variables continuas:\n\nsimular_colisionador_2 <- function(n = 10){\n x <- rnorm(n, 100, 20) \n y <- rnorm(n, 100, 20)\n z <- rbern(n, p = 0.92 * ((x + y) > 220) + 0.05) \n tibble(x, z, y)\n}\nsims_colisionador <- simular_colisionador_2(1000)\n\n\\(X\\) y \\(Y\\) son independientes:\n\nggplot(sims_colisionador, aes(x = x, y = y)) + geom_point()\n\n\n\n\n\n\n\n\nSin embargo, si condicionamos a un valor de \\(Z\\), \\(X\\) y \\(Y\\) ya no son independientes:\n\nggplot(sims_colisionador, aes(x = x, y = y, group = z, colour = factor(z))) + \n geom_point() + geom_smooth(method = \"lm\", se = FALSE) \n\n`geom_smooth()` using formula = 'y ~ x'\n\n\n\n\n\n\n\n\n\nY vemos que condicional a \\(Z\\), \\(X\\) y \\(Y\\) están correlacionadas, aunque no hay relación causal entre \\(X\\) y \\(Y\\).\n\n5.7.1 Ejemplos de colisionadores\nExisten muchos ejemplos de colisionadores en análisis de datos. Algunos ejemplos se deben a sesgo de selección (puedes dibujar diagramas para cada uno de estos):\n\nPodemos observar correlaciones entre habilidades que en realidad son independientes si observamos muestras de estudiantes seleccionados por un examen de admisión (por ejemplo, para entrar es necesario tener alta habilidad atlética y/o alta habilidad académica).\nEntre los artículos científicos publicados (ver McElreath (2020)), aquellos que son más tomados por las noticias son los menos confiables. Esta correlación puede aparecer aunque no exista relación en proyectos científicos entre confiabilidad e interés de los medios, pues lo que se fondea o publica puede tener dos razones: ser trabajo muy confiable, o ser trabajo que “está de moda” o atrae la atención de los medios.\n\nPero también puede ser consecuencia de condicionar a variables endógenos (que resultan ser colisionadores), y ocurren como parte del procesamiento o construcción de modelos. Un ejemplo interesante de McElreath (2020) es el siguiente:\n\nNos interesa saber si la edad influye en la felicidad o bienestar de las personas.\nAlgún investigador puede pensar que es necesario controlar por sí las personas están casadas o no, por ejemplo, para “quitar” ese efecto o algo así.\nEsto puede ser mala idea si consideramos que un diagrama apropiado puede ser \\(F \\rightarrow Matrim \\leftarrow Edad\\), que se basa en las observaciones de que personas más felices generalmente tienen mayor posibilidad de casarse, y también conforme pasa el tiempo, hay más oportunidades para casarse.\nEsto induce una correlación no causal entre edad y felicidad dentro de los grupos de casados y no casados, y puede llevar a conclusiones incorrectas.",
+ "text": "5.7 Colisionador o causas alternativas\nEn este caso, a \\(Z\\) también le llamamos un colisionador. Este es el caso que puede ser más difícil de entender en un principio. Consiste de la siguiente estructura:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n X\n Y\n Z\n edge [minlen = 3]\n X -> Z\n Y -> Z\n}\n\", width = 200, height = 50)\n\n\n\n\n\n\n\nEn este caso \\(X\\) y \\(Y\\) son independientes. Tanto \\(X\\) como \\(Y\\) influyen en \\(Z\\).\nSin embargo, si condicionamos a \\(Z\\) entonces \\(X\\) y \\(Y\\) están asociados.\n\nPor ejemplo, si observamos que el pasto está mojado, entonces saber que no llovió implica que probablemente se encendieron los aspersores.\nComo la conjunta se factoriza como:\n\\[p(x,y,z) = p(x)p(y)p(z|x,y)\\] Entonces integrando sobre \\(Z\\):\n\\[p(x,y) = \\int p(x,y,z)dz = p(x)p(y)\\int p(z|x,y)\\, dz\\] pero \\(p(z|x,y)\\) integra uno porque es una densidad, de forma que \\(x\\) y \\(y\\) son independientes.\nMostramos un ejemplo simulado:\n\nsimular_colisionador <- function(n = 10){\n x <- rbern(n, 0.5) \n y <- rbinom(n, 2, 0.7)\n z <- rbern(n, p = 0.1 + 0.7 * x * (y > 1)) \n tibble(x, z, y)\n}\nsims_colisionador <- simular_colisionador(50000)\n\n\\(X\\) y \\(Y\\) son independientes:\n\nsims_colisionador|> select(x, y) |> \n count(x, y) |> \n group_by(x) |> \n mutate(p_cond = n / sum(n)) |>\n select(x, y, p_cond) |> \nggplot(aes(x = y, y = p_cond, fill = factor(x))) +\n geom_col(position = \"dodge\") +\n labs(subtitle = \"Condicional de Y dada X\")\n\n\n\n\n\n\n\ncor(sims_colisionador |> select(x,y))\n\n x y\nx 1.00000e+00 3.52394e-05\ny 3.52394e-05 1.00000e+00\n\n\nSin embargo, si condicionamos a \\(Z\\), que puede tomar los valores 0 o 1:\n\nsims_colisionador |> \n count(x, y, z) |> \n group_by(x, z) |> \n mutate(p_cond = n / sum(n)) |>\n select(x, y, z, p_cond) |> \nggplot(aes(x = y, y = p_cond, fill = factor(x))) +\n geom_col(position = \"dodge\") + facet_wrap(~ z) +\n labs(subtitle = \"Condicional de Y dada X y Z\")\n\n\n\n\n\n\n\n\nY vemos que la condicional de \\(Y\\) dada \\(Z\\) y \\(X\\) depende de \\(X\\) y de \\(Z\\).\nLas correlaciones condicionales, por ejemplo, no son cero:\n\nprint(\"Dado Z = 0\")\n\n[1] \"Dado Z = 0\"\n\ncor(sims_colisionador |> filter(z == 0) |> select(x,y)) |> round(3)\n\n x y\nx 1.000 -0.273\ny -0.273 1.000\n\nprint(\"Dado Z = 1\")\n\n[1] \"Dado Z = 1\"\n\ncor(sims_colisionador |> filter(z == 1) |> select(x,y)) |> round(3)\n\n x y\nx 1.000 0.348\ny 0.348 1.000\n\n\nOtro ejemplo con variables continuas:\n\nsimular_colisionador_2 <- function(n = 10){\n x <- rnorm(n, 100, 20) \n y <- rnorm(n, 100, 20)\n z <- rbern(n, p = 0.92 * ((x + y) > 220) + 0.05) \n tibble(x, z, y)\n}\nsims_colisionador <- simular_colisionador_2(1000)\n\n\\(X\\) y \\(Y\\) son independientes:\n\nggplot(sims_colisionador, aes(x = x, y = y)) + geom_point()\n\n\n\n\n\n\n\n\nSin embargo, si condicionamos a un valor de \\(Z\\), \\(X\\) y \\(Y\\) ya no son independientes:\n\nggplot(sims_colisionador, aes(x = x, y = y, group = z, colour = factor(z))) + \n geom_point() + geom_smooth(method = \"lm\", se = FALSE) \n\n`geom_smooth()` using formula = 'y ~ x'\n\n\n\n\n\n\n\n\n\nY vemos que condicional a \\(Z\\), \\(X\\) y \\(Y\\) están correlacionadas, aunque no hay relación causal entre \\(X\\) y \\(Y\\).\n\n5.7.1 Ejemplos de colisionadores\nExisten muchos ejemplos de colisionadores en análisis de datos. Algunos ejemplos se deben a sesgo de selección (puedes dibujar diagramas para cada uno de estos):\n\nPodemos observar correlaciones entre habilidades que en realidad son independientes si observamos muestras de estudiantes seleccionados por un examen de admisión (por ejemplo, para entrar es necesario tener alta habilidad atlética y/o alta habilidad académica).\nEntre los artículos científicos publicados (ver McElreath (2020)), aquellos que son más tomados por las noticias son los menos confiables. Esta correlación puede aparecer aunque no exista relación en proyectos científicos entre confiabilidad e interés de los medios, pues lo que se fondea o publica puede tener dos razones: ser trabajo muy confiable, o ser trabajo que “está de moda” o atrae la atención de los medios.\n\nPero también puede ser consecuencia de condicionar a variables endógenos (que resultan ser colisionadores), y ocurren como parte del procesamiento o construcción de modelos. Un ejemplo interesante de McElreath (2020) es el siguiente:\n\nNos interesa saber si la edad influye en la felicidad o bienestar de las personas.\nAlgún investigador puede pensar que es necesario controlar por sí las personas están casadas o no, por ejemplo, para “quitar” ese efecto o algo así.\nEsto puede ser mala idea si consideramos que un diagrama apropiado puede ser \\(F \\rightarrow Matrim \\leftarrow Edad\\), que se basa en las observaciones de que personas más felices generalmente tienen mayor posibilidad de casarse, y también conforme pasa el tiempo, hay más oportunidades para casarse.\nEsto induce una correlación no causal entre edad y felicidad dentro de los grupos de casados y no casados, y puede llevar a conclusiones incorrectas.",
"crumbs": [
"5Modelos gráficos y causalidad"
]
@@ -364,7 +364,7 @@
"href": "05-dags.html#razonamiento-de-descendientes",
"title": "5 Modelos gráficos y causalidad",
"section": "5.8 Razonamiento de descendientes",
- "text": "5.8 Razonamiento de descendientes\nCondicionar a un descendiente puede entenderse como “condicionar parcialmente” o “débilmente” a los padres de ese descendiente.\nPor ejemplo, condicionar a un colisionador también produce dependencias condicionales:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n X\n Y\n Z\n A\n edge [minlen = 3]\n X -> Z\n Y -> Z\n Z -> A\n}\n\", width = 200, height = 50)\n\n\n\n\n\n\nEn este caso,\n\n\\(X\\) y \\(Y\\) son independientes\n\\(X\\) y \\(Y\\) son dependientes si condicionamos a \\(A\\).\n\nDependiendo de la naturaleza de la asociación entre el colisionador \\(Z\\) y su descendiente \\(A\\), esta dependencia puede ser más fuerte o más débil.\nPor ejemplo, en nuestro ejemplo donde el pasto mojado es un colisionador entre cuánta agua dieron los aspersores y cuánta lluvia cayó, un descendiente del pasto mojado es el estado de las plantas del jardín. Aunque los aspersores trabajan independientemente de la lluvia, si observamos que las plantas se secaron entonces lluvia y aspersores están correlacionados: por ejemplo, si noto que los aspersores están descompuestos, entonces concluimos que no hubo lluvia.\n\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n X [label = lluvia]\n Y [label = aspersores]\n Z [label = humedad]\n A [label = plantas]\n edge [minlen = 3]\n X -> Z\n Y -> Z\n Z -> A\n}\n\", width = 200, height = 50)\n\n\n\n\n\n\nEjemplo\n\nsimular_desc <- function(n = 10){\n x <- rbern(n, 0.5) \n y <- rbinom(n, 2, 0.7)\n z <- rbern(n, p = 0.1 + 0.7 * x * (y > 1)) \n a <- rbern(n, p = 0.5 + 0.5 * z)\n tibble(x, z, y, a)\n}\nsims_colisionador <- simular_desc(50000)\n# No hay correlación\ncor(sims_colisionador$x, sims_colisionador$y)\n\n[1] 0.0005184629\n\n\nSin embargo,\n\ncor(sims_colisionador |> filter(a ==0) |> select(x,y))\n\n x y\nx 1.0000000 -0.2758999\ny -0.2758999 1.0000000\n\n\n\ncor(sims_colisionador |> filter(a ==1) |> select(x,y))\n\n x y\nx 1.000000 0.111952\ny 0.111952 1.000000\n\n\n\n\n5.8.1 Ejemplo: dependencias de colisionador\nVerificamos que en nuestro modelo de Santa Clara, efectivamente nuestro modelo no implica ninguna dependencia no condicional entre sensibilidad de la prueba y prevalencia. Eso debería ser claro de la simulación, pero de todas formas lo checamos\n\nlibrary(cmdstanr)\nmod_sc <- cmdstan_model(\"./src/sclara.stan\")\nprint(mod_sc)\n\ndata {\n int<lower=0> N;\n int<lower=0> n;\n int<lower=0> kit_pos;\n int<lower=0> n_kit_pos;\n int<lower=0> kit_neg;\n int<lower=0> n_kit_neg;\n}\n\nparameters {\n real<lower=0, upper=1> theta; //seroprevalencia\n real<lower=0, upper=1> sens; //sensibilidad\n real<lower=0, upper=1> esp; //especificidad\n}\n\ntransformed parameters {\n real<lower=0, upper=1> prob_pos;\n\n prob_pos = theta * sens + (1 - theta) * (1 - esp);\n\n}\nmodel {\n // modelo de número de positivos\n n ~ binomial(N, prob_pos);\n // modelos para resultados del kit\n kit_pos ~ binomial(n_kit_pos, sens);\n kit_neg ~ binomial(n_kit_neg, esp);\n // iniciales para cantidades no medidas\n theta ~ beta(1.0, 10.0);\n sens ~ beta(2.0, 1.0);\n esp ~ beta(2.0, 1.0);\n}\n\n\nEn este caso, no pondremos información acerca de positivos en la prueba:\n\ndatos_lista <- list(N = 0, n = 0,\n kit_pos = 103, n_kit_pos = 122,\n kit_neg = 399, n_kit_neg = 401)\najuste <- mod_sc$sample(data = datos_lista, refresh = 1000, iter_sampling = 400)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 Iteration: 1 / 1400 [ 0%] (Warmup) \nChain 1 Iteration: 1000 / 1400 [ 71%] (Warmup) \nChain 1 Iteration: 1001 / 1400 [ 71%] (Sampling) \nChain 1 Iteration: 1400 / 1400 [100%] (Sampling) \nChain 1 finished in 0.0 seconds.\nChain 2 Iteration: 1 / 1400 [ 0%] (Warmup) \nChain 2 Iteration: 1000 / 1400 [ 71%] (Warmup) \nChain 2 Iteration: 1001 / 1400 [ 71%] (Sampling) \nChain 2 Iteration: 1400 / 1400 [100%] (Sampling) \nChain 2 finished in 0.0 seconds.\nChain 3 Iteration: 1 / 1400 [ 0%] (Warmup) \nChain 3 Iteration: 1000 / 1400 [ 71%] (Warmup) \nChain 3 Iteration: 1001 / 1400 [ 71%] (Sampling) \nChain 3 Iteration: 1400 / 1400 [100%] (Sampling) \nChain 3 finished in 0.0 seconds.\nChain 4 Iteration: 1 / 1400 [ 0%] (Warmup) \nChain 4 Iteration: 1000 / 1400 [ 71%] (Warmup) \nChain 4 Iteration: 1001 / 1400 [ 71%] (Sampling) \nChain 4 Iteration: 1400 / 1400 [100%] (Sampling) \nChain 4 finished in 0.0 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.0 seconds.\nTotal execution time: 0.5 seconds.\n\nsims <- ajuste$draws(c(\"theta\", \"sens\", \"esp\"), format = \"df\")\nresumen <- ajuste$summary(c(\"theta\"))\n\n\nggplot(sims, aes(x = theta, y = sens)) + geom_point() +\n scale_x_sqrt()\n\n\n\n\n\n\n\n\nNo vemos ninguna asocación entre estas dos variables.\nSin embargo, al condicionar al valor de Positivos, creamos una relación que no podemos interpretar como casual. En este caso particular supondremos prácticamente fija la sensibilidad para ver solamente lo que sucede en el colisionador de especificidad y número de positivos (la especificidad en este ejemplo es más crítica):\n\ndatos_lista <- list(N = 3300, n = 50,\n kit_pos = 1030000, n_kit_pos = 1220000, # números grandes para que esté practicamente\n# fija la sensibilidad\n kit_neg = 399, n_kit_neg = 401)\najuste <- mod_sc$sample(data = datos_lista, refresh = 1000, iter_sampling = 400)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 Iteration: 1 / 1400 [ 0%] (Warmup) \nChain 1 Iteration: 1000 / 1400 [ 71%] (Warmup) \nChain 1 Iteration: 1001 / 1400 [ 71%] (Sampling) \nChain 1 Iteration: 1400 / 1400 [100%] (Sampling) \nChain 1 finished in 0.0 seconds.\nChain 2 Iteration: 1 / 1400 [ 0%] (Warmup) \nChain 2 Iteration: 1000 / 1400 [ 71%] (Warmup) \nChain 2 Iteration: 1001 / 1400 [ 71%] (Sampling) \nChain 2 Iteration: 1400 / 1400 [100%] (Sampling) \nChain 2 finished in 0.0 seconds.\nChain 3 Iteration: 1 / 1400 [ 0%] (Warmup) \nChain 3 Iteration: 1000 / 1400 [ 71%] (Warmup) \nChain 3 Iteration: 1001 / 1400 [ 71%] (Sampling) \nChain 3 Iteration: 1400 / 1400 [100%] (Sampling) \nChain 3 finished in 0.0 seconds.\nChain 4 Iteration: 1 / 1400 [ 0%] (Warmup) \nChain 4 Iteration: 1000 / 1400 [ 71%] (Warmup) \nChain 4 Iteration: 1001 / 1400 [ 71%] (Sampling) \nChain 4 Iteration: 1400 / 1400 [100%] (Sampling) \nChain 4 finished in 0.0 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.0 seconds.\nTotal execution time: 0.5 seconds.\n\nsims <- ajuste$draws(c(\"theta\", \"sens\", \"esp\"), format = \"df\")\nresumen <- ajuste$summary(c(\"theta\"))\n\n\nggplot(sims, aes(x = theta, y = esp)) + geom_point() \n\n\n\n\n\n\n\n\nY vemos que condiconando al colisionador, obtenemos una relación fuerte entre prevalencia y especificidad de la prueba: necesitaríamos más datos de especificidad para obtener una estimación útil.\n\nLa razón de que la especificidad es más importante en este ejemplo es que la prevalencia es muy baja al momento del estudio, y los falsos positivos pueden introducir más error en la estimación\nTambién repetimos nótese que el análisis correcto de estos datos no se puede hacer con intervalos separados para cada cantidad, sino que debe examinarse la conjunta de estos parámetros.\n\n\nCon estas tres estructuras elementales podemos entender de manera abstracta la existencia o no de asociaciones entre nodos de cualquier gráfica dirigida.",
+ "text": "5.8 Razonamiento de descendientes\nCondicionar a un descendiente puede entenderse como “condicionar parcialmente” o “débilmente” a los padres de ese descendiente.\nPor ejemplo, condicionar a un colisionador también produce dependencias condicionales:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n X\n Y\n Z\n A\n edge [minlen = 3]\n X -> Z\n Y -> Z\n Z -> A\n}\n\", width = 200, height = 50)\n\n\n\n\n\n\nEn este caso,\n\n\\(X\\) y \\(Y\\) son independientes\n\\(X\\) y \\(Y\\) son dependientes si condicionamos a \\(A\\).\n\nDependiendo de la naturaleza de la asociación entre el colisionador \\(Z\\) y su descendiente \\(A\\), esta dependencia puede ser más fuerte o más débil.\nPor ejemplo, en nuestro ejemplo donde el pasto mojado es un colisionador entre cuánta agua dieron los aspersores y cuánta lluvia cayó, un descendiente del pasto mojado es el estado de las plantas del jardín. Aunque los aspersores trabajan independientemente de la lluvia, si observamos que las plantas se secaron entonces lluvia y aspersores están correlacionados: por ejemplo, si noto que los aspersores están descompuestos, entonces concluimos que no hubo lluvia.\n\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n X [label = lluvia]\n Y [label = aspersores]\n Z [label = humedad]\n A [label = plantas]\n edge [minlen = 3]\n X -> Z\n Y -> Z\n Z -> A\n}\n\", width = 200, height = 50)\n\n\n\n\n\n\nEjemplo\n\nsimular_desc <- function(n = 10){\n x <- rbern(n, 0.5) \n y <- rbinom(n, 2, 0.7)\n z <- rbern(n, p = 0.1 + 0.7 * x * (y > 1)) \n a <- rbern(n, p = 0.5 + 0.5 * z)\n tibble(x, z, y, a)\n}\nsims_colisionador <- simular_desc(50000)\n# No hay correlación\ncor(sims_colisionador$x, sims_colisionador$y)\n\n[1] 0.001412209\n\n\nSin embargo,\n\ncor(sims_colisionador |> filter(a ==0) |> select(x,y))\n\n x y\nx 1.0000000 -0.2798845\ny -0.2798845 1.0000000\n\n\n\ncor(sims_colisionador |> filter(a ==1) |> select(x,y))\n\n x y\nx 1.0000000 0.1127725\ny 0.1127725 1.0000000\n\n\n\n\n5.8.1 Ejemplo: dependencias de colisionador\nVerificamos que en nuestro modelo de Santa Clara, efectivamente nuestro modelo no implica ninguna dependencia no condicional entre sensibilidad de la prueba y prevalencia. Eso debería ser claro de la simulación, pero de todas formas lo checamos\n\nlibrary(cmdstanr)\nmod_sc <- cmdstan_model(\"./src/sclara.stan\")\nprint(mod_sc)\n\ndata {\n int<lower=0> N;\n int<lower=0> n;\n int<lower=0> kit_pos;\n int<lower=0> n_kit_pos;\n int<lower=0> kit_neg;\n int<lower=0> n_kit_neg;\n}\n\nparameters {\n real<lower=0, upper=1> theta; //seroprevalencia\n real<lower=0, upper=1> sens; //sensibilidad\n real<lower=0, upper=1> esp; //especificidad\n}\n\ntransformed parameters {\n real<lower=0, upper=1> prob_pos;\n\n prob_pos = theta * sens + (1 - theta) * (1 - esp);\n\n}\nmodel {\n // modelo de número de positivos\n n ~ binomial(N, prob_pos);\n // modelos para resultados del kit\n kit_pos ~ binomial(n_kit_pos, sens);\n kit_neg ~ binomial(n_kit_neg, esp);\n // iniciales para cantidades no medidas\n theta ~ beta(1.0, 10.0);\n sens ~ beta(2.0, 1.0);\n esp ~ beta(2.0, 1.0);\n}\n\n\nEn este caso, no pondremos información acerca de positivos en la prueba:\n\ndatos_lista <- list(N = 0, n = 0,\n kit_pos = 103, n_kit_pos = 122,\n kit_neg = 399, n_kit_neg = 401)\najuste <- mod_sc$sample(data = datos_lista, refresh = 1000, iter_sampling = 400)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 Iteration: 1 / 1400 [ 0%] (Warmup) \nChain 1 Iteration: 1000 / 1400 [ 71%] (Warmup) \nChain 1 Iteration: 1001 / 1400 [ 71%] (Sampling) \nChain 1 Iteration: 1400 / 1400 [100%] (Sampling) \nChain 1 finished in 0.0 seconds.\nChain 2 Iteration: 1 / 1400 [ 0%] (Warmup) \nChain 2 Iteration: 1000 / 1400 [ 71%] (Warmup) \nChain 2 Iteration: 1001 / 1400 [ 71%] (Sampling) \nChain 2 Iteration: 1400 / 1400 [100%] (Sampling) \nChain 2 finished in 0.0 seconds.\nChain 3 Iteration: 1 / 1400 [ 0%] (Warmup) \nChain 3 Iteration: 1000 / 1400 [ 71%] (Warmup) \nChain 3 Iteration: 1001 / 1400 [ 71%] (Sampling) \nChain 3 Iteration: 1400 / 1400 [100%] (Sampling) \nChain 3 finished in 0.0 seconds.\nChain 4 Iteration: 1 / 1400 [ 0%] (Warmup) \nChain 4 Iteration: 1000 / 1400 [ 71%] (Warmup) \nChain 4 Iteration: 1001 / 1400 [ 71%] (Sampling) \nChain 4 Iteration: 1400 / 1400 [100%] (Sampling) \nChain 4 finished in 0.0 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.0 seconds.\nTotal execution time: 0.5 seconds.\n\nsims <- ajuste$draws(c(\"theta\", \"sens\", \"esp\"), format = \"df\")\nresumen <- ajuste$summary(c(\"theta\"))\n\n\nggplot(sims, aes(x = theta, y = sens)) + geom_point() +\n scale_x_sqrt()\n\n\n\n\n\n\n\n\nNo vemos ninguna asocación entre estas dos variables.\nSin embargo, al condicionar al valor de Positivos, creamos una relación que no podemos interpretar como casual. En este caso particular supondremos prácticamente fija la sensibilidad para ver solamente lo que sucede en el colisionador de especificidad y número de positivos (la especificidad en este ejemplo es más crítica):\n\ndatos_lista <- list(N = 3300, n = 50,\n kit_pos = 1030000, n_kit_pos = 1220000, # números grandes para que esté practicamente\n# fija la sensibilidad\n kit_neg = 399, n_kit_neg = 401)\najuste <- mod_sc$sample(data = datos_lista, refresh = 1000, iter_sampling = 400)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 Iteration: 1 / 1400 [ 0%] (Warmup) \nChain 1 Iteration: 1000 / 1400 [ 71%] (Warmup) \nChain 1 Iteration: 1001 / 1400 [ 71%] (Sampling) \nChain 1 Iteration: 1400 / 1400 [100%] (Sampling) \nChain 1 finished in 0.0 seconds.\nChain 2 Iteration: 1 / 1400 [ 0%] (Warmup) \nChain 2 Iteration: 1000 / 1400 [ 71%] (Warmup) \nChain 2 Iteration: 1001 / 1400 [ 71%] (Sampling) \nChain 2 Iteration: 1400 / 1400 [100%] (Sampling) \nChain 2 finished in 0.0 seconds.\nChain 3 Iteration: 1 / 1400 [ 0%] (Warmup) \nChain 3 Iteration: 1000 / 1400 [ 71%] (Warmup) \nChain 3 Iteration: 1001 / 1400 [ 71%] (Sampling) \nChain 3 Iteration: 1400 / 1400 [100%] (Sampling) \nChain 3 finished in 0.0 seconds.\nChain 4 Iteration: 1 / 1400 [ 0%] (Warmup) \nChain 4 Iteration: 1000 / 1400 [ 71%] (Warmup) \nChain 4 Iteration: 1001 / 1400 [ 71%] (Sampling) \nChain 4 Iteration: 1400 / 1400 [100%] (Sampling) \nChain 4 finished in 0.0 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 0.0 seconds.\nTotal execution time: 0.5 seconds.\n\nsims <- ajuste$draws(c(\"theta\", \"sens\", \"esp\"), format = \"df\")\nresumen <- ajuste$summary(c(\"theta\"))\n\n\nggplot(sims, aes(x = theta, y = esp)) + geom_point() \n\n\n\n\n\n\n\n\nY vemos que condiconando al colisionador, obtenemos una relación fuerte entre prevalencia y especificidad de la prueba: necesitaríamos más datos de especificidad para obtener una estimación útil.\n\nLa razón de que la especificidad es más importante en este ejemplo es que la prevalencia es muy baja al momento del estudio, y los falsos positivos pueden introducir más error en la estimación\nTambién repetimos nótese que el análisis correcto de estos datos no se puede hacer con intervalos separados para cada cantidad, sino que debe examinarse la conjunta de estos parámetros.\n\n\nCon estas tres estructuras elementales podemos entender de manera abstracta la existencia o no de asociaciones entre nodos de cualquier gráfica dirigida.",
"crumbs": [
"5Modelos gráficos y causalidad"
]
@@ -388,5 +388,85 @@
"crumbs": [
"5Modelos gráficos y causalidad"
]
+ },
+ {
+ "objectID": "06-calculo-do.html",
+ "href": "06-calculo-do.html",
+ "title": "6 Identificación y cálculo-do",
+ "section": "",
+ "text": "6.1 Cambiando el proceso generador de datos\nComenzamos con el ejemplo más simple de una variable confusora:\ngrViz(\"\n digraph {\n node [shape = plaintext];\n X [label = 'X'];\n Y [label = 'Y'];\n U [label = 'U'];\n X -> Y;\n U-> X ;\n U -> Y;\n {rank = same; X; Y;}\n }\n \")\nNos interesa estimar el efecto causal de \\(X\\) sobre \\(Y\\). Sucede que en muchas ocasiones existen variables como \\(U\\) que son causas comunes de \\(X\\) y \\(Y\\). Como vimos, esto implica que no podemos simplemente ver la correlación entre \\(X\\) y \\(Y\\) para entender el efecto de \\(X\\) sobre \\(Y\\), pues una causa común de variación conjunta entre estas dos variables. Esta variable \\(U\\) puede ser observada o no.\nEste tipo de confusores ocurren muchas veces en datos observacionales (es decir, de un proceso o sistema que funcione sin intervención de los investigadores). Por ejemplo, si un estudio observa que aquellos que se aplicaron (voluntariamente) un tratamiento \\(X\\), tienen menor riesgo de hospitalización \\(Y\\) por cierta enfermadad. Sin embargo, se observa también que aquellos que se aplicaron el tratamiento tienen menos riesgo de tener accidentes viales. Esto indica que la observación de la reducción de riesgo de hospitalización entre los que escogieron el tratamiento probablemente se debe al menos en parte a una variable confusora (por ejemplo, qué tipo de actividades hacen, qué tan cautelosos son, etc.)",
+ "crumbs": [
+ "6Identificación y cálculo-do"
+ ]
+ },
+ {
+ "objectID": "06-calculo-do.html#cambiando-el-proceso-generador-de-datos",
+ "href": "06-calculo-do.html#cambiando-el-proceso-generador-de-datos",
+ "title": "6 Identificación y cálculo-do",
+ "section": "",
+ "text": "6.1.1 Experimentación\nCuando es posible, podemos proponer generar nuevos datos donde alteramos el proceso generador. Una forma muy efectiva y útil, que es muy conveniente cuando es posible, es controlar la asignación del tratamiento. Si en el diagrama anterior, diseñamos un estudio donde observamos a un grupo de personas para las cuales el tratamiento se asignó de acuerdo a un proceso aleatorio, entonces el nuevo diagrama para este nuevo proceso generador es:\n\ngrViz(\"\n digraph {\n node [shape = plaintext];\n X [label = 'X'];\n Y [label = 'Y'];\n R\n U [label = 'U'];\n R -> X\n X -> Y;\n U -> Y;\n {rank = same; X; Y;}\n }\n \")\n\n\n\n\n\nNótese que:\n\nLa variable \\(R\\) no puede ser endógena (es decir, ninguna flecha del sistema puede incidir en ella), pues se utiliza un dado o algo totalmente no relacionado para asignar el tratamiento. Por ejemplo, también podríamos asignar el tratamiento utilizando la segunda letra del apellido de las personas.\nNo puede existir una flecha de \\(U\\) a \\(X\\), pues nada en \\(X\\) responde a cambios en \\(X\\), qué solo depende del proceso de aleatorización \\(R\\).\n\nEn este caso, no es necesario estratificar por ninguna variable y podemos proponer directamente un modelo estadístico para \\(Y\\) en función de \\(X\\) que nos permita estimar el efecto causal de \\(X\\) sobre \\(Y\\).\n\n\n\n\n\n\nExperimentos\n\n\n\nEsto describe la idea básica de un experimento simple: es una herramienta para modificar el proceso generador de datos que nos permite identificar efectos causales de manera relativamente simple.\nCuando es posible hacer experimentos de calidad, esta puede ser la mejor forma de estimar efectos causales.\n\n\nEn muchos casos, sin embargo, no es posible hacer experimentos de calidad. Hay varias diversas razones, por ejemplo cuando se trata de experimentos que involucran personas:\n\nNo es ético aleatorizar: es totalmente inaceptable asignar aleatoriamente a personas a un tratamientos como fumar 20 cigarros al día, o aleatorizar a niños a recibir educación o no.\nAleatorización imposible o imperfecta: no es posible lograr un control total sobre la asignación del tratamiento, y la adherencia al tratamiento asignado de las personas puede variar (por ejemplo, uso de tapabocas en escuelas). A lo más podemos considerar los efectos de una política que intenta tratar a una selección aleatoria de individuos (IIT, o intent-to-treat).\n\nAsí que muchas preguntas causales no están sujetas a modificaciones del proceso generador de datos mediante aleatorización, y es necesario recurrir a otras estrategias.",
+ "crumbs": [
+ "6Identificación y cálculo-do"
+ ]
+ },
+ {
+ "objectID": "06-calculo-do.html#el-operador-do",
+ "href": "06-calculo-do.html#el-operador-do",
+ "title": "6 Identificación y cálculo-do",
+ "section": "6.2 El operador do",
+ "text": "6.2 El operador do\nRegresamos al diagrama original donde \\(U\\) es una causa común de \\(X\\) y \\(Y\\), y que no tenemos recursos o no es posible hacer un experimento. ¿Existe algún procedimiento estadístico que nos permita estimar el efecto causal de \\(X\\) sobre \\(Y\\)?\nEscribiremos la distribución condicional de la respuesta \\(Y\\) dada una manipulación de \\(X\\) como sigue (es decir, en la situación experimental):\n\\[p(Y| do(X=x))\\]\nEsto significa: ¿cómo se distribuye la \\(Y\\) dado que intervenimos en la población completa (aunque podemos también considerar subpoblaciones más adelante) para poner en \\(X=x\\)? En primer lugar, notemos que esto no es lo mismo que la distribución condicional usual\n\\[p(Y|X=x),\\] que siempre podemos estimar directamente de los datos, y no es la que nos interesa. En el siguiente ejemplo vemos la distinción entre las dos distribuciones:\n\nEjemplo (Pearl)\nSupongamos que tenemos el siguiente modelo del diagrama causal:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n\n edge [minlen = 3]\n T -> A\n T -> Z\n \n \n}\n\")\n\n\n\n\n\n\ndonde \\(T\\) es la temperatura, \\(A\\) son las unidades de agua embotellada vendidas y \\(Z\\) es la actividad de los mosquitos (medido con muestreo, por ejemplo).\nNo interesa contestar la pregunta: ¿qué tanto influyen las ventas de agua embotellada en la actividad de los mosquitos? Del diagrama, sabemos que no hay ningún camino causal de \\(Z\\) a \\(A\\), por lo que nuestra respuesta debería ser igual a 0.\nSin embargo, sabemos que estas dos variables están asociadas (por el análisis de DAGs), de manera que describir cómo cambia \\(p(Z|A)\\) cuando condicionamos a distintos valores de \\(A\\) no responde nuestra pregunta. La distribución \\(p(Z|do(A = a))\\) nos dice cómo se distribuye \\(Z\\) cuando manipulamos \\(a\\) artificialmente. Por ejemplo, si cerramos todas las tiendas un día haciendo \\(do(A=0)\\), veríamos que esta variable no tiene efecto sobre la actividad de mosquitos, por ejemplo comparado con \\(do(A = 10000)\\).\nIlustramos la diferencia entre \\(p(Y|X)\\) y \\(p(Y|do(X))\\) simulando del ejemplo anterior. Supondremos que sólo consideramos un día del año a lo largo de varios años, para no modelar el comportamiento cíclo de la temperatura:\n\nsimular_t <- function(n = 10, dia = 150){\n # simular un año, alrededor del día 160 (en junio)\n t_maxima <- rnorm(n, 28, 2)\n mosquitos <- rpois(n, 250 + 10 * (t_maxima - 28))\n a_unidades <- rnorm(n, 20000 + 2000 * (t_maxima - 28), 2000)\n tibble(t_maxima, a_unidades, mosquitos)\n}\nset.seed(128)\nsimular_dias <- simular_t(50)\n\nSi simulamos, vemos que \\(mosquitos\\) y \\(unidades\\) son dependientes, pues tenemos un camino abierto dado por la bifurcación en temperatura:\n\nggplot(simular_dias, aes(x = a_unidades, y = mosquitos)) + geom_point() +\n geom_smooth(method = \"loess\", method.args = list(degree = 1)) +\n xlab(\"Ventas de agua embotellada\")\n\n`geom_smooth()` using formula = 'y ~ x'\n\n\n\n\n\n\n\n\n\nSabemos que esta asociación no es causal, pues no hay caminos causales entre estas variables dos variables, pero que hay una dependencia debido a la bifurcación en \\(T\\). La gráfica muestra que la media condicional \\(E[M|A=a]\\) depende fuertemente de \\(a\\), lo que quiere decir que \\(p(m|a)\\) depende de \\(a\\) fuertemente.\n\n\nUna intervención simple\nEn este caso, nos interesaría saber qué sucede si alteramos artificalmente el número de botellas de agua vendidas (puedes imaginar distintas maneras de hacer esto).\nCuando hacemos esto, quitamos las aristas que van hacia \\(A\\), pues \\(A\\) ya no está determinado por el proceso generador de datos. Tenemos entonces la nueva gráfica:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n A\n edge [minlen = 3]\n U_t -> T\n T -> Z\n U_m -> Z\n{ rank = same; A; Z }\n}\n\")\n\n\n\n\n\n\nEn esta nueva gráfica, \\(A\\) y \\(Z\\) son independientes, que es la respuesta correcta. Como cambiamos la gráfica, su proceso generador es diferente al original de los datos observados. Sin embargo, en este ejemplo puedes ver por qué es claro que el cambio que hicimos (manipular \\(A\\) en lugar de que esté determinado por su proceso generador original) no cambia el modelo de \\(Z\\), de manera que podemos simular de nuestro nuevo proceso generador donde manipulamos \\(A\\):\n\nsimular_cirugia <- function(n = 10, a_unidades = a_unidades){\n # simular un año, alrededor del día 160 (en junio)\n t_maxima <- rnorm(n, 28, 2)\n #### cirugía #########\n # ahora a_unidades es fijado por nosotros:\n # a_unidades <- rnorm(n, 20000 + 2000 * (t_maxima - 28), 2000)\n a_unidades <- a_unidades\n ######################\n mosquitos <- rpois(n, 250 + 10 * (t_maxima - 28))\n tibble(t_maxima, a_unidades, mosquitos)\n}\n\nY ahora simulamos y graficamos \\(p(Z|do(A=a))\\) para distintos valores de \\(a\\):\n\nset.seed(128)\nsimular_dias_2 <- map_df(seq(10000, 30000, 1000),\n \\(u) simular_cirugia(50, a_unidades = u))\n\n\nggplot(simular_dias_2, aes(x = a_unidades, y = mosquitos)) +\n geom_point() + geom_smooth()\n\n`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = \"cs\")'\n\n\n\n\n\n\n\n\n\ny vemos, como esperaríamos, que no hay relación entre unidades de agua embotellada y mosquitos.",
+ "crumbs": [
+ "6Identificación y cálculo-do"
+ ]
+ },
+ {
+ "objectID": "06-calculo-do.html#cálculo-do-de-pearl",
+ "href": "06-calculo-do.html#cálculo-do-de-pearl",
+ "title": "6 Identificación y cálculo-do",
+ "section": "6.3 Cálculo-do de Pearl",
+ "text": "6.3 Cálculo-do de Pearl\nEl cálculo do nos da reglas para operar con probabilidades que incluyen nuestro operador do de intervención. En este ejemplo particular, veremos cómo es el argumento:\nNótese que al intervenir \\(A\\) hemos modificado el proceso generador. Si la conjunta original tiene distribución \\(p\\), escribimos \\(p_m\\) para la conjunta de la gráfica modificada, de manera que \\(p(Z|do(A)) = p_m(Z|A)\\): con esto podemos pasar de una pregunta causal (lado izquierdo con operador do) a una estadpística (lado derecho).\nAunque intuitivamente vimos cómo simular de esta distribución arriba, especificamos abajo qué reglas son las que nos permiten hacer esto: ¿cómo calculamos \\(p_m\\)?\nEn primer lugar, consideremos la marginal \\(p_m(T)\\). Esta marginal es invariante a nuestra cirugía, pues la arista \\(T\\to A\\) que eliminamos \\(T\\) no afecta el proceso que determina \\(T\\). De modo que la marginal del proceso modificado es igual a la marginal observada:\n\\[p_m(T) = p(T)\\] En segundo lugar, tenemos que\n\\[p_m(Z|T=t,A=a) = p(Z|T=t,A=a),\\] Pues el proceso por el cual \\(Z\\) responde a \\(T\\) y \\(A\\) es el mismo, no importa si \\(A\\) fue modificada artificalmente o no.\nJuntamos estos argumentos. Primero, por definición,\n\\[p(Z|do(A=a)) = p_m(Z|A=a)\\]\nPor la regla de probabilidad total, podemos condicionar todo a \\(T\\) y marginalizar. La segunda igualdad la obtenemos por la independencia entre \\(T\\) y \\(Z\\) en nuestra gráfica modificada (están \\(d\\) separadas):\n\\[p_m(z|a) = \\int p_m(z|a,t)p_m(t|a)dt = \\int p_m(z|a,t)p_m(t)dt\\] En segunda igualdad, nótese que cambiamos \\(p_m(t|a) = p_m(t)\\), lo cual podemos verificar pues en la gráfica modificada \\(A\\) y \\(T\\) están \\(d\\)-separados, lo que implica que son condicionalmente independientes.\nFinalmente, las últimas dos distribuciones podemos extraerlas de los datos, como explicamos arriba \\(p_m(z|t,a) = p(z|t,a)\\) y \\(p_m(t) = p(t),\\) y terminamos con la fórmula:\n\\[p(z|do(a))=p_m(z|a) = \\int p(z|a,t)p(t)dt \\]\nLas dos distribuciones de la derecha están en el contexto de \\(p\\), el proceso generador de datos original. Así que podemos estimarlas de los datos observados.\n\nEste argumento justifica el proceso que hicimos arriba: simulamos primero \\(T\\) con su proceso generador, y después simulamos \\(Z\\) condicional a \\(A\\) y \\(T\\) según el proceso generador original, el cual no depende de \\(A\\) en este ejemplo.\n\nEn el caso de arriba, simulamos de la distribución para entender cómo se distribuía \\(Z\\) dependiendo de modificaciones a \\(A\\). Muchas veces nos interesa calcular solamente la esperanza condicional, es decir, cuál es el valor esperado de la variable de interés dado el nivel intervenido, es decir:\n\\(E(Z|do(A=a)) = E_m(Z|A =a),\\)\nque mostramos arriba con la línea ajustada. También quisiéramos calcular contrastes particulares, como qué pasaría si las ventas de agua las aumentamos en 10 mil unidades:\n\\[E(Z|do(A=30000)) - E(Z|do(A=20000)),\\] que podemos calcular de manera simple con simulación:\n\nsimular_contraste <- map_df(c(20000, 30000),\n \\(u) simular_cirugia(1000, a_unidades = u)) |> \n group_by(a_unidades) |> \n summarise(media_mosquitos = mean(mosquitos))\nsimular_contraste\n\n# A tibble: 2 × 2\n a_unidades media_mosquitos\n <dbl> <dbl>\n1 20000 250.\n2 30000 249.\n\n\nY vemos que no hay diferencia entre las dos medias.\n\nEjemplo\nAhora hagamos otro ejemplo donde hay una relación causal que queremos estimar. Imaginemos una ciudad en donde temperaturas altas producen desabasto de agua en algunos hogares, debido a un aumento del riego y uso de agua en general. Nos interesa estimar el efecto del desabasto en las compras de agua embotellada. Nuestro diagrama ahora es:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n\n edge [minlen = 3]\n U_t -> T\n T -> A\n T -> D\n D -> A\n U_a -> A\n U_d -> D\n\n{ rank = same; A; D }\n\n}\n\")\n\n\n\n\n\n\n\nsimular_t <- function(n = 10, dia = 150){\n # simular un año, alrededor del día 160 (en junio)\n t_maxima <- rnorm(n, 28, 2)\n u <- rnorm(n, 0, 1)\n desabasto_agua <- 1/(1 + exp(-(t_maxima - 28) + u))\n unidades <- rnorm(n, 20000 + 2000 * (t_maxima - 28) + 8000*desabasto_agua, 2000)\n tibble(t_maxima, unidades, desabasto_agua)\n}\nset.seed(128)\nsimular_dias <- simular_t(150)\n\n\nggplot(simular_dias, aes(x = desabasto_agua, y = unidades)) + \n geom_point() + geom_smooth()\n\n`geom_smooth()` using method = 'loess' and formula = 'y ~ x'\n\n\n\n\n\n\n\n\n\nLa correlación parece muy fuerte, sin embargo, sabemos que hay un camino no causal de asociación entre estas dos variables.\nIgual que en ejemplo anterior, vamos a intervenir teóricamente en el desabasto de agua. Después de la cirugía, nuestro diagrama modificado es:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n\n edge [minlen = 3]\n U_t -> T\n T -> A\n D -> A\n U_a -> A\n{ rank = same; A; D }\n\n}\n\")\n\n\n\n\n\n\nAhora queremos calcular \\(p(a|do(d)) = p_m(a|d)\\) en función de los datos. Siguiendo el mismo argumento que en el ejemplo anterior, sabemos que tenemos que estratificar o condicionar a \\(T\\) para poder usar nuestro proceso generador de observaciones, y obtenemos:\n\\[p(a|do(d))=p_m(a|d) = \\int p(a|d,t)p(t)dt \\] Aunque a veces es posible calcular analíticamente el lado derecho analíticamente, podemos simular como hicimos en los ejemplos anteriores:\n\nsimular_cirugia <- function(n = 10, da = 0){\n # simular un año, alrededor del día 160 (en junio)\n t_maxima <- rnorm(n, 28, 2)\n ### cirugía ####\n #u <- rnorm(n, 0, 1) \n desabasto_agua <- da\n ######\n unidades <- rnorm(n, 20000 + 2000 * (t_maxima - 28) + 8000*desabasto_agua, 2000)\n tibble(t_maxima, unidades, desabasto_agua)\n}\nset.seed(128)\nsimular_dias_c <- map_df(seq(0, 1, 0.1), \\(da) simular_cirugia(1000, da = da))\n\n\nggplot(simular_dias_c, aes(x = desabasto_agua, y = unidades)) + \n geom_point() + geom_smooth()\n\n`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = \"cs\")'\n\n\n\n\n\n\n\n\n\nPodemos también resumir promediando:\n\nefecto_verdadero_desabasto <- simular_dias_c |> \n group_by(desabasto_agua) |> \n summarise(media_unidades = mean(unidades)) |> \n rename(desabasto = desabasto_agua)\nggplot(efecto_verdadero_desabasto,\n aes(x = desabasto, y = media_unidades)) + \n geom_point() + geom_smooth()\n\n`geom_smooth()` using method = 'loess' and formula = 'y ~ x'\n\n\n\n\n\n\n\n\n\nY este es el efecto causal del desabasto de agua. No tenemos medidas de incertidumbre pues conocemos todos los parámetros de los modelos. La media condicional parece ser lineal, así que podríamos resumir con un modelo lineal:\n\n# Modelo 1 (con datos de intervención)\nlm(unidades ~ desabasto_agua, simular_dias_c)\n\n\nCall:\nlm(formula = unidades ~ desabasto_agua, data = simular_dias_c)\n\nCoefficients:\n (Intercept) desabasto_agua \n 19831 8272 \n\n\nAproximadamente, cada incremento en puntos porcentuales de 10% en desabasto incrementa las ventas en unas 800 unidades. Compara con el análisis donde no estratificamos o controlamos por la temperatura:\n\n# Modelo 2\nlm(unidades ~ desabasto_agua, simular_dias)\n\n\nCall:\nlm(formula = unidades ~ desabasto_agua, data = simular_dias)\n\nCoefficients:\n (Intercept) desabasto_agua \n 14102 19491 \n\n\nOtra forma de estratificar es ajustando un modelo que incluye la variable de temperatura. Podríamos hacer\n\n# Modelo 3\nlm(unidades ~ desabasto_agua + t_maxima, simular_dias)\n\n\nCall:\nlm(formula = unidades ~ desabasto_agua + t_maxima, data = simular_dias)\n\nCoefficients:\n (Intercept) desabasto_agua t_maxima \n -35030 8648 1948",
+ "crumbs": [
+ "6Identificación y cálculo-do"
+ ]
+ },
+ {
+ "objectID": "06-calculo-do.html#fórmula-de-ajuste",
+ "href": "06-calculo-do.html#fórmula-de-ajuste",
+ "title": "6 Identificación y cálculo-do",
+ "section": "6.4 Fórmula de ajuste",
+ "text": "6.4 Fórmula de ajuste\nEn resumen, tenemos la primera regla de Pearl de inferencia causal:\n\n\n\n\n\n\nFórmula de ajuste (Pearl)\n\n\n\nSea \\(G\\) donde los padres de \\(X\\) son \\(Z_1,Z_2\\). El efecto causal total de \\(X\\) en \\(Y\\) se puede calcular como\n\\[p(y|do(x)) = \\int p(y|x, z_1,z_2) p(z_1,z_2)\\, dz_1dz_2\\] Es decir, condicionamos al valor de \\(x\\) y todos los padres de \\(X\\) para calcular \\(p(y|x,z_1,z_2)\\), y después marginalizamos sobre los padres.\n\n\nEsta fórmula se extiende a más de dos padres \\(Z_1,Z_2,Z_3,\\ldots, Z_k\\).\n\n\n\n\n\n\nTip\n\n\n\nA este proceso se llama de diferentes maneras en distintos contextos:\n\nEstamos calculando el efecto causal estratificando por las variables \\(z\\).\nControlamos por las variables \\(z\\) para calcular el efecto causal.\n\n\n\nPodemos pensar en esta fórmula de dos maneras: en primer lugar, si estamos modelando toda nuestra gráfica causal, podemos simular de la conjunta de la gráfica mutilada:\n\nFijando el nivel del tratamiento \\(T\\)\nSimulando \\(p(z_1,z_2,\\ldots, z_k)\\) de nuestro modelo completo (y tomar sólo los valores de las \\(z\\)’s).\nUsar \\(t\\) y las \\(z\\) simuladas para simular \\(y\\).\nAl final, nótese que nos quedan simulaciones de \\(p_m(y|t)\\) (marginalizamos sobre las \\(z\\)).\n\nEl otro enfoque busca sólo construir modelos para la parte que nos interesa:\n\nConstruir un modelo separado para \\(p(z_1, z_2,\\ldots, z_k) = p(z)\\) (que puede ser difícil si tenemos muchas variables) a partir los datos. Podemos también simular tomando al azar esta variables de nuestros datos.\nConstruir un modelo \\(p(y|t, z)\\) para simular la \\(y\\) a partir de los datos.\nMarginalizar sobre las \\(z\\)’s para quedarnos con \\(p_m(y|t)\\)\n\nFinalmente, si tenemos un modelo \\(p(y| t, z)\\) podemos también investigar cómo se comporta \\(E[y|t_2,z] - E[y|t_1,z]\\) para distintos combinaciones de valores de \\(Z\\).\nNota 1: Con este principio podemos resolver algunos problemas, pero no todos. Veremos que en algunos casos existen padres que no son observados, por ejemplo, no es posible condicionar para usar la fórmula de ajuste y es necesario desarrollar otras estrategias.\nNota 2: En regresión lineal, cuando incluímos una variable en el modelo (que consideramos una variable control), estamos estratificando por ella: por ejemplo, en el modelo lineal \\(U\\sim N(m_u(d,t), \\sigma_u)\\), donde\n\\[m_u = \\beta_0 +\\beta_1 d + \\beta_2 t\\] Estamos calculando un estimador para cada valor de \\(T=t\\), que es:\n\\[m_u = (\\beta_0 + \\beta_2 t) + \\beta_1 d = \\gamma_0 + \\gamma_1 d\\] Esta es una de las maneras más simples de obtener el efecto de \\(d\\) estratificando por, o controlando por \\(t\\), siempre y cuando los modelos lineales sean apropiados.\nNótese que en este último caso, tenemos que el efecto de \\(d\\) no depende de las covariables, de forma que no es necesario hacer el promedio sobre la conjunta, es decir, suponemos que el efecto causal es el mismo independientemente de los valores de las variables de control. Sin embargo, este no siempre es el caso.\nNota 3 Sin nuestro modelo \\(p(y|t,z)\\) es lineal, y nos interesa calcular el efecto causal promedio de la variable \\(t\\), no es necesario promediar por la conjunta de \\(p(z)\\). Bajo estas condiciones, el efecto causal promedio está simplemente por el coeficiente de \\(t\\) en el modelo lineal. Sin embargo, si este no es el caso, entonces para estimar el efecto causal promedio es necesario promediar apropiadamente según la fórmula de ajuste.",
+ "crumbs": [
+ "6Identificación y cálculo-do"
+ ]
+ },
+ {
+ "objectID": "06-calculo-do.html#bloqueando-puertas-traseras",
+ "href": "06-calculo-do.html#bloqueando-puertas-traseras",
+ "title": "6 Identificación y cálculo-do",
+ "section": "6.5 Bloqueando puertas traseras",
+ "text": "6.5 Bloqueando puertas traseras\nEn las partes anteriores vimos que estratificando por los padres de la variable de tratamiento \\(X\\) podemos construir un estimador del efecto de \\(X\\) sobre otra variable \\(Y\\), pasando de una distribución observacional a una conceptualmente experimental (dado que los supuestos causales sean aproximadamente correctos).\nSin embargo, esta aplicación de la fórmula de ajuste no funciona si existen padres que no fueron observados, y por tanto no podemos estratificar por ellos. El siguiente método (ajuste por “puerta trasera”) nos da una técnica adicional que podemos usar dado ciertos tipos de estructura en nuestro modelo causal, y presenta una mejoría sobre la fórmula de ajuste simple (veremos también por ejemplo, que a veces podemos usar menos variables que padres de la variable de interés). Nótese que una vez más, este criterio sólo depende de la gráfica causal \\(G\\) asociada a nuestro modelo, y no los modelos locales que utilizemos para modelar la condicional de cada nodo.\n\n\n\n\n\n\nAjuste de puerta trasera (Pearl)\n\n\n\nSi tenemos dos variables \\(T\\) y \\(Y\\) en una gráfica \\(G\\), un conjunto \\(Z\\) de variables satisface el criterio de puerta trasera relativo a \\(T\\) y \\(Y\\) cuando \\(Z\\) bloquea cualquier camino entre \\(T\\) y \\(Y\\) que tenga una arista que incida en \\(T\\), y ninguna variable de \\(Z\\) es descendiente de \\(T\\).\nEn tal caso, podemos utilizar la fórmula de ajuste, pero en lugar de estratificar por los padres de \\(T\\), estratificamos por las variables en \\(Z\\)\n\n\nLa idea es:\n\nQueremos bloquear todos los caminos no causales entre \\(T\\) y \\(Y\\).\nQueremos no perturbar todos los caminos dirigidos de \\(T\\) a \\(Y\\) (caminos causales).\nNo queremos activar caminos no causales entre \\(T\\) y \\(Y\\) al condicionar.\n\nCumplimos 1 al estratificar por variables que bloquean los caminos que son causas de \\(T\\), pues estos caminos no son causales y distorsionan la relación entre \\(T\\) y \\(Y\\). Al mismo tiempo, no bloqueamos caminos causales porque ningúna variable de \\(Z\\) es descendiente de \\(T\\), de modo que se satisface el criterio 2 (todos los caminos causales comienzan con \\(T\\to\\)). Finalmente, al excluir descendientes de \\(T\\) también implica que no condicionamos a colisionadores del tipo \\(T\\to \\cdots \\to Z_1\\gets Y\\), pues esto activa un camino no causal entre \\(T\\) y \\(Y\\) (se cumple 3).\n\nEjemplo (Pearl)\nConsideramos primero este ejemplo simple, donde queremos evaluar la efectividad de un tratamiento en cierta enfermedad. Los datos que tenemos disponibles son si una persona recibió 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. Nuestros supuestos están en la siguiente gráfica:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2, rankdir = LR]\n node [shape=plaintext]\n Trata\n Res\n node [shape = circle]\n NSE\n Peso\n U\n edge [minlen = 3]\n NSE -> Peso\n NSE -> Trata\n Trata -> Res\n Peso -> Res\n U -> NSE\n U -> Peso\n}\n\")\n\n\n\n\n\n\nObservamos que no podemos directamente usar la fórmula de ajuste pues NSE no es una variable observada.\nEn esta circunstancia no podríamos identificar el efecto causal, pues existen un caminos abiertos no causales. Quizá el tratamiento no es muy efectivo, y parece ser bueno pues fue aplicado a personas con menor peso que las que no recibieron el tratamiento, a través del efecto de NSE. Sin embargo, supón que tuviéramos disponible la variable Peso:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2, rankdir = LR]\n node [shape=plaintext]\n Trata\n Res\n Peso\n node [shape = circle]\n NSE\n U\n edge [minlen = 3]\n NSE -> Peso\n NSE -> Trata\n Trata -> Res\n Peso -> Res\n U -> NSE\n U -> Peso\n}\n\")\n\n\n\n\n\n\nEn este caso, todavía no podemos aplicar la fórmula original de ajuste pues no conocemos \\(NSE\\). Sin embargo, podemos bloquear los caminos no causales estratificando por Peso, y entonces podemos usar el criterio de puerta trasera para identificar el efecto del tratamiento, aún cuando no tengamos NSE.\n\n\nEjemplo\nPrimero consideramos un modelo generador:\n\ninv_logit <- function(x) 1 / (1 + exp(-x))\nsimular_bd <- function(n = 10){\n nse <- sample(c(0, 1), n, replace = TRUE)\n peso <- rnorm(n, 70 - 7 * nse, 12 + 2 * nse)\n trata <- rbinom(n, 1, 0.8 * nse + 0.2 * (1 - nse))\n p_trata <- inv_logit(1 * trata - 0.2 * (peso - 70))\n res <- rbinom(n, 1, p_trata)\n tibble(nse, peso, trata, res)\n}\ndatos_bd <- simular_bd(10000)\nhead(datos_bd)\n\n# A tibble: 6 × 4\n nse peso trata res\n <dbl> <dbl> <int> <int>\n1 1 71.9 0 0\n2 0 45.0 0 1\n3 0 73.5 0 0\n4 0 66.1 0 1\n5 1 49.4 1 1\n6 0 69.0 1 1\n\n\nVeamos qué sucede si cruzamos tratamiento con resultado (es una muestra grande y el error de estimación no es importante):\n\ndatos_bd |> \n count(trata, res) |>\n group_by(trata) |> \n mutate(p = n / sum(n)) |> \n filter(res == 1) |> \n ungroup() |> \n mutate(dif = p - lag(p))\n\n# A tibble: 2 × 5\n trata res n p dif\n <int> <int> <int> <dbl> <dbl>\n1 0 1 2678 0.533 NA \n2 1 1 3686 0.741 0.208\n\n\nSabemos que esta diferencia en respuesta puede estar confundida por un camino no causal. El verdadero efecto casual podemos calcularlo en nuestras simulaciones como sigue a partir de nuestro modelo (igualmente, usamos una muestra muy grande):\n\nsimular_efecto <- function(n = 10, peso = NULL){\n # cómo es la población\n nse <- sample(c(0, 1), n, replace = TRUE)\n if(is.null(peso)){\n peso <- rnorm(n, 70 - 7 * nse, 12 + 2 * nse)\n }\n # asignar al azar\n trata <- rbinom(n, 1, 0.5)\n p_trata <- inv_logit(1 * trata - 0.2 * (peso - 70))\n res <- rbinom(n, 1, p_trata)\n tibble(nse, peso, trata, res)\n}\nsims_efecto <- simular_efecto(20000)\nresumen <- sims_efecto |> \n count(trata, res) |>\n group_by(trata) |> \n mutate(p = n / sum(n)) |> \n filter(res == 1) |> \n ungroup() |> \n mutate(dif = p - lag(p))\ndif_real <- resumen$dif[2]\nresumen\n\n# A tibble: 2 × 5\n trata res n p dif\n <int> <int> <int> <dbl> <dbl>\n1 0 1 5929 0.590 NA \n2 1 1 6996 0.703 0.113\n\n\nLa estimación ingenua del cruce simple es mucho más grande que el verdadero efecto.\nPodemos también calcular el efecto para un peso particular:\n\nsims_efecto <- simular_efecto(20000, peso = 70)\nres_70 <- sims_efecto |> \n count(trata, res) |>\n group_by(trata) |> \n mutate(p = n / sum(n)) |> \n filter(res == 1) |> \n ungroup() |> \n mutate(dif = p - lag(p))\ndif_70 <- res_70$dif[2]\nres_70\n\n# A tibble: 2 × 5\n trata res n p dif\n <int> <int> <int> <dbl> <dbl>\n1 0 1 5002 0.500 NA \n2 1 1 7344 0.735 0.235\n\n\nSuponiendo nuestro diagrama, queremos estimar estratificando por peso. Podríamos usar un sólo modelo logístico, pero pueden ser más simples los cálculos si construimos nuestro modelo en stan. En este caso, podríamos calcular las diferencias para un peso particular, por ejemplo 70 kg (en lugar de modelar estaturas para producir una estimación de diferencia promedio).\nUsaremos una muestra de 2 mil personas:\n\nmod_trata <- cmdstan_model(\"./src/trata-backdoor.stan\")\nprint(mod_trata)\n\ndata {\n int<lower=0> N;\n vector[N] trata;\n array[N] int res;\n vector[N] peso;\n\n}\n\ntransformed data {\n real media_peso;\n\n // centrar\n media_peso = mean(peso);\n}\n\nparameters {\n real gamma_0;\n real gamma_1;\n real gamma_2;\n}\n\ntransformed parameters {\n vector[N] p_logit_res;\n\n p_logit_res = gamma_0 + gamma_1 * trata + gamma_2 * (peso - media_peso);\n\n}\n\nmodel {\n // modelo de resultado\n res ~ bernoulli_logit(p_logit_res);\n gamma_0 ~ normal(0, 2);\n gamma_1 ~ normal(0, 1);\n gamma_2 ~ normal(0, 0.2);\n\n\n}\ngenerated quantities {\n real dif_trata;\n real p_trata;\n real p_no_trata;\n\n real peso_sim = 70;\n {\n array[2000] int res_trata;\n array[2000] int res_no_trata;\n for(k in 1:2000){\n res_trata[k] = bernoulli_rng(\n inv_logit(gamma_0 + gamma_1 * 1 +\n gamma_2 * (peso_sim - media_peso)));\n res_no_trata[k] = bernoulli_rng(\n inv_logit(gamma_0 + gamma_1 * 0 +\n gamma_2 * (peso_sim - media_peso)));\n }\n p_trata = mean(res_trata);\n p_no_trata = mean(res_no_trata);\n }\n dif_trata = p_trata - p_no_trata;\n}\n\n\n\nset.seed(915)\ndatos_bd <- simular_bd(2000)\ndatos_lista <- list(N = nrow(datos_bd),\n trata = datos_bd$trata, res = datos_bd$res,\n peso = datos_bd$peso)\najuste <- mod_trata$sample(data = datos_lista, refresh = 1000)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 1 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 1 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 1 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 1 finished in 1.9 seconds.\nChain 2 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 2 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 2 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 2 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 2 finished in 1.9 seconds.\nChain 3 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 3 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 3 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 3 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 3 finished in 1.9 seconds.\nChain 4 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 4 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 4 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 4 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 4 finished in 2.0 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 1.9 seconds.\nTotal execution time: 8.2 seconds.\n\nsims <- ajuste$draws( format = \"df\")\nresumen <- ajuste$summary(c( \"dif_trata\"))\n\n\nresumen |> select(variable, mean, q5, q95)\n\n# A tibble: 1 × 4\n variable mean q5 q95\n <chr> <dbl> <dbl> <dbl>\n1 dif_trata 0.214 0.162 0.268\n\nsims |> select(dif_trata) |> \n ggplot(aes(x = dif_trata)) + geom_histogram() +\n geom_vline(xintercept = dif_70, colour = \"red\")\n\nWarning: Dropping 'draws_df' class as required metadata was removed.\n\n\n`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.\n\n\n\n\n\n\n\n\n\nY obtenemos una estimación correcta del efecto en 70 kg. Podríamos también calcular el efecto en distintos pesos (nuestro estimador es una curva), promediar estimando una distribución de pesos modelada, o tomar una distribución fija de pesos para modelar (cada una de estas estrategias tiene propósitos diferentes).\nSi queremos tener un efecto promedio, podemos modelar los pesos. Otra estrategia es promediar sobre los valores observados de la muestra. Nótese que esto ignora una parte de la incertidumbre proveniente de la muestra particular usada.\n\nmod_trata <- cmdstan_model(\"./src/trata-backdoor-promedio.stan\")\nprint(mod_trata)\n\ndata {\n int<lower=0> N;\n vector[N] trata;\n array[N] int res;\n vector[N] peso;\n\n}\n\ntransformed data {\n real media_peso;\n\n // centrar\n media_peso = mean(peso);\n}\n\nparameters {\n real gamma_0;\n real gamma_1;\n real gamma_2;\n}\n\ntransformed parameters {\n vector[N] p_logit_res;\n\n p_logit_res = gamma_0 + gamma_1 * trata + gamma_2 * (peso - media_peso);\n\n}\n\nmodel {\n // modelo de resultado\n res ~ bernoulli_logit(p_logit_res);\n gamma_0 ~ normal(0, 2);\n gamma_1 ~ normal(0, 1);\n gamma_2 ~ normal(0, 0.2);\n\n\n}\ngenerated quantities {\n real dif_trata;\n real p_trata;\n real p_no_trata;\n vector[N] probs;\n\n for(i in 1:N){\n probs[i] = 1.0 / N;\n }\n\n {\n array[2000] int res_trata;\n array[2000] int res_no_trata;\n for(k in 1:2000){\n real peso_sim = peso[categorical_rng(probs)];\n res_trata[k] = bernoulli_rng(\n inv_logit(gamma_0 + gamma_1 * 1 +\n gamma_2 * (peso_sim - media_peso)));\n res_no_trata[k] = bernoulli_rng(\n inv_logit(gamma_0 + gamma_1 * 0 +\n gamma_2 * (peso_sim - media_peso)));\n }\n p_trata = mean(res_trata);\n p_no_trata = mean(res_no_trata);\n }\n dif_trata = p_trata - p_no_trata;\n\n}\n\n\n\ndatos_lista <- list(N = nrow(datos_bd),\n trata = datos_bd$trata, res = datos_bd$res,\n peso = datos_bd$peso)\najuste <- mod_trata$sample(data = datos_lista, refresh = 1000)\n\nRunning MCMC with 4 sequential chains...\n\nChain 1 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 1 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 1 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 1 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 1 finished in 10.9 seconds.\nChain 2 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 2 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 2 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 2 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 2 finished in 10.9 seconds.\nChain 3 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 3 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 3 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 3 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 3 finished in 10.9 seconds.\nChain 4 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 4 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 4 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 4 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 4 finished in 10.9 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 10.9 seconds.\nTotal execution time: 43.9 seconds.\n\nsims <- ajuste$draws(c(\"dif_trata\"), format = \"df\")\n\n\nresumen <- ajuste$summary(c( \"dif_trata\"))\nresumen |> select(variable, mean, q5, q95)\n\n# A tibble: 1 × 4\n variable mean q5 q95\n <chr> <dbl> <dbl> <dbl>\n1 dif_trata 0.111 0.0805 0.141\n\nsims |> select(dif_trata) |> \n ggplot(aes(x = dif_trata)) + geom_histogram() +\n geom_vline(xintercept = dif_real, colour = \"red\")\n\nWarning: Dropping 'draws_df' class as required metadata was removed.\n\n\n`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.\n\n\n\n\n\n\n\n\n\nY recuperamos nuevamente el efecto verdadero que mostramos arriba.",
+ "crumbs": [
+ "6Identificación y cálculo-do"
+ ]
+ },
+ {
+ "objectID": "06-calculo-do.html#reglas-del-cálculo-do-opcional",
+ "href": "06-calculo-do.html#reglas-del-cálculo-do-opcional",
+ "title": "6 Identificación y cálculo-do",
+ "section": "6.6 Reglas del cálculo-do (opcional)",
+ "text": "6.6 Reglas del cálculo-do (opcional)\nExisten tres axiomas básicos del cálculo-do de las que se derivan los demás resultados, como veremos en el siguiente ejemplo del criterio de la puerta delantera.\nAntes de verlas, un resumen rápido de las reglas es el siguiente:\n\nLa regla 1 nos dice que las distribuciones asociadas a intervenciones satisfacen también la equivalencia de \\(d\\)-separación e independencia condicional: si \\(Y\\) y \\(Z\\) están \\(d\\)-separadas dado en la gráfica manipulada, entonces \\(p(y | do(x), z) = p(y|do(x))\\).\nLa regla 2 es el criterio de la puerta trasera: si condicionamos a variables \\(W\\) que bloquean toda puerta trasera de \\(X\\) a \\(Y\\), podemos cambiar \\(do(x)\\) por \\(x\\): \\(p(y | do(x), w) = p(y | x, w)\\).\nLa regla 3 expresa que si no hay caminos causales de \\(X\\) a \\(Y\\), entonces \\(p(y|do(x)) = p(y)\\).\n\n\n\n\n\n\n\nCompletitud (Shpitser, Pearl)\n\n\n\nSi un efecto causal es identificable (puede expresarse en términos de cantidades observacionales), entonces puede derivarse una estrategia de identificación a partir de las tres reglas del cálculo-do.\n\n\nNota: esto no excluye que bajo ciertas hipótesis adicionales a las de nuestra gráfica causal (por ejemplo cómo se comportan las distribuciones particulares qeu componen el modelo), sea posible identificar efectos causales con otros medios que van más allá del cálculo-do.\nCon más generalidad, abajo están estas reglas (donde condicionamos a más variables o hacemos más intervenciones, y afinamos las condiciones):\nDenotamos por \\(G_m\\) la gráfica mutilada por \\(do(x)\\), donde quitamos todas las aristas que entran en \\(X\\). Los tres axiomas son:\nRegla 1 Ignorar observaciones: Si \\(Y\\) y \\(Z\\) están \\(d\\)-separados por \\(X\\) y \\(W\\) en \\(G_m\\),\n\\[ p(y|do(x), z, w) = p(y|do(x), w)\\] O en otras palabras, si \\(p_m\\) es la conjunta para \\(G_m\\),\n\\[p_m(y|x,z,w) = p_m(y|x, w)\\] es cierto si \\(Y\\) y \\(Z\\) están \\(d\\)-separados por \\(X\\) y \\(W\\) en \\(G_m\\) (condicionalmente independientes). Así que esta regla es independencia condicional dado \\(d\\)-separación, pero para la gráfica intervenida.\nRegla 2 Usando observaciones como intervenciones:\nSi \\(Y\\) y \\(Z\\) están \\(d\\)-separados por \\(X\\) y \\(W\\) en \\(G_m\\) quitándole todas las aristas que salen de \\(Z\\), entonces\n\\[ p(y|do(x), do(z), w) = p(y|do(x), z, w)\\] Regla 3 Ignorar intervenciones:\nSi \\(Z\\) y \\(Y\\) están \\(d\\)-separadas por \\(X\\) y \\(W\\) en la gráfica \\(G_m\\) donde además quitamos cualquier arista a \\(Z\\) si \\(Z\\) no es antecesor de \\(W\\) en \\(G_m\\), entonces:\n\\[ p(y|do(x), do(z), w) = p(y|do(x), w)\\]",
+ "crumbs": [
+ "6Identificación y cálculo-do"
+ ]
+ },
+ {
+ "objectID": "06-calculo-do.html#el-criterio-de-puerta-delantera",
+ "href": "06-calculo-do.html#el-criterio-de-puerta-delantera",
+ "title": "6 Identificación y cálculo-do",
+ "section": "6.7 El criterio de puerta delantera",
+ "text": "6.7 El criterio de puerta delantera\nEn algunos casos, puede ser que no sea posible bloquear algún camino no causal con variables observadas. Un ejemplo clásico es el de la discusión acerca de la relación de fumar con cáncer de pulmón. Algunos estadísticos plantearon que los estudios de asociación entre fumar y cáncer de pulmón podrían tener efectos gravemente confundidos, por ejemplo, por aspectos genéticos que hacen a una persona propensa a fumar al mismo tiempo que aumenta su probabilidad de fumar:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n F\n C\n node [shape = circle]\n U\n edge [minlen = 3]\n U -> F\n U -> C\n F -> C\n{rank= same; C; F}\n}\n\")\n\n\n\n\n\n\nEn este caso, el efecto de fumar (\\(F\\)) sobre cáncer (\\(C\\)) no es identificable pues no podemos condicionar a la variable de Genotipo (\\(U\\)). Supongamos que tenemos una medida adicional, que es la cantidad de depósitos de alquitrán den los pulmones de los pacientes. Este es es afectado por \\(F\\), y a su vez, el alquitrán incrementa la probabilidad de cáncer:\n\n\nCódigo\ngrViz(\"\ndigraph {\n graph [ranksep = 0.2]\n node [shape=plaintext]\n F\n C\n A\n node [shape = circle]\n U\n edge [minlen = 3]\n U -> F\n U -> C\n F -> A\n A -> C\n{rank= same; C; F; A}\n}\n\")\n\n\n\n\n\n\nLa idea es primero estimar el efecto de \\(F\\) sobre \\(A\\), y después estimar el efecto de \\(A\\) sobre \\(C\\). La “composición” de estos dos efectos, dado el diagrama, debe darnos el estimador correcto. Primero consideramos el efecto de \\(F\\) sobre \\(A\\), y tenemos que (regla 2)\n\\[p(a|do(f)) = p(a|f),\\] La igualdad se debe a que una vez que condicionamos a \\(F\\) no hay puertas traseras entre \\(F\\) y \\(A\\) (pues no condicionamos a \\(C\\)). Esta dependencia causal la podemos entonces estimar de los datos.\nEl efecto de \\(A\\) sobre \\(C\\) también es identificable, pues el camino no causal se bloquea cuando condicionamos a \\(A\\), de forma que por la fórmula de ajuste:\n\\[p(c|do(a)) = \\int p(c|a, f') p(f')\\, df'\\]\nAhora encadenamos estas dos ecuaciones:\n\\[p(c|do(f)) = \\int p(c|do(a))p(a|f)\\,da\\]\nque equivale en simulación a: dado un valor de \\(F\\), simulamos \\(A=a\\) con nuestro modelo ajustado con datos naturales. Ahora intervenimos \\(A\\) con el valor a que obtuvimos y simulamos \\(C\\). Sin embargo, para hacer este último paso con datos naturales, necesitamos usar el criterio de puerta trasera como explicamos arriba: simulamos entonces \\(f´\\) de \\(p(f)\\), y después simulamos \\(C\\) en función de \\(a\\) y \\(f´\\) (con una distribución construida a partir de datos).\nRequerimos en este caso construir y estimar la condicional \\(p(c|a, f)\\) basado en los datos.\nEn fórmula, en general, se escribe como:\n\n\n\n\n\n\nCriterio de fuerta delantera (Pearl)\n\n\n\nDecimos que un conjunto de variables \\(A\\) satisface el criterio de puerta delantera en relación a las variables \\(F\\) y \\(C\\) cuando:\n\n\\(A\\) intercepta todos las cadenas dirigidos de \\(F\\) a \\(C\\)\nNo hay ningún camino activo de puerta trasera de \\(F\\) a \\(A\\)\nTodos los caminos de puerta trasera de \\(A\\) a \\(C\\) están bloqueados por \\(F\\).\n\nSi \\(A\\) satisface el criterio de puerta delantera en relación a \\(F\\) y \\(C\\), entonces el efecto causal de \\(F\\) en \\(C\\) es identificable y está dado por la fórmula:\n\\[p(c|do(f)) = \\int \\left [ \\int p(c|a,f´)p(f´)\\,df´ \\right ] p(a|f)\\,da\\]\n\n\nTodas estas cantidades puede estimarse de los datos.\n\nEjemplo: proceso generador\nAntes de aplicar este nuevo procedimiento, describamos el proceso generador que utilizaremos:\n\n# simular distribución natural\nsimular_fd <- function(n = 10, efecto_a = 0.3){\n ## causa común\n u <- rnorm(n, 0, 1);\n # cantidad que fuma\n f <- exp(rnorm(n, 1 + 0.2 * u, 0.1))\n # acumulación de alquitrán\n a <- rnorm(n, 4 * f, 2)\n # probabilidad de cancer\n p_c <- inv_logit(-6 + efecto_a * a + 2 * u)\n c <- rbinom(n, 1, p_c)\n tibble(f, a, c, u)\n}\n# simular datos intervenidos (suponiendo que conocemos todo)\nsim_int_f <- function(n = 100, do_f = 0.3, efecto_a = 0.3){\n a <- rnorm(n, 4 * do_f, 2)\n u <- rnorm(n, 0, 1)\n p_c <- inv_logit(-6 + efecto_a * a + 2 * u)\n c <- rbinom(n, 1, p_c)\n tibble(do_f = do_f, media_c = mean(c))\n}\n\n\nset.seed(4481)\nsims_fd <- simular_fd(5000)\nsims_fd_1 <- simular_fd(10000)\nqplot(sims_fd$f, sims_fd$a)\n\nWarning: `qplot()` was deprecated in ggplot2 3.4.0.\n\n\n\n\n\n\n\n\n\n¿Cómo se ve la relación de fumador con cáncer? En esta gráfica mostramos también el valor de la variable no observada \\(U\\). Nótese que parte de la correlación positiva que existe es debido a esta variable \\(U\\).\n\nggplot(sims_fd, aes(x = f, y = c, colour = u)) + \n geom_jitter() + scale_colour_continuous(type = \"viridis\")\n\n\n\n\n\n\n\n\nAhora veamos cómo se ve el efecto de \\(F\\) sobre \\(C\\) y también cómo se ve el cruce de \\(F\\) y \\(C\\) en los datos naturales:\n\nsims_1 <- map_df(seq(1, 4, 0.5), ~ sim_int_f(100000, .x))\n\nsims_1 |> \n ggplot() + geom_line(aes(x = do_f, y = media_c)) +\n geom_smooth(data = sims_fd_1, aes(x = f, y = c), method = \"loess\", span = 0.3, se = FALSE, colour =\"red\") + xlab(\"Grado de tabaquismo\") +\n xlim(c(1,4))\n\n`geom_smooth()` using formula = 'y ~ x'\n\n\nWarning: Removed 376 rows containing non-finite values (`stat_smooth()`).\n\n\n\n\n\n\n\n\n\nEn efecto causal promedio de fumar, en cada nivel, sobre la incidencia de cáncer de pulmón, suponiendo nuestro proceso generador. Nótese que la relación no es tan fuerte como observamos en los datos naturales (en rojo). Esto se debe a que en los datos naturales, las personas existe una causa común entre no fumar y prevenir cáncer de pulmón.\n\n\nEjemplo: estimación con puerta delantera\nVeamos cómo sería la estimación si tuviéramos datos disponible, y si es que podemos recuperar el efecto correcto dados los datos observados y la técnica de puerta delantera.\nNótese que sólo necesitamos \\(p(c|a, f), p(a|f)\\) y \\(p(f)\\). Estos son modelos estadísticos con el que podemos identificar el efecto que nos interesa. Una vez que los estimemos, podemos usar simulación:\n\nFijamos una \\(f\\).\nSimulamos una \\(a\\) del modelo \\(p(a|f)\\)\nPara calcular \\(\\int p(c|a,f')p(f')\\), tenemos que simular un valor \\(f'\\) de la marginal de \\(p(f)\\), y luego, sustituir junto la \\(a\\) de 1 para simular una \\(c\\) de \\(p(c|a, f')\\).\nConsideramos solamente \\(c\\) y \\(f\\) para resumir el efecto.\n\n\nset.seed(481)\nsims_fd <- simular_fd(2000)\nmod_front_door <- cmdstan_model(\"./src/front-door.stan\")\nprint(mod_front_door)\n\ndata {\n int<lower=0> N;\n int<lower=0> n_f;\n vector[N] f;\n vector[N] a;\n array[N] int<lower=0, upper=1> c;\n array[n_f] real do_f;\n\n}\n\ntransformed data {\n real media_a;\n real media_f;\n\n media_a = mean(a);\n media_f = mean(f);\n}\n\nparameters {\n real<lower=0> alpha;\n real alpha_a;\n real<lower=0> alpha_f;\n real int_a;\n real beta_0;\n real<lower=0> beta_1;\n real<lower=0> beta;\n real<lower=0> a_f;\n real<lower=0> b_f;\n real<lower=0> sigma_a;\n real<lower=0> sigma_f;\n\n}\n\n\n\ntransformed parameters {\n\n\n}\n\nmodel {\n f ~ gamma(a_f, b_f);\n a ~ normal(beta * f, sigma_a);\n c ~ bernoulli_logit(int_a + alpha_a * a + alpha_f * f);\n alpha_a ~ normal(0, 1);\n alpha_f ~ normal(0, 1);\n int_a ~ normal(0, 3);\n sigma_a ~ normal(0, 1);\n sigma_f ~ normal(0, 0.1);\n alpha ~ normal(0, 1);\n beta ~ normal(0, 1);\n beta_0 ~ normal(0, 3);\n beta_1 ~ normal(0, 1);\n\n}\ngenerated quantities {\n array[n_f] real mean_c;\n\n for(i in 1:n_f){\n array[2000] real res_sim;\n for(j in 1:2000){\n real a_sim = normal_rng(beta * (do_f[i]), sigma_a);\n real f_sim = gamma_rng(a_f, b_f);\n res_sim[j] = bernoulli_rng(inv_logit(int_a + alpha_a * a_sim + alpha_f * f_sim));\n }\n mean_c[i] = mean(res_sim);\n }\n\n}\n\n\n\ndo_f <- seq(1, 4, 0.1)\nn_f <- length(do_f)\nsims <- mod_front_door$sample(data = list(N = nrow(sims_fd),\n f = sims_fd$f, a = sims_fd$a,\n c = sims_fd$c, do_f = do_f, n_f = n_f),\n init = 0.01, step_size = 0.01, \n refresh = 1000,\n parallel_chains = 4)\n\nRunning MCMC with 4 parallel chains...\n\nChain 1 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 2 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 3 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 4 Iteration: 1 / 2000 [ 0%] (Warmup) \nChain 4 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 4 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 3 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 3 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 1 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 1 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 2 Iteration: 1000 / 2000 [ 50%] (Warmup) \nChain 2 Iteration: 1001 / 2000 [ 50%] (Sampling) \nChain 4 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 4 finished in 42.8 seconds.\nChain 3 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 3 finished in 43.1 seconds.\nChain 2 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 2 finished in 44.3 seconds.\nChain 1 Iteration: 2000 / 2000 [100%] (Sampling) \nChain 1 finished in 44.5 seconds.\n\nAll 4 chains finished successfully.\nMean chain execution time: 43.7 seconds.\nTotal execution time: 44.7 seconds.\n\n\n\nsims_efecto_tbl <- sims$draws(\"mean_c\", format = \"df\") |> \n pivot_longer(cols = contains(\"mean_c\"), values_to = \"media_c\") |> \n separate(name, c(\"nom\", \"id\"), \n sep = \"[\\\\[\\\\]]\", convert = TRUE, extra = \"drop\") |> \n left_join(tibble(f = do_f) |> \n mutate(id = seq_along(f))) \nresumen_tbl <- sims_efecto_tbl |> \n group_by(id, f) |> \n summarise(media = mean(media_c), \n q5 = quantile(media_c, 0.05),\n q95 = quantile(media_c, 0.95))\n\n\nggplot(resumen_tbl) + \n geom_linerange(aes(x= f, ymax = q95, ymin = q5), colour = \"red\") + \n geom_point(aes(x = f, y = media), colour = \"red\") +\n geom_line(data = sims_1, aes(x = do_f, y = media_c)) +\n xlab(\"Nivel de tabaquismo\") + ylab(\"Prop afectada\")\n\n\n\n\n\n\n\n\nY parece que hemos obtenido una estimación razonable del efecto causal de fumar sobre cáncer. Recordemos también que debemos ser cuidadosos al comparar intervalos que salen del mismo modelo por su nivel de traslape.\nPor ejemplo, si quisiéramos calcular contrastes con el nivel 2 de tabaquismo:\n\nefecto_2 <- sims_efecto_tbl |> filter(f == 2) |> \n select(.draw, efecto_2 = media_c)\ncomp_tbl <- left_join(sims_efecto_tbl, efecto_2) |> \n mutate(dif_2 = media_c - efecto_2)\n\nJoining with `by = join_by(.draw)`\n\ncomp_tbl |> group_by(f) |> \n summarise(media = mean(dif_2), q5 = quantile(dif_2, 0.05),\n q95 = quantile(dif_2, 0.95)) |> \nggplot() + geom_linerange(aes(x= f, ymax = q95, ymin = q5)) + geom_point(aes(x = f, y = media)) +\n xlab(\"Nivel de tabaquismo\") + ylab(\"Prop afectada\")\n\n\n\n\n\n\n\n\nNota: nótese como en este ejemplo hemos evitado incluir en nuestro modelo la variable no observada \\(U\\), gracias al procedimiento de puerta delantera descrito arriba.\nEs posible sin embargo intentar un modelo completo bayesiano, sin necesidad de recordar la fórmula. El procedimiento, que es más difícil de ajustar: considera una variable latente \\(U\\) no observada, y es necesario definir cómo puede ser su relación con sus descendientes. Es necesario más cuidado en definir formas funcionales e iniciales apropiadas para que los muestreadores funcionen apropiadamente.",
+ "crumbs": [
+ "6Identificación y cálculo-do"
+ ]
}
]
\ No newline at end of file
