Exercise 2 (20 minutes)

You practice: becoming the simulator — setting a confounder's strength and sign yourself, predicting the bias from the formula, and checking whether OLS agrees. When you're done, assist some of the others …

In exercise 1 you adjusted a real estimate — but you could never check whether the adjusted number was right, because nobody knows the true effect of the Bali attack. Here you do know it, because you set it yourself.

1. Copy in this simulation. It is the lecture's toy example, all three steps in one chunk. Three numbers in it are the knobs you will turn: the true effect (1 * bali), how strongly the confounder affects the outcome (0.3 * age, \(\beta_{C \rightarrow Y}\)), and how strongly it drives selection into treatment (-0.1 * age, \(\beta_{D \rightarrow C}\)).

pacman::p_load(tidyverse, estimatr, modelsummary)

set.seed(831983) # So everyone gets the SAME "random" numbers
toydat <- tibble(
  age = rnorm(n = 500, mean = 40, sd = 5)) %>%     # Step 1: age
  mutate(
    # Step 2: who gets treated (older people are interviewed earlier)
    int_day = (-0.1 * age) + rnorm(n = 500, mean = 0, sd = 1),
    bali = case_when(
      int_day < mean(int_day) ~ 0,   # interviewed before the attack
      int_day >= mean(int_day) ~ 1), # interviewed after the attack
    # Step 3: the outcome. THE TRUTH: the Bali effect is 1
    xeno = (1 * bali) + (0.3 * age) + rnorm(n = 500, mean = 0, sd = 1))

Stuck?

• Did pacman::p_load(...) run without red errors? Fix that first.
• Getting different numbers from the ones below? Check you kept set.seed(831983) and re-ran the whole chunk, from set.seed() down — never just a part of it — the seed is what makes everyone's "random" data identical.
• Odd results? Session → Restart R, then rerun from the top.
• Still stuck? Ask in the t-R-ouble forum on Absalon.

2. Verify the formula. Fit the three regressions you need, and confirm that \(\tilde\beta = \beta + (\beta_{C \rightarrow Y} \times \beta_{D \rightarrow C})\) actually holds.

Hint

You need three models: xeno ~ bali (biased), xeno ~ bali + age (adjusted), and — the one students forget — age ~ bali, which measures the imbalance.

R solution — don't peek too early ;) !

ols_bi   <- lm_robust(xeno ~ bali, data = toydat)        # biased estimate
ols_mult <- lm_robust(xeno ~ bali + age, data = toydat)  # adjusted estimate
ols_dc   <- lm_robust(age ~ bali, data = toydat)         # the imbalance

b_bi <- coef(ols_bi)["bali"]
b_mu <- coef(ols_mult)["bali"]
b_cy <- coef(ols_mult)["age"]
b_dc <- coef(ols_dc)["bali"]

# Does the decomposition close?
c(bivariate = b_bi, adjusted_plus_OVB = b_mu + (b_cy * b_dc))
#         bivariate.bali adjusted_plus_OVB.bali 
#                   0.17                   0.17

The bivariate estimate is , the adjusted estimate is , and the omitted variable bias is . Given that you set the true effect to 1, which estimate got closer? The adjusted one — the bivariate estimate is badly biased downwardsThe bivariate oneBoth are equally closeNeither is close

3. Predict before you compute. Now make the confounder three times stronger on the outcome — change 0.3 * age to 0.9 * age — but change nothing else. What happens to the bivariate estimate? It gets closer to the truth of 1It stays about the sameIt becomes even more biased — so much that it turns negativeIt doubles to about 2

R solution — check your prediction

set.seed(831983) # So everyone gets the SAME "random" numbers
toy_strong <- tibble(
  age = rnorm(n = 500, mean = 40, sd = 5)) %>%     # Step 1: age
  mutate(
    # Step 2: who gets treated (older people are interviewed earlier)
    int_day = (-0.1 * age) + rnorm(n = 500, mean = 0, sd = 1),
    bali = case_when(
      int_day < mean(int_day) ~ 0,   # interviewed before the attack
      int_day >= mean(int_day) ~ 1), # interviewed after the attack
    # Step 3: the outcome. THE TRUTH: the Bali effect is 1
    xeno = (1 * bali) + (0.9 * age) + rnorm(n = 500, mean = 0, sd = 1)) #<< 0.9 instead of 0.3

lm_robust(xeno ~ bali, data = toy_strong)        # bivariate
lm_robust(xeno ~ bali + age, data = toy_strong)  # adjusted

The bivariate estimate is now about -1.67 — it has the wrong sign. A strong enough confounder does not merely shrink an effect; it can make a genuinely positive effect look negative. The adjusted estimate still lands near 1.14.

4. Now break one arrow. Change 0.3 * age to 0 * age — age no longer affects xenophobia at all — but leave -0.1 * age in Step 2, so the two groups are still badly imbalanced on age. What is the omitted variable bias now?

   With 0 * age, the bivariate estimate is and the imbalance \(\beta_{D \rightarrow C}\) is still about years. The bias is essentially zerostill large and negativelarger than beforeimpossible to compute.

Why?

Because the bias is a product: \(\beta_{C \rightarrow Y} \times \beta_{D \rightarrow C}\). Multiply a big imbalance by zero and you get zero. An imbalanced variable is only a confounder if it also affects the outcome. This is why "the groups differ on X" is not, by itself, an argument that your estimate is biased — you have to make the case for both arrows.

5. Flip the selection. Go back to 0.3 * age, but change -0.1 * age in Step 2 to 0.1 * age — now older people are more likely to end up in the treatment group. Predict the direction first, then run it.

   The bivariate estimate becomes , which is too large — the bias is now positivetoo small — the bias is still negativeexactly rightunchanged

R solution

set.seed(831983) # So everyone gets the SAME "random" numbers
toy_flip <- tibble(
  age = rnorm(n = 500, mean = 40, sd = 5)) %>%     # Step 1: age
  mutate(
    # Step 2: who gets treated (older people are interviewed earlier)
    int_day = (0.1 * age) + rnorm(n = 500, mean = 0, sd = 1),
    bali = case_when(
      int_day < mean(int_day) ~ 0,   # interviewed before the attack
      int_day >= mean(int_day) ~ 1), # interviewed after the attack
    # Step 3: the outcome. THE TRUTH: the Bali effect is 1
    xeno = (1 * bali) + (0.3 * age) + rnorm(n = 500, mean = 0, sd = 1)) #<< Step 2 now uses 0.1

lm_robust(xeno ~ bali, data = toy_flip)
lm_robust(xeno ~ bali + age, data = toy_flip)

Both arrows are now positive (\(+ \times + = +\)), so the bias is positive and the raw estimate overstates the true effect of 1. Same data-generating machinery, one sign flipped, and the error changes direction — exactly as the sign table on the slides predicted.

6. The last knob is the interesting one. Change Step 2 to 0 * age: interview timing no longer depends on age at all. Fit both models and compare them.

R solution — don't peek too early ;) !

set.seed(831983) # So everyone gets the SAME "random" numbers
toy_rand <- tibble(
  age = rnorm(n = 500, mean = 40, sd = 5)) %>%     # Step 1: age
  mutate(
    # Step 2: who gets treated (older people are interviewed earlier)
    int_day = (0 * age) + rnorm(n = 500, mean = 0, sd = 1),
    bali = case_when(
      int_day < mean(int_day) ~ 0,   # interviewed before the attack
      int_day >= mean(int_day) ~ 1), # interviewed after the attack
    # Step 3: the outcome. THE TRUTH: the Bali effect is 1
    xeno = (1 * bali) + (0.3 * age) + rnorm(n = 500, mean = 0, sd = 1)) #<< Step 2 now uses 0

modelsummary(
  list("Bivariate" = lm_robust(xeno ~ bali, data = toy_rand),
       "Adjusted"  = lm_robust(xeno ~ bali + age, data = toy_rand)),
  coef_rename = c("bali" = "Exposed to Bali", "age" = "Age"),
  stars = TRUE, gof_map = c("nobs"),
  output = "kableExtra"
)

	Bivariate	Adjusted
(Intercept)	12.064***	−0.262
	(0.108)	(0.392)
Exposed to Bali	1.044***	1.044***
	(0.158)	(0.091)
Age		0.308***
		(0.010)
Num.Obs.	500	500
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

Write one sentence as a # comment saying what is special about this case — then compare:

Model answer

With 0 * age in Step 2 the two models give essentially the same treatment estimate (≈ 1.04), and both are close to the true effect of 1. Nothing is confounded, because age no longer predicts who gets treated — the imbalance term \(\beta_{D \rightarrow C}\) is approximately zero, so the bias is zero no matter how strongly age affects xenophobia.

This is what a randomized controlled trial looks like (Lecture 6). Randomisation sets \(\beta_{D \rightarrow C}\) to zero for every confounder at once — measured or not — which is why an RCT needs no controls, and why it outranks anything you can do with a regression.

7. Bonus, for the fast. Everything above kept set.seed(831983). Re-run the default simulation with three or four different seeds and collect the adjusted estimate each time. How close to the true value of 1 does it get, and why is it never exactly 1?

R solution

for (seed in c(1, 42, 2024, 831983)) { # Repeat the lecture's simulation with four seeds
  set.seed(seed)
  toy <- tibble(age = rnorm(n = 500, mean = 40, sd = 5)) %>%
    mutate(int_day = (-0.1 * age) + rnorm(n = 500, mean = 0, sd = 1),
           bali = case_when(int_day < mean(int_day) ~ 0,
                            int_day >= mean(int_day) ~ 1),
           xeno = (1 * bali) + (0.3 * age) + rnorm(n = 500, mean = 0, sd = 1))
  print(coef(lm_robust(xeno ~ bali + age, data = toy))["bali"]) # The adjusted estimate
}

The adjusted estimate scatters around 1 without ever landing exactly on it — that is sampling error (Lecture 3), not bias. Bias means being systematically off in one direction; sampling error means bouncing around the truth. Multiple OLS removed the bias here; only a larger sample would shrink the bounce.

8. Discuss with your neighbour. In this exercise you always knew the truth. Work through:

   • In real research, which of the three quantities in the formula can you actually estimate from your data, and which would you have to argue for on theoretical grounds?
   • You adjust for a confounder and your estimate barely moves. Give two different explanations for that.
   • Someone claims: "I controlled for everything in the data set, so my estimate is causal." What is wrong with that sentence?

Talking points

• For an observed confounder you can estimate both arrows — \(\beta_{C \rightarrow Y}\) from the adjusted model and \(\beta_{D \rightarrow C}\) from a balance table. For an unobserved one you can estimate neither; you can only reason about the plausible sign of each arrow and thus the likely direction of the bias. That reasoning is a normal, respectable part of a discussion section.
• Two explanations: (i) the variable genuinely was not much of a confounder — one of the two arrows is near zero, as in task 4; or (ii) it is a confounder, but its bias is being offset by another confounder pulling the opposite way. A stable coefficient is evidence, not proof.
• "Everything in the data set" is not "everything". The sentence quietly assumes the data set happens to contain every variable that affects both treatment and outcome — which is an assumption about the world, not a fact about the data, and it is untestable. The honest version is: "I adjusted for the confounders I could observe, and here is why I think the remaining ones are unlikely to overturn this."