13.1 Standardization as an alternative to IP weighting

• Average treatment effect (ATE), or average causal effect (ACE) of smoking cessation on weight gain
• Causal contrast to assess: $E\left[Y^{a=1,c=0}\right]-E\left[Y^{a=0,c=0}\right]$
• Marginal weight gain difference had everyone stopped smoking vs had no one stopped smoking (everyone was treated vs everyone was untreated)
• Ignoring loss-to-follow-up from 1971 (baseline) through 1982
• Ignoring time-varying nature of the treatment and potential time-varying confounding
• Assuming conditional exchangeability based on age, sex, race, education, smoking intensity and duration, physical activity, exercise, and baseline weight

13.1 Standardization as an alternative to IP weighting

• Same causal question
• G-formula for ATE estimation
• NHEFS data from 1629 cigarette smokers aged 25-74 years who had a baseline visit and a follow-up visit
• Assumptions recap in the presentation for Chapter 12

13.1 Standardization as an alternative to IP weighting

• $\hat{E}[Y|A=1]$ - $\hat{E}[Y|A=0]$
• Does not have a causal interpretation because exchangeability does not hold
  • Quitters and non-quitters are different in respect to characteristics affecting weight gain

lm(data = data, formula = wt82_71 ~ qsmk) %>% broom::tidy(., conf.int = T) %>% select(term, conf.low, estimate, conf.high)

## # A tibble: 2 x 4
##   term        conf.low estimate conf.high
##   <chr>          <dbl>    <dbl>     <dbl>
## 1 (Intercept)     1.54     1.98      2.43
## 2 qsmk            1.66     2.54      3.43

13.1 Standardization as an alternative to IP weighting

Assume

• Exchangeability: $Y^a\perp \perp A$
  • Treated individuals (A=1) had they been untreated (A=0) had the same probability of the potential outcome (PO)
    • Conditional exchangeability: $Y^a\perp \perp A|L$, where L closes all back-door paths between A and Y
    • When accounted for all variables in L, treated individuals (A=1) had they been untreated (A=0) had the same probability of the potential outcome
• Positivity
  • Positive probability of observing each level of treatment in each strata of L
  • $Pr(A=a|L=l>0)$ for all $a$ and $l$
• Consistency
  • We observe PO - the one under actually received treatment
  • $Pr[Y^{a=1}|A=1]=Pr[Y=1|A=1]$
  • Well-defined intervention ➡️ well-defined (de)confounding
• No model misspesification

13.1 Standardization as an alternative to IP weighting

• Adjustment (standardization) formula
  • ${\Sigma }_{l}E[Y|A=a,L=l]Pr[L=l]$
  • if conditional exchangeability holds given L
  • in this example we also condition on C=0: ${\Sigma }_{l}E[Y|A=a,C=0,L=l]Pr[L=l]$
• To compute the standardized mean outcome in the uncensored treated, we need to compute the mean outcomes in the uncensored treated in each stratum $l$ of the confounders $L$
  • $E[Y|A=1,C=0,L=l]$

13.2 Estimating the mean outcome via modeling

• High-dimensional data with few observation ➡️ nonparametric modeling is not going to work
• Parametric model:
  • $\hat{E}[Y|A=a,C=0,L=l]Pr[L=l]$

13.3 Standardizing the mean outcome to the confounder distribution

n <- nrow(data)
fit_qsmk <- glm(data = data, qsmk ~ 1, family = binomial())
p_qsmk <- predict(fit_qsmk, type = "response")
p_qsmk <- p_qsmk[[1]]
# fit a regression of the outcome on the exposure and confounders
Q <- glm(data = data, formula = wt82_71 ~ qsmk + sex + race + age + I(age*age) + education + smokeintensity + I(smokeintensity*smokeintensity) + smokeyrs + I(smokeyrs*smokeyrs) + exercise + active + wt71 + I(wt71*wt71), family = gaussian())
# create predicted Y_a for all observations
data %<>% mutate(
  new_qsmk = rbinom(n = n, size = 1, prob = p_qsmk)
)
data %<>% mutate(
  wt82_71_pred = predict(Q, type = "response")
)
# generate data sets where a is set to 0 and 1
data_qsmk0 <- data %>% mutate(qsmk = 0) 
data_qsmk1 <- data %>% mutate(qsmk = 1)

# compute PO for all observations, sets qsmk=0
data_qsmk0 %<>% mutate(
  wt82_71_po = predict(Q, newdata = data_qsmk0, type = "response")
)
# compute PO for all observations, sets qsmk=1
data_qsmk1 %<>% mutate(
  wt82_71_po = predict(Q, newdata = data_qsmk1, type = "response")
)
data_new <- bind_rows(data_qsmk0, data_qsmk1)
glm(wt82_71_po ~ qsmk, family = gaussian(link = "identity"), data = data_new) %>% broom::tidy()

## # A tibble: 2 x 5
##   term        estimate std.error statistic   p.value
##   <chr>          <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)     1.75    0.0718      24.3 6.89e-120
## 2 qsmk            3.46    0.102       34.1 6.41e-217

library(tidymodels)
set.seed(972188635)
boots <- bootstraps(data, times = 1e3, apparent = FALSE)
gform <- function(data) {
  glm(formula = wt82_71~qsmk + sex + race + age + I(age*age) + education + smokeintensity + I(smokeintensity*smokeintensity) + smokeyrs + I(smokeyrs*smokeyrs) + exercise + active + wt71 + I(wt71*wt71), family = gaussian(), data = data)
}
p_qsmk <- glm(data = data, qsmk~1, family = binomial())
qsmk_pred <- predict(p_qsmk, type = "response")
n <- nrow(data)
qsmk_pred <- qsmk_pred[[1]]
boot_models <- boots %>% 
  mutate(
    model = map(.x = splits, ~gform(data = .x)),
    splits = map(splits, ~ as_tibble(.x)),
    # simulate qsmk with observed frequency as if it was randomized
    splits = map(splits, ~mutate(.x, qsmk = rbinom(n = n, size = 1, prob = qsmk_pred))),
    # simulate P.O.s
    splits = map(splits, ~mutate(.x, po = rnorm(n = n, mean = predict(Q, newdata = .x), sd = predict(Q, newdata = .x, se.fit = T)$se.fit))),
    gform = map(splits, ~glm(po ~ qsmk, family = gaussian(link = "identity"), data = .x)),
    gform_tidy = map(gform, ~broom::tidy(.x))
    )

# percentile method for 95% CI
pct_res <- boot_models %>% 
  select(id, gform_tidy) %>% 
  bind_rows(.id = "boot") %>% 
  unnest(cols = c(gform_tidy)) %>% 
  filter(term == "qsmk")

13.3 Standardizing the mean outcome to the confounder distribution

p <- pct_res %>% ggplot(aes(x = estimate)) +
  geom_histogram(bins = 100) +
  theme_light() +
  theme_minimal() +
  geom_vline(aes(xintercept = quantile(estimate, probs = 0.5)), color = "cyan") +
  geom_vline(aes(xintercept = quantile(estimate, probs = 0.025)), color = "lightblue") +
  geom_vline(aes(xintercept = quantile(estimate, probs = 0.975)), color = "lightblue")

pct_res %>% summarise(
    Q2.5 = quantile(estimate, probs = 0.025), Q50 = quantile(estimate, probs = 0.5), Q95.7 = quantile(estimate, probs = 0.975)
  )

## # A tibble: 1 x 3
##    Q2.5   Q50 Q95.7
##   <dbl> <dbl> <dbl>
## 1  3.14  3.47  3.81

# M. Hernan: 3.5 (2.6-4.5)

13.4 IP weighting or standardization?

• Better precision
• More flexible
• Allows for detailed mediation and interaction analyses
• Both IPTW and coputation of g-formula can be misspecidied
• Doubly robust methods

13.5 How seriously do we take our estimates?

• Identifiability conditions
• Measurement bias
• Selection bias
• Model misspecifications
• Keep calm and remain skeptical 🤓

References

Hernán MA, Robins JM (2020). Causal Inference: What If. Boca Raton: Chapman & Hall/CRC (v. 30mar21)