El bueno, el feo y el bayesiano

Listening

Introducción

Cuenta Matheus Facure Alves en este capítulo de Causal Inference for the Brave and True que hay buenos y malos controles en esto de la inferencia causal. Y no le falta razón, tenemos al bueno y al feo, pero yo quiero hablar del otro, no del malo, sino del bayesiano

Ejemplo

Para ilustrar lo de buenos y malos (yo llamaré feos) confounders Matheus usa unos datos simulados. Yo voy a usar el mismo proceso de generación de datos.

Usa el ejemplo de que eres un científico de datos en un equipo de recobro en una fintech. Y que tu tarea es estimar el impacto que tiene enviar un mail pidiendo a la gente que negocie su deuda. La variable repuesta es el importe de por pagar de los clientes morosos.

El conjunto de datos luce tal que así

#  algunas funciones para configurar ggplot y gggdag. El que quiera que visite mi github

source(here::here("aux_functions/ggdag-mask.R"))
source(here::here("aux_functions/setup.R"))
library(tidyverse)

library(dagitty)
library(ggdag)
library(ggokabeito)

library(cmdstanr)
library(brms)
library(posterior) # para cosas como rvar

options(brms.backend = "cmdstanr")

La generación de los datos, copiada de Matheus es el siguiente

semilla <- 44
n <-  5000

email <- rbinom(n = n, size = 1, prob = 0.5)
credit_limit <-  rgamma(n = n, shape = 6, scale = 200)
risk_score <- rbeta(n= n, shape1 = credit_limit, shape2 = mean(credit_limit) )

opened <- rnorm(n = n , mean = 5  + 0.001 * credit_limit - 4 * risk_score, sd = 2 )
opened <- (opened > 4) * email

agreement <- rnorm(n = n, mean = 30 +(-0.003*credit_limit - 10*risk_score), sd = 7) * 2 * opened
agreement <- as.integer(agreement  > 40)

payments <- floor(rnorm(n = n , mean = 500 + 0.16*credit_limit - 40*risk_score + 11*agreement + email, sd = 75) / 10) *10

data <-  data.frame( payments, email, opened, agreement, credit_limit, risk_score)

write_csv(data, here::here("data/collections_email.csv"))

El bueno y sus amigos

Sabemos que email ha sido generado aleatoriamente y por tanto podríamos estimar el efecto causal de forma simple

m1 <- lm(payments ~ email, data = data)
summary(m1)

#> 
#> Call:
#> lm(formula = payments ~ email, data = data)
#> 
#> Residuals:
#>    Min     1Q Median     3Q    Max 
#>   -312    -72     -2     68    511 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   669.17       2.13  314.59   <2e-16 ***
#> email           2.86       2.99    0.96     0.34    
#> ---
#> Signif. codes:  
#> 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 106 on 4998 degrees of freedom
#> Multiple R-squared:  0.000184,   Adjusted R-squared:  -1.62e-05 
#> F-statistic: 0.919 on 1 and 4998 DF,  p-value: 0.338

confint(m1, level = 0.89)

#>              5.5 %  94.5 %
#> (Intercept) 665.77 672.567
#> email        -1.91   7.633

Vaya, resulta que la estimación no es muy precisa.

Podríamos pensar que las variables risk_score y credit_limit influyen en payments, al fin y al cabo, no se le va a dar un préstamo muy alto a quien tenga un alto risk_score. Así el dag podría ser este

email_dag <- dagify(
  payments ~ email + credit_limit + risk_score,
  exposure = "email",
  outcome = "payments",
  labels = c(
    payments = "payments",
    email = "Envio mail",
    credit_limit = "credit_limit",
    risk_score = "risk_score"
  ),
  coords = list(
    x = c(
      payments = 5,
      email = 4,
      credit_limit = 4,
      risk_score = 4
    ),
    y = c(
      payments = 2,
      email = 2,
      credit_limit = 3,
      risk_score = 1
    ))
)

email_dag |>
  tidy_dagitty() |>
  node_status() |>
  ggplot(
    aes(x, y, xend = xend, yend = yend, color = status)
  ) +
  geom_dag_edges() +
  geom_dag_point() +
  geom_dag_label_repel() +
  scale_color_okabe_ito(na.value = "grey90") +
  theme_dag() +
  theme(legend.position = "none") +
  coord_cartesian(clip = "off")

Pues podríamos considerar que los amigos del bueno serían estas dos variables, risk_score y credit_limit. Que sabemos que no son necesarias para obtener una estimación insesgada del efecto causal, puesto que no son variables de confusión. No obstante son amigos del bueno porque al estar relacionadas con la respuesta (si nuestro modelo causal es correcto), incrementan la precisión de la estimación , tal y como vemos en el resultado del modelo siguiente.

Amigos frecuentistas

m2 <- lm(payments ~ email + credit_limit + risk_score , data = data)
summary(m2)

#> 
#> Call:
#> lm(formula = payments ~ email + credit_limit + risk_score, data = data)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -275.21  -50.39    0.36   51.11  251.63 
#> 
#> Coefficients:
#>               Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   522.7745     9.4144   55.53  < 2e-16 ***
#> email           5.1121     2.1308    2.40  0.01647 *  
#> credit_limit    0.1784     0.0077   23.17  < 2e-16 ***
#> risk_score   -143.0931    37.0293   -3.86  0.00011 ***
#> ---
#> Signif. codes:  
#> 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 75.3 on 4996 degrees of freedom
#> Multiple R-squared:  0.491,  Adjusted R-squared:  0.491 
#> F-statistic: 1.61e+03 on 3 and 4996 DF,  p-value: <2e-16

confint(m2, level = 0.89)

#>                  5.5 %   94.5 %
#> (Intercept)   507.7259 537.8232
#> email           1.7061   8.5183
#> credit_limit    0.1661   0.1907
#> risk_score   -202.2836 -83.9025

Amigos bayesianos

m4 <- brm(
  payments ~ email + credit_limit + risk_score ,
  data = data,
  seed = 44,
  chains = 4,
  iter = 3000,
  warmup = 1000,
  cores = 4,
  file = here::here("brms_stan_models/email2")
)

# hqy que evitar el exceso de decimales dando impresión de falsa exactitud

round(posterior_summary(m4), 2) 

#>                 Estimate Est.Error      Q2.5     Q97.5
#> b_Intercept       522.90      9.35    504.68    540.90
#> b_email             5.11      2.15      0.93      9.34
#> b_credit_limit      0.18      0.01      0.16      0.19
#> b_risk_score     -143.52     36.64   -213.93    -71.11
#> sigma              75.33      0.74     73.90     76.80
#> Intercept         670.63      1.05    668.56    672.65
#> lprior            -10.92      0.01    -10.93    -10.91
#> lp__           -28710.05      1.56 -28713.90 -28708.00

posteriors <- m4 |> as_draws()
bayesplot::mcmc_areas(posteriors, pars = c("b_email"), prob = 0.89)

El malo feo

Ahora bien, y si nuestro modelo causal es de otra manera. Tal y como bien cuenta Carlos aquí hay un modelo mental (a priori) de como funcionan las causas. Y ese modelo es nuestra respuesta a la pregunta causal inversa, que no es otra que “¿qué causas han causado este efecto?”. A eso repondemos con por ejemplo un DAG causal. Pero ese DAG no es más que nuestra asunción, puede que no sea correcto. Otra cosa es medir el efecto de la intervención, que es lo que verás en todo lo que leas sobre inferencia causal.

Bueno, dicho esto, supongamos que nuestro modelo de como funcionan las cosas es este

email_dag_full <- dagify(
  payments ~ email + agreement + opened + credit_limit + risk_score,
  agreement ~ email + opened + credit_limit + risk_score,
  opened ~ email + credit_limit + risk_score,

  exposure = "email",
  outcome = "payments",
  labels = c(
    payments = "payments",
    email = "Envio mail",
    agreement = "agreement",
    opened = "opened",
    credit_limit = "credit_limit",
    risk_score = "risk_score"
  )
  ,
  coords = list(
    x = c(
      payments = 2,
      email = 0,
      credit_limit = 2,
      risk_score = 4,
      opened = 2,
      agreement = 1.3
    ),
    y = c(
      payments = 0,
      email = 3,
      credit_limit = 3,
      risk_score = 3,
      opened = 2,
      agreement = 1
    ))
)

curvatures = c(-0.2,-0.3, 0, 0.3, 0, 0,
               -0.2, 0, 0, 0.2, 0, 0.2)

email_dag_full |>
  tidy_dagitty() |>
  node_status() |>
  ggplot(
    aes(x, y, xend = xend, yend = yend, color = status)
  ) +
  geom_dag_edges_arc(curvature = curvatures) +
  geom_dag_point() +
  geom_dag_label(colour= "black", size = 4, alpha = 0.8) +
  scale_color_okabe_ito(na.value = "grey90") +
  theme_dag() +
  theme(legend.position = "none") +
  coord_cartesian(clip = "off")

Matheus lo explica mejor que yo, pero lo que subyace es que abrir o no un mail depende claramente de si se lo han enviado, y seguramente quien ha recibido el mail y no lo ha abierto tiene comportamiento distinto de quién si lo ha abierto en cuánto al payments. También es plausible que distintos valores en risk_score o credit_limit influyan tanto en abrir o no un email de la fintech, como en la probabilidad de intentar llegar a uno acuerdo y también en en si pagas o no.

Ahora bien, dado este DAG, podríamos pensar en que para obtener el efecto total de email no hay que condicionar por nada, porque si condicionamos por opened al ser este una variable que hace de mediator se “descuenta” ese efecto global que pasa por opened en el coeficiente de email. Entonces Matheus dice que son malos controles si condicionamos por opened, agreement , credit_limit y ris_score , (se condiciona por estos dos últimos porque si no, se abre un path no causal ).

Pero claro, yo no diría que son malos, sino que depende de qué efecto quieres estimar. Para el efecto total no hay que condicionar por ellos, pero si parea ver el efecto directo. La librería dagitty nos ayuda en eso

Para el efecto total

ggdag_adjustment_set(email_dag_full, effect = "total") + theme_dag()

Para efecto directo

ggdag_adjustment_set(email_dag_full, effect = "direct") + theme_dag()

Asi si hacemos el modelo condicionando.

m_email_ugly <- lm(payments ~ email + agreement + opened + credit_limit + risk_score,
                  data = data)

summary(m_email_ugly)

#> 
#> Call:
#> lm(formula = payments ~ email + agreement + opened + credit_limit + 
#>     risk_score, data = data)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -276.04  -50.55   -0.06   51.68  256.19 
#> 
#> Coefficients:
#>               Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   519.8975     9.4158   55.22  < 2e-16 ***
#> email           0.9214     2.6941    0.34  0.73237    
#> agreement      16.2390     4.0962    3.96  7.5e-05 ***
#> opened         -1.3044     3.7686   -0.35  0.72926    
#> credit_limit    0.1784     0.0077   23.18  < 2e-16 ***
#> risk_score   -137.0530    36.9849   -3.71  0.00021 ***
#> ---
#> Signif. codes:  
#> 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 75.2 on 4994 degrees of freedom
#> Multiple R-squared:  0.493,  Adjusted R-squared:  0.493 
#> F-statistic:  972 on 5 and 4994 DF,  p-value: <2e-16

confint(m_email_ugly, level = 0.89)

#>                  5.5 %   94.5 %
#> (Intercept)   504.8464 534.9485
#> email          -3.3851   5.2278
#> agreement       9.6914  22.7867
#> opened         -7.3284   4.7195
#> credit_limit    0.1661   0.1907
#> risk_score   -196.1726 -77.9334

Y vemos que la estimación del efecto directo de email es prácticamente 1. Y si nos vamos a como se han simulado los datos vemos que efectivamente el efecto directo es 1. Condicionar por las variables que hacen de mediadoras hacen que en el coeficiente que tenemos de email lo que se tenga sea el efecto que queda de email una vez descontado el efecto que pasa a través de las otras variables.

El bayesiano

Pero ¿y si ajustamos el dag completo a la vez?. En la inferencia causal que se suele enseñar se suele ajustar un solo modelo por pregunta causal. En este caso sería payments ~ email + credit_limit + risk_score para el efecto total y payments ~ email + agreement + opened + credit_limit + risk_score para el efecto directo.

Veamos como sería en inferencia bayesiana. Ajustamos todas las relaciones del DAG anterior

# Especificar la fórmula con múltiples partes y distribuciones
f1 <- bf(
  payments ~ email + agreement + opened + credit_limit + risk_score,
  family = gaussian() # Asumimos que payments es continua, por lo tanto usamos la distribución normal
) +
  bf(
    agreement ~ email + opened + credit_limit + risk_score,
    family = bernoulli("logit") # agreement es binaria (0 o 1), entonces usamos una distribución Bernoulli
  ) +
  bf(
    opened ~ email + credit_limit + risk_score,
    family = bernoulli("logit") # opened es binaria (0 o 1), entonces usamos también una distribución Bernoulli
  ) +
  set_rescor(FALSE) # No correlación entre las ecuaciones

# Ajustar el modelo usando brms, tarda un rato, unos 13 minutos en mi pc

m5 <- brm(
  f1,
  data = data,
  seed = 44,
  chains = 4,
  iter = 4000,
  warmup = 1000,
  cores = 4,
  file = here::here("brms_stan_models/email_flbi"),
  file_refit = "on_change"
)

round(posterior_summary(m5), digits = 2)

#>                           Estimate Est.Error      Q2.5
#> b_payments_Intercept        519.78      9.50    501.27
#> b_agreement_Intercept       -19.75     21.92    -82.29
#> b_opened_Intercept          -14.09      9.55    -39.43
#> b_payments_email              0.93      2.69     -4.33
#> b_payments_agreement         16.28      4.11      8.16
#> b_payments_opened            -1.38      3.75     -8.72
#> b_payments_credit_limit       0.18      0.01      0.16
#> b_payments_risk_score      -136.76     37.15   -210.85
#> b_agreement_email           -16.60     52.62   -142.40
#> b_agreement_opened           38.94     53.40      7.90
#> b_agreement_credit_limit      0.00      0.00      0.00
#> b_agreement_risk_score       -3.56      2.10     -7.57
#> b_opened_email               14.98      9.54      7.61
#> b_opened_credit_limit         0.00      0.00      0.00
#> b_opened_risk_score          -4.08      1.49     -6.94
#> sigma_payments               75.17      0.75     73.70
#> Intercept_payments          670.61      1.07    668.49
#> Intercept_agreement         -19.64     18.11    -72.31
#> Intercept_opened             -7.17      4.70    -19.66
#> lprior                      -22.51      2.77    -29.41
#> lp__                     -31352.41      3.69 -31360.90
#>                              Q97.5
#> b_payments_Intercept        538.50
#> b_agreement_Intercept        -5.34
#> b_opened_Intercept           -6.63
#> b_payments_email              6.22
#> b_payments_agreement         24.33
#> b_payments_opened             6.05
#> b_payments_credit_limit       0.19
#> b_payments_risk_score       -64.98
#> b_agreement_email            41.21
#> b_agreement_opened          174.15
#> b_agreement_credit_limit      0.00
#> b_agreement_risk_score        0.51
#> b_opened_email               40.38
#> b_opened_credit_limit         0.00
#> b_opened_risk_score          -1.11
#> sigma_payments               76.66
#> Intercept_payments          672.68
#> Intercept_agreement          -6.43
#> Intercept_opened             -3.54
#> lprior                      -18.69
#> lp__                     -31346.60

Y vemos que el coeficiente b_payments_email es de nuevo la estimación correcta del efecto directo. En brms cuando se ajustan varios modelos a la vez, los coeficientes se suelen especificar como b_variable_respuesta_variable_explicativa.

Pero ahora vamos a la parte divertida. Una vez ajustado este modelo, ¿puedo usarlo para estimar el efecto total de email?

Pues si, si recordamos la inferencia causal no es más que el efecto de la intervención, así que hagamos una intervención. Recorramos el DAG partiendo de que todo el mundo tiene email=1 y usemos las posterioris para ir obteniendo las variables intermedias. Y luego hagamos lo mismo con email=0. En esta parte es dónde usaremos la función rvar de la librería posterior que nos va permitir trabajar con la matriz de las muestras mcmc en la posterior de forma sencilla.

Lo primero es obtener las posteriores.

posteriores <-  as_tibble(m5)

post_rvars <- as_draws_rvars(posteriores)

Podemos ver por ejemplo la posterior del coeficiente de email en el modelo de payments. Al ser de tipo rvar nos la muestra como media y un más menos de la desviación típica, pero en ese objeto están los 12 mil valores de la posteriori

post_rvars$b_payments_email

#> rvar<12000>[1] mean ± sd:
#> [1] 0.93 ± 2.7

Vamos a ir recorriendo el DAG. Simulamos el valor de la variable opened suponiendo que email = 1 para todos los individuos y luego suponemos que email = 0 para todos. Vamos utilizando las posterioris obtenidas en el modelo. Para email = 0 no hace falta hacer nada, puesto que si no envías el email no se puede abrir, por tanto opened cuando email = 0 es 0.

p1_opened <- with(post_rvars,
                  b_opened_Intercept +
                    b_opened_email* 1 +  # intervencion email = 1
                    b_opened_credit_limit * data$credit_limit +
                    b_opened_risk_score* data$risk_score)

 
 
# p0_opened <- with(post_rvars,
#                   b_opened_Intercept +
#                     b_opened_email* 0 +  # intervencion email = 0
#                     b_opened_credit_limit * data$credit_limit +
#                     b_opened_risk_score* data$risk_score
# 
# )

En p1_opened tiene de longitud 5000, hay una variable aleatoria para cada observación de los datos

p1_opened |> length()

#> [1] 5000

Y para cada observación tengo su variable aleatorio con 12000 valores.

p1_opened[1:5]

#> rvar<12000>[5] mean ± sd:
#> [1]  0.33667 ± 0.057  -0.00085 ± 0.071 
#> [3]  0.12994 ± 0.050   0.10567 ± 0.052 
#> [5]  0.08905 ± 0.056

Pero opened es una variable dicotómica con valores 0 o 1 , por tanto voy a simular los valores de la variable. Como el modelo me lo devuelve en escala logit lo paso a probabilidades y simulo los valores

# creo una version de rbinom que funcione con el tipo de dato rvar
rvar_rbinom <-  rfun(rbinom)

opened_sim_1 <-  rvar_rbinom(nrow(data),1,inv_logit_scaled(p1_opened))



opened_sim_1[1:5]

#> rvar<12000>[5] mean ± sd:
#> [1] 0.59 ± 0.49  0.49 ± 0.50  0.54 ± 0.50  0.53 ± 0.50 
#> [5] 0.53 ± 0.50

Seguimos recorriendo el DAG, de nuevo , una parte es que todo el mundo ha recibido un email y otra que nadie lo ha recibido.

Ahora predecimos la variable agreement, para cuando email = 1 usamos lo estimado en el paso anterior para la variable opened

# predecir agreement

p1_agreement_lineal <-
  with(post_rvars,
       b_agreement_Intercept +
         b_agreement_email* 1 +  # intervencion email = 1
         b_agreement_opened * opened_sim_1  + 
         b_agreement_credit_limit * data$credit_limit +
         b_agreement_risk_score* data$risk_score

  )

#  para email = 0, opened es 0 
p0_agreement_lineal <-
  with(post_rvars,
       b_agreement_Intercept +
         b_agreement_email* 0 +  # intervencion email = 0
         b_agreement_opened * 0 + # si no hay mail no lo puede abrir
         b_agreement_credit_limit * data$credit_limit +
         b_agreement_risk_score* data$risk_score

  )

# simulamos de nuevo los valores dicotómcios

agreement_sim_1 <-  rvar_rbinom(nrow(data),1, inv_logit_scaled(p1_agreement_lineal))
agreement_sim_0 <-   rvar_rbinom(nrow(data),1, inv_logit_scaled(p0_agreement_lineal))

Vemos claramente que parea los primeros 5 individuos , cuando intervenimos diciendo que todos han recibido el email, algunos tienen hasta una prob de 0.30 en media en agreement

agreement_sim_1[1:5]

#> rvar<12000>[5] mean ± sd:
#> [1] 0.25 ± 0.43  0.30 ± 0.46  0.29 ± 0.45  0.30 ± 0.46 
#> [5] 0.29 ± 0.46

agreement_sim_0[1:5]

#> rvar<12000>[5] mean ± sd:
#> [1] 0.000000 ± 0.0000  0.000083 ± 0.0091 
#> [3] 0.000083 ± 0.0091  0.000000 ± 0.0000 
#> [5] 0.000000 ± 0.0000

Vamos a la última parte con la variable payments. Hacemos lo mismo, utilizamos los valores estimados en los pasos anteriores.

El siguiente paso es importante, puesto que puedo calcular dos cosas para cada observación. La distribución de su valor medio esperado cuando email = 1 o email = 0 y también puedo calcular la distribución de sus valores predichos.

La distribución de la media para cada individuo es mucho menos dispersa que la distribución de sus valores predichos. Es lo que se conoce como la posterior average prediction distribution en el primer caso o la posterior predictive distribution en el segundo.

Calculamos la distribución de las medias

## payments, media esperada
p1_payments <-
  with(post_rvars,
       b_payments_Intercept +
         b_payments_email* 1 +  # intervencion email = 1
         b_payments_opened * opened_sim_1  +
         b_payments_agreement * agreement_sim_1 +
         b_payments_credit_limit * data$credit_limit +
         b_payments_risk_score* data$risk_score

  )

p0_payments <-
  with(post_rvars,
       b_payments_Intercept +
         b_payments_email* 0 +  # intervencion email = 0
         b_payments_opened * 0  +
         b_payments_agreement * agreement_sim_0 +
         b_payments_credit_limit * data$credit_limit +
         b_payments_risk_score* data$risk_score

  )

y ahora la distribución de los valores predichos. Simplemente se simulan valores normales para cada observación, de forma que la media es la estimada en el paso anterior pero ponemos como desviación típica la posterior de la desviación típica de payments que nos ha dado el modelo bayesiano.

# posterior_predict hay que incorporar la varianza de payments

var_rnorm <-  rfun(rnorm)
p1_payments_pp <-
  with(post_rvars,
       var_rnorm(nrow(data), b_payments_Intercept +
         b_payments_email* 1 +  # intervencion email = 1
         b_payments_opened * opened_sim_1  +
         b_payments_agreement * agreement_sim_1 +
         b_payments_credit_limit * data$credit_limit +
         b_payments_risk_score* data$risk_score,
         sigma_payments)

  )

p0_payments_pp <-
  with(post_rvars,
       var_rnorm(nrow(data),b_payments_Intercept +
         b_payments_email* 0 +  # intervencion email = 0
         b_payments_opened * 0  +
         b_payments_agreement * agreement_sim_0 +
         b_payments_credit_limit * data$credit_limit +
         b_payments_risk_score* data$risk_score,
         sigma_payments)

  )

y aquí vemos como para el primer individuo ambas distribuciones (cuando email = 1) tienen media similar, pero diferente desviación típica.

p1_payments[1]

#> rvar<12000>[1] mean ± sd:
#> [1] 741 ± 7.3

p1_payments_pp[1]

#> rvar<12000>[1] mean ± sd:
#> [1] 742 ± 76

Pues con esto ya nos podemos calcular para cada observación el CATE (conditional average treatment effect), pero con la ventaja de que podemos hacerlo sobre la distribución de las medias o sobre la distribución de los valores predichos.

cate_ind_posterior_predict <-  p1_payments_pp - p0_payments_pp
cate_ind_posterior_average_predict <-  p1_payments - p0_payments

data_with_cate <-  data.frame(data, cate_ind_posterior_average_predict, cate_ind_posterior_predict) |>
  rownames_to_column(var = "id")



head(data_with_cate)

#>   id payments email opened agreement credit_limit
#> 1  1      820     1      1         0       1652.1
#> 2  2      730     1      0         0        954.3
#> 3  3      650     1      1         1       1226.8
#> 4  4      530     1      0         0       1160.9
#> 5  5      510     1      1         0       1149.8
#> 6  6      760     1      0         0        720.8
#>   risk_score cate_ind_posterior_average_predict
#> 1     0.5659                          4.2 ± 7.4
#> 2     0.4668                          5.1 ± 7.6
#> 3     0.5058                          4.9 ± 7.6
#> 4     0.4946                          5.0 ± 7.6
#> 5     0.4957                          5.0 ± 7.6
#> 6     0.3634                          6.7 ± 8.0
#>   cate_ind_posterior_predict
#> 1                  4.6 ± 107
#> 2                  6.9 ± 107
#> 3                  6.5 ± 107
#> 4                  4.8 ± 107
#> 5                  5.0 ± 106
#> 6                  5.7 ± 106

Y ahora vamos a pintar un poco. Vemos que la distribución a nivel individual de los valores predichos del cate es muy dispersa.

sim <-  data_with_cate |> sample_n(10)

sim |>
  ggplot() +
  ggdist::stat_halfeye(aes( y = id, xdist = cate_ind_posterior_predict), fill ="lightblue")

y podríamos ver la distribución del cate en la muestra. Podemos tomar medias de los cate individuales. Y vemos que de esta forma podemos recuperar el efecto “total”.

data_with_cate |>
  summarise(
    cate_medio = rvar_mean(cate_ind_posterior_predict)) |>
  ggplot() +
  ggdist::stat_halfeye(aes( xdist = cate_medio), fill ="lightblue") +
    labs (title= "Efecto total usando modelo bayesiano")

Nota final

En inferencia causal podemos hacer un modelo para cada pregunta causal, o si usamos esta forma de hacerlo en bayesiano, podemos tener un solo modelo global y utilizando las posterioris responder a varias preguntas causales.

Adenda

Al tener las posterioris, podemos de forma sencilla tener distribuciones de los efectos condicionando por las variables que queramos. Por ejemplo podemos ver como cambia el efecto según combinaciones de credit_limit y risk_score

data_with_cate |>
  sample_n(300) |>
  mutate(
    credit_limit_cut = cut_number(credit_limit, 3),
    risk_cut  = cut_number(risk_score, 3)
    )|>
  group_by(credit_limit_cut, risk_cut) |>
  summarise(
    cate_medio = rvar_mean(cate_ind_posterior_average_predict)
  ) |>
  ggplot() +
  ggdist::stat_halfeye(aes(y = credit_limit_cut, xdist = cate_medio, fill = risk_cut), alpha = 0.4)