ch06_scm.qmd

# Structural Causal Models {#SCM}

```{r }
#| include: false
library(dplyr)
library(tidyr)
library(tidybayes)
library(rethinking)
library(brms)
library(skimr)
library(ggplot2)
library(dagitty)
library(ggdag)
library(ggdist)
library(patchwork)
library(paletteer)
```

Some options to facilitate the computations

```{r}
#  For execution on a local, multicore CPU with excess RAM
options(mc.cores = parallel::detectCores())
#  To avoid recompilation of unchanged Stan programs
rstan::rstan_options(auto_write = TRUE)
```

The default theme used by `ggplot2`

```{r}
# The default theme used by ggplot2
ggplot2::theme_set(ggthemes::theme_fivethirtyeight())
ggplot2::theme_update(title = element_text(color = "midnightblue"))
```

This is a ver y important point, in the intro to chapter 6.

> Regression will not sort it out. Regression is indeed an oracle, but a cruel one. It speaks in riddle and delights in punishing us for asking bad questions. The selection-distortion effect can happen inside of a multiple regression, becase the fact of adding a predictor induces statistical selection within the model, a phenomenon that goes by the unhelpful name of **collider bias**. This can mislead us into believing, for axample, that there is a negative associaiton between newswothiness and trustworthiness in general, ehen in fact it is just a consequence of conditioning on some variable.

## Multicollinearity

Multicollinearity means a very strong association between 2 or more predictor variables.

### Multicollinear legs

Create the data

```{r}
set.seed(6)
dataLegs <- 
  tibble(height   = rnorm(100, mean = 10, sd = 2),
         leg_prop = runif(100, min = 0.4, max = 0.5)) |> 
  mutate(leg_left  = leg_prop * height + rnorm(100, mean = 0, sd = 0.02),
         leg_right = leg_prop * height + rnorm(100, mean = 0, sd = 0.02))
```

which has the following correlations

```{r}
dataLegs |>
   cor() |>
   ggcorrplot::ggcorrplot(type = "lower", lab = TRUE, digits = 2,
                          show.legend = FALSE,
                          title = "Correlations between variables") +
  scale_fill_paletteer_c("pals::warmcool", direction = -1)
```

```{r ch06_fit06_01}
tictoc::tic(msg = sprintf("run time of %s, use the cache.", "110 secs."))
fit06_01 <- xfun::cache_rds({brm(data = dataLegs, 
      family = gaussian,
      height ~ 1 + leg_left + leg_right,
      prior = c(prior(normal(10, 100), class = Intercept),
                prior(normal(2, 10), class = b),
                prior(exponential(1), class = sigma)),
      iter = 2000, warmup = 1000, chains = 4, cores = detectCores(),
      seed = 6)},
  file = "ch06_fit06_01")
tictoc::toc()
```

```{r}
posterior_summary(fit06_01) |> round(digits = 3)
```

```{r}
tidybayes::get_variables(fit06_01)
```

```{r}
tidybayes::gather_draws(fit06_01, b_Intercept, b_leg_left, b_leg_right, sigma) |>
  mean_hdi() |>
  ggplot(aes(x = .value, xmin = .lower, xmax = .upper, y = .variable, color = .variable)) +
  geom_pointinterval() +
  ggrepel::geom_text_repel(aes(label = round(.value, 2))) +
  scale_color_paletteer_d("Manu::Kereru") +
  labs(title = "Leg model",
       x = "value", y = NULL, color = NULL)
```

### Multicollinear milk

```{r}
data(milk)
dataMilk <- milk
dataMilk <- dataMilk |>
  mutate(K = as.vector(scale(kcal.per.g)),
         `F` = as.vector(scale(perc.fat)),
         L = as.vector(scale(perc.lactose)))

```

```{r}
plotMilk <- list()
plotMilk <- within(plotMilk, {
  fun_diag <- function(data, mapping, ...){
    ggplot(data = data, mapping = mapping) +
      geom_density(linewidth = 1, color = "pink")}
  fun_upper <- function(data, mapping) {
    ggplot(data = data, mapping = mapping) +
      geom_text(size = 1, color = "blue")}
  fun_lower <- function(data, mapping) {
    ggplot(data = data, mapping = mapping) +
      stat_density2d(linewidth = 1/3, color = "blue") +
      geom_point(size = 1, color = "blue") +
      geom_smooth(method = "loess")}
  plot <- GGally::ggpairs(
    data = dataMilk,
    columns = c("K", "F", "L"),
    title = "Milk example",
    diag = list(continuous = fun_diag),
    lower = list(continuous = fun_lower))
})
plotMilk$plot
```

```{r ch06_fit06_03}
tictoc::tic(msg = sprintf("run time of %s, use the cache.", "80 secs."))
fit06_03 <- xfun::cache_rds({
  # k regressed on f
  brm(data = dataMilk, 
      family = gaussian,
      K ~ 1 + `F`,
      prior = c(prior(normal(0, 0.2), class = Intercept),
                prior(normal(0, 0.5), class = b),
                prior(exponential(1), class = sigma)),
      iter = 2000, warmup = 1000, chains = 4, cores = detectCores(),
      seed = 6)},
  file = "ch06_fit06_03")
tictoc::toc()
```

```{r ch06_fit06_04}
tictoc::tic(msg = sprintf("run time of %s, use the cache.", "90 secs."))
fit06_04 <- xfun::cache_rds({
  # k regressed on f
  brm(data = dataMilk, 
      family = gaussian,
      K ~ 1 + L,
      prior = c(prior(normal(0, 0.2), class = Intercept),
                prior(normal(0, 0.5), class = b),
                prior(exponential(1), class = sigma)),
      iter = 2000, warmup = 1000, chains = 4, cores = detectCores(),
      seed = 6)},
  file = "ch06_fit06_04")
tictoc::toc()
```

and the coefficients are

```{r}
posterior_summary(fit06_03) |> round(digits = 3)
```

```{r}
posterior_summary(fit06_04) |> round(digits = 2)
```

and the multivariate which shows that each variable has now a much larger variance caused by the colinearity.

```{r ch06_fit06_05}
tictoc::tic(msg = sprintf("run time of %s, use the cache.", "60 secs."))
fit06_05 <- xfun::cache_rds({
  brm(data = dataMilk, 
      family = gaussian,
      K ~ 1 + `F` + L,
      prior = c(prior(normal(0, 0.2), class = Intercept),
                prior(normal(0, 0.5), class = b),
                prior(exponential(1), class = sigma)),
      iter = 2000, warmup = 1000, chains = 4, cores = detectCores(),
      seed = 6)},
  file = "ch06_fit06_05")
tictoc::toc()
```

## Post-treatment bias

```{r}
# how many plants would you like?
set.seed(7)
dataPlants <- tibble(
  h0        = rnorm(100, mean = 10, sd = 2), 
  treatment = rep(0:1, each = 100 / 2),
  fungus    = rbinom(100, size = 1, prob = .5 - treatment * 0.4),
  h1        = h0 + rnorm(100, mean = 5 - 3 * fungus, sd = 1))
skimr::skim(dataPlants)
```

### A prior is born

> If we center our prior for $p$ on 1, that implies an expectation of no change in height. That is less than we know. But we would allow $p$ to be less than 1, in case the experiment ges wrong. We also want to ensure $p>0$.

Therefore we use $p$ with log-normal distribution

$$
\begin{align*}
h_{1,i} &\sim \mathcal{N}(\mu_i, \sigma) \\
\mu_i &= h_{0,i} \times p \\
p &\sim \mathcal{LogNormal}(0, 0.25) \\
\sigma &\sim \mathcal{Exp}(1)
\end{align*}
$$

```{r ch06_fit06_06}
tictoc::tic(msg = sprintf("run time of %s, use the cache.", "90 secs."))
fit06_06 <- xfun::cache_rds({
  out <- brm(
    data = dataPlants, 
    family = gaussian,
    h1 ~ 0 + h0,
    prior = c(prior(lognormal(0, 0.25), class = b, lb = 0),
              prior(exponential(1), class = sigma)),
    iter = 2000, warmup = 1000, chains = 4, cores = detectCores(),
    seed = 6)
  out <- brms::add_criterion(out, criterion = c("waic", "loo"))
  out},
  file = "ch06_fit06_06")
tictoc::toc()
```

```{r}
brms::posterior_summary(fit06_06) |> round(digits = 2)
```

So the increase is 1.38 relative to $h_0$.

Now including the treatment and fungus we have

$$
\begin{align*}
h_{1,i} &\sim \mathcal{N}(\mu_i, \sigma) \\
\mu_i &= h_{0,i} \times p \\
p &\sim \alpha + \beta_1 treatment_i + \beta_2 fungus_i \\
\alpha &\sim \mathcal{LogNormal}(0, 0.25) \\
\beta_1 &\sim \mathcal{N}(0, 0.5) \\
\beta_2 &\sim \mathcal{N}(0, 0.5) \\
\sigma &\sim \mathcal{Exp}(1)
\end{align*}
$$

```{r ch06_fit06_07}
tictoc::tic(msg = sprintf("run time of %s, use the cache.", "80 secs."))
fit06_07 <- xfun::cache_rds({
  out <- brm(
    data = dataPlants, 
    family = gaussian,
    bf(h1 ~ h0 * (a + t * treatment + f * fungus),
       a + t + f ~ 1, 
       nl = TRUE),
    prior = c(prior(lognormal(0, 0.2), nlpar = a, lb = 0),
                prior(normal(0, 0.5), nlpar = t),
                prior(normal(0, 0.5), nlpar = f),
                prior(exponential(1), class = sigma)),
    iter = 2000, warmup = 1000, chains = 4, cores = detectCores(),
    seed = 6)
  out <- brms::add_criterion(out, criterion = c("waic", "loo"))
  out},
  file = "ch06_fit06_07")
tictoc::toc()
```

```{r}
brms::posterior_summary(fit06_07) |> round(digits = 2)
```

Now the effect of the treatment is almost non existent.

### Blocked by consequence

The problem is that the fungus is part of a chain between the treatment and the growth.

```{r}
ggdag::dagify(h1 ~ h0, h1 ~ `F`, `F` ~ `T`) |>
    ggdag::ggdag(layout = "sugiyama", node_size = 8, text_col = "yellow") +
    ggdag::theme_dag_blank(
      panel.background = element_rect(fill = "snow2", color = "snow2"))
```

so now we redo the model but without the fungus effect.

$$
\begin{align*}
h_{1,i} &\sim \mathcal{N}(\mu_i, \sigma) \\
\mu_i &= h_{0,i} \times p \\
p &\sim \alpha + \beta_1 treatment_i \\
\alpha &\sim \mathcal{LogNormal}(0, 0.25) \\
\beta_1 &\sim \mathcal{N}(0, 0.5) \\
\sigma &\sim \mathcal{Exp}(1)
\end{align*}
$$

```{r ch06_fit06_08}
tictoc::tic(msg = sprintf("run time of %s, use the cache.", "80 secs."))
fit06_08 <- xfun::cache_rds({
  out <- brm(
    data = dataPlants, 
    family = gaussian,
    bf(h1 ~ h0 * (a + t * treatment),
      a + t ~ 1, nl = TRUE),
    prior = c(prior(lognormal(0, 0.2), nlpar = a, lb = 0),
                prior(normal(0, 0.5), nlpar = t),
                prior(exponential(1), class = sigma)),
    iter = 2000, warmup = 1000, chains = 4, cores = detectCores(),
    seed = 6)
  out <- brms::add_criterion(out, criterion = c("waic", "loo"))},
  file = "ch06_fit06_08")
tictoc::toc()
```

```{r}
posterior_summary(fit06_08) |> round(digits = 3)
```

and we now see more treatment effect.

## Collider bias

### Collider of false sorrow

```{r}
dataHappy <- rethinking::sim_happiness(seed = 1977, N_years = 1000)
# select age > 17 and rescale to [0, 1] and create indexed factor
# creating factor makes it easer with brms
dataHappy_gt17 <- dataHappy |>
  filter(age > 17) |>
  mutate(A = scales::rescale(age, to = c(0, 1)),
         mid = factor(married + 1, labels = c("single", "married")))
# glimpse(dataHappy_gt17)
```

```{r ch06_fit06_09}
tictoc::tic(msg = sprintf("run time of %s, use the cache.", "80 secs."))
fit06_09 <- xfun::cache_rds({
  out <- brm(
    data = dataHappy_gt17, 
    family = gaussian,
    happiness ~ 0 + mid + A,
    prior = c(prior(normal(0, 1), class = b, coef = midmarried),
                prior(normal(0, 1), class = b, coef = midsingle),
                prior(normal(0, 2), class = b, coef = A),
                prior(exponential(1), class = sigma)),
    iter = 2000, warmup = 1000, chains = 4, cores = detectCores(),
    seed = 6)
  out <- brms::add_criterion(out, criterion = c("waic", "loo"))},
  file = "ch06_fit06_09")
tictoc::toc()
```

```{r}
posterior_summary(fit06_09)
```

The fit finds that the effect of age on happiness is negative

now lets do it without the marriage factor

```{r ch06_fit06_10}
tictoc::tic(msg = sprintf("run time of %s, use the cache.", "90 secs."))
fit06_10 <- xfun::cache_rds({
  out <- brm(
    data = dataHappy_gt17,
    family = gaussian,
    happiness ~ 1 + A,
    prior = c(prior(normal(0, 1), class = Intercept),
                prior(normal(0, 2), class = b),
                prior(exponential(1), class = sigma)),
    iter = 2000, warmup = 1000, chains = 4, cores = detectCores(),
    seed = 6)
  out <- brms::add_criterion(out, criterion = c("waic", "loo"))},
  file = "ch06_fit06_10")
tictoc::toc()
```

```{r}
posterior_summary(fit06_10)
```

Now the age has no effect on happiness! When we include marriage, we include a spurious association.

### The haunted DAG

## Confronting counfounding

See Overthinking box in section 6.4.3. Confounding occurs when

$$
Pr(Y \mid X) \neq Pr(Y \mid do(X))
$$

## Summary