AE 20: Causal Inference from a Randomized Experiment

Packages

library(tidyverse)
library(tidymodels)

Data

set.seed(42)
n <- 46
wheat_trial <- tibble(
  id                 = 1:n,
  age                = round(runif(n, 4, 22), 1),
  baseline_wheat_ige = round(rexp(n, rate = 1/30), 1),
  vwg_oit            = c(rep(1, 23), rep(0, 23))
) |>
  mutate(
    max_dose_mg = pmax(0, round(
      150 + 950 * vwg_oit - 3 * baseline_wheat_ige + 12 * age + rnorm(n, sd = 90)
    ))
  )

glimpse(wheat_trial)

Rows: 46
Columns: 5
$ id                 <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, …
$ age                <dbl> 20.5, 20.9, 9.2, 18.9, 15.6, 13.3, 17.3, 6.4, 15.8,…
$ baseline_wheat_ige <dbl> 13.3, 7.1, 30.8, 32.4, 38.6, 17.1, 90.6, 14.9, 10.6…
$ vwg_oit            <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ max_dose_mg        <dbl> 1331, 1391, 1126, 960, 1197, 1175, 1052, 1184, 1384…

Variables:

Variable	Description
`id`	Unique participant identifier
`age`	Participant age (years)
`baseline_wheat_ige`	Baseline wheat IgE level (kU/L) — a measure of allergy severity
`vwg_oit`	Treatment indicator: 1 = VWG OIT, 0 = Placebo
`max_dose_mg`	Maximum tolerated wheat protein dose (mg) at the year-1 food challenge

Exercise 1: Study design

With your teammate, fill in the following based on the paper and today’s lecture:

Treatment indicator (\(Z\)): vwg_oit — 1 if the participant received VWG OIT, 0 if placebo
Treatment group (\(Z = 1\)): Participants receiving vital wheat gluten (VWG) oral immunotherapy
Control group (\(Z = 0\)): Participants receiving a matched placebo
Outcome (\(Y\)): max_dose_mg — maximum tolerated wheat protein dose (mg) at the year-1 food challenge
Covariates (\(X\)): age, baseline_wheat_ige (and in the real study: sex, threshold dose at baseline)

Exercise 2: Exploratory data analysis

a. Create a visualization comparing the distribution of max_dose_mg across the two treatment groups. What do you notice?

ggplot(
  wheat_trial,
  aes(
    x = factor(vwg_oit, labels = c("Placebo", "VWG OIT")),
    y = max_dose_mg
  )
) +
  geom_boxplot() +
  labs(
    x     = "Treatment group",
    y     = "Max tolerated dose (mg)",
    title = "Max wheat dose by treatment group"
  )

The VWG OIT group has a substantially higher median max tolerated dose than the placebo group, and the distributions show little overlap. This suggests a strong treatment effect.

b. Compute the mean and standard deviation of max_dose_mg for each treatment group.

wheat_trial |>
  group_by(vwg_oit) |>
  summarise(
    n         = n(),
    mean_dose = mean(max_dose_mg),
    sd_dose   = sd(max_dose_mg)
  )

# A tibble: 2 × 4
  vwg_oit     n mean_dose sd_dose
    <dbl> <int>     <dbl>   <dbl>
1       0    23      234.    156.
2       1    23     1201.    144.

Exercise 3: Modeling

a. Fit a linear regression model predicting max_dose_mg from vwg_oit.

wheat_fit <- linear_reg() |>
  fit(max_dose_mg ~ vwg_oit, data = wheat_trial)

tidy(wheat_fit)

# A tibble: 2 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)     234.      31.3      7.48 2.29e- 9
2 vwg_oit         967.      44.2     21.9  2.79e-25

b. Interpret the slope coefficient for vwg_oit in context. Since this is a randomized experiment, use causal language in your interpretation (e.g., “causes,” “leads to”).

Receiving VWG OIT causes an estimated increase of approximately 950 mg in max tolerated wheat protein dose compared to placebo, on average. Because treatment was randomly assigned, we can attribute this difference to the treatment itself rather than to any confounding variable.

Exercise 4: Inference

Option A: Bootstrap confidence interval

a. Calculate the observed fit using specify() |> fit().

observed_fit <- wheat_trial |>
  specify(max_dose_mg ~ vwg_oit) |>
  fit()

observed_fit

# A tibble: 2 × 2
  term      estimate
  <chr>        <dbl>
1 intercept     234.
2 vwg_oit       967.

b. Generate 1000 bootstrap samples and fit a model to each one.

set.seed(1120)

boot_fits <- wheat_trial |>
  specify(max_dose_mg ~ vwg_oit) |>
  generate(reps = 1000, type = "bootstrap") |>
  fit()

boot_fits

# A tibble: 2,000 × 3
# Groups:   replicate [1,000]
   replicate term      estimate
       <int> <chr>        <dbl>
 1         1 intercept     225.
 2         1 vwg_oit       942.
 3         2 intercept     246.
 4         2 vwg_oit       932.
 5         3 intercept     269.
 6         3 vwg_oit       943.
 7         4 intercept     201.
 8         4 vwg_oit      1021.
 9         5 intercept     277.
10         5 vwg_oit       955.
# ℹ 1,990 more rows

c. Compute the 95% confidence interval using the percentile method.

get_confidence_interval(
  boot_fits,
  point_estimate = observed_fit,
  level          = 0.95,
  type           = "percentile"
)

# A tibble: 2 × 3
  term      lower_ci upper_ci
  <chr>        <dbl>    <dbl>
1 intercept     173.     299.
2 vwg_oit       877.    1053.

d. Interpret the confidence interval for vwg_oit in context. Use causal language.

We are 95% confident that receiving VWG OIT causes an increase of between approximately [lower bound] mg and [upper bound] mg in max tolerated wheat protein dose compared to placebo, on average. Because the interval does not include 0, this is consistent with a statistically significant treatment effect.

Option B: Hypothesis test

a. State the null and alternative hypotheses.

\(H_0\): The slope predicting max tolerated dose from VWG OIT treatment is 0 — VWG OIT has no effect on max tolerated dose, \(\beta_1 = 0\).

\(H_A\): The slope is not 0 — VWG OIT has a significant effect on max tolerated dose, \(\beta_1 \neq 0\).

b. Calculate the observed fit and simulate the null distribution using 1000 permutations.

observed_fit <- wheat_trial |>
  specify(max_dose_mg ~ vwg_oit) |>
  fit()

set.seed(24601)
null_dist <- wheat_trial |>
  specify(max_dose_mg ~ vwg_oit) |>
  hypothesize(null = "independence") |>
  generate(reps = 1000, type = "permute") |>
  fit()

c. Visualize the null distribution and shade the p-value region.

visualize(null_dist) +
  shade_p_value(obs_stat = observed_fit, direction = "two-sided")

d. Calculate the p-value.

null_dist |>
  get_p_value(obs_stat = observed_fit, direction = "two-sided")

# A tibble: 2 × 2
  term      p_value
  <chr>       <dbl>
1 intercept       0
2 vwg_oit         0

e. State your conclusion in context.

Since \(p < 0.05\), we reject \(H_0\). We have sufficient evidence of a statistically significant treatment effect: VWG OIT causes a significant increase in maximum tolerated wheat protein dose compared to placebo.

Exercise 5: Further covariate adjustment

a. Fit a model that adds baseline_wheat_ige as an additional predictor. How does the vwg_oit estimate compare to the model in Exercise 3?

wheat_fit_adj <- linear_reg() |>
  fit(max_dose_mg ~ vwg_oit + baseline_wheat_ige, data = wheat_trial)

tidy(wheat_fit_adj)

# A tibble: 3 × 5
  term               estimate std.error statistic  p.value
  <chr>                 <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)          312.      27.8       11.2  2.22e-14
2 vwg_oit              947.      34.2       27.7  4.57e-29
3 baseline_wheat_ige    -1.84     0.328     -5.62 1.31e- 6

The vwg_oit estimate is very similar to the unadjusted model. This is expected: because treatment was randomly assigned, the two groups should already be balanced on baseline_wheat_ige, so including it should not substantially change the treatment effect estimate.

b. Why might we include baseline_wheat_ige in our model even if the two groups are already balanced on this variable?

Even when groups are balanced on a covariate, including it in the model can reduce residual variance and produce more precise (narrower confidence interval) estimates of the treatment effect. In a small trial like this (n = 46), such precision gains can be meaningful.