Week 3 seminar

Author

Joshua Loftus

Note

This is the student version of the notebook. The code you write goes where a line says # ... or has a ____. A code chunk with a blank starts with #| eval: false, so the file renders as it stands; delete that line once you have filled the blank. Working through the chunks one at a time in RStudio (Ctrl/Cmd+Enter) is unaffected. Answers are revealed in the seminar.

library(tidyverse)       # should've already install.packages("tidyverse")
library(palmerpenguins)  # might need to install.packages("palmerpenguins")
library(broom)           # might need to install.packages("broom")
library(yardstick)       # might need to install.packages("yardstick")
theme_set(theme_minimal())

Before the seminar

The lecture proved one result, the inverse propensity weighting theorem, which you prove again with the slides closed (exercise 4). It stated the optimal threshold theorem and the adjustment theorem and left their proofs to you (exercises 1 and 3), and it stated the confounding gap without deriving it (exercise 2). Do all four on paper before the seminar; the seminar’s problems are variations on them. You may use, without proof, two facts from ST102: the tower property, \(\mathbb E[W] = \mathbb E\big[\mathbb E[W \mid X]\big]\), and that if \(V\) is independent of \(X\), then \(\mathbb E[k(X, V) \mid X = x] = \mathbb E[k(x, V)]\) for any function \(k\) (freeze \(X\) at \(x\); the distribution of \(V\) does not move). A third fact, taking out what is known: \(\mathbb E[f(X)\,W \mid X] = f(X)\,\mathbb E[W \mid X]\). Throughout, assume that every expectation that appears exists: finite means, and finite variances wherever a variance or covariance appears.

WarningHints

Each exercise has hints, collapsed below it. Do not open a hint unless you have been stuck for at least a few minutes, and open them in order. Exercise 4’s proof is on the lecture slides. Solutions to exercises 1 to 3 are shown in the seminar if you ask.

Exercise 1 (optimal threshold theorem). Let \(p(x) = P(Y = 1 \mid X = x)\). A false positive (predict \(1\), truth \(0\)) costs \(c_{FP} \ge 0\), a false negative (predict \(0\), truth \(1\)) costs \(c_{FN} \ge 0\), not both are \(0\), and a correct prediction costs nothing. A rule \(\delta\) assigns a prediction \(\delta(x) \in \{0, 1\}\) to every \(x\).

  1. At a fixed \(x\), write the expected cost of predicting \(1\) and the expected cost of predicting \(0\).

  2. Show that predicting \(1\) has the smaller expected cost exactly when \(p(x) > c^\star = c_{FP}/(c_{FP} + c_{FN})\). What happens at \(p(x) = c^\star\)?

  3. Let \(\delta^\star\) be the rule from (b). Show that \(\mathbb E\big[\text{cost}(\delta^\star(X), Y)\big] \le \mathbb E\big[\text{cost}(\delta(X), Y)\big]\) for every rule \(\delta\).

  4. With equal costs, what is \(c^\star\), and what is the probability of a wrong label at \(x\) under \(\delta^\star\)?

Given \(X = x\), \(Y\) is \(1\) with probability \(p(x)\). If you predict \(1\), when is that a mistake, and what does the mistake cost?

  1. Compare the two expected costs and solve the inequality for \(p(x)\). “Not both \(0\)” is what lets you divide.
  1. Condition on \(X\) first (the tower property), as in week 2’s exercise 1. What do you know about the inner expectation, at every \(x\)?

Exercise 2 (confounding gap). The lecture’s DAG made linear. \(Z\), \(\varepsilon_T\) and \(\varepsilon_Y\) are independent and normal (so that conditioning on \(T = 0\) and on \(T = 1\) makes sense), \(\varepsilon_T\) and \(\varepsilon_Y\) have mean \(0\), \(\mathrm{Var}(T) > 0\), and \[ T = \alpha Z + \varepsilon_T, \qquad Y = \tau T + \beta Z + \varepsilon_Y. \]

  1. Under \(do(T = t)\), which variables keep their distributions? Show that \(\mathbb E[Y \mid do(T = t)] = \tau t + \beta\, \mathbb E Z\), so that moving \(t\) from \(0\) to \(1\) changes the mean of \(Y\) by \(\tau\).

  2. As observed, show that \(\mathbb E[Y \mid T = t] = \tau t + \beta\, \mathbb E[Z \mid T = t]\), and hence \[ \mathbb E[Y \mid T = 1] - \mathbb E[Y \mid T = 0] = \tau + \beta\,\big(\mathbb E[Z \mid T = 1] - \mathbb E[Z \mid T = 0]\big). \] Where did you use independence?

  3. Show that the population slope of the least-squares line of \(Y\) on \(T\) alone, \(\mathrm{Cov}(T, Y)/\mathrm{Var}(T)\), is \(\tau + \beta\,\delta\) with \(\delta = \mathrm{Cov}(Z, T)/\mathrm{Var}(T)\). Which quantity in week 2’s omitted-variable bias plays the role of \(\delta\)?

  4. When is the gap zero?

  1. Delete \(T\)’s equation and set \(T = t\). Is \(Z\) upstream or downstream of \(T\)?
  1. Take \(\mathbb E[\,\cdot \mid T = t]\) of \(Y\)’s equation term by term. \(T\) is a function of \(Z\) and \(\varepsilon_T\); what does that say about \(\varepsilon_Y\) and \(T\)?
  1. Expand \(\mathrm{Cov}(T, Y)\) using \(Y\)’s equation. \(\mathrm{Cov}(T, \varepsilon_Y) = 0\) for the reason you found in (b).

Exercise 3 (adjustment theorem). The lecture stated this one; prove it. Structural model: \(Z\), \(\varepsilon_T\) and \(\varepsilon_Y\) independent, \(T = h(Z, \varepsilon_T)\), \(Y = g(T, Z, \varepsilon_Y)\); under \(do(T = t)\) the equation for \(T\) is deleted, \(T\) is set to \(t\), every other equation stays, and the noises keep their distribution. Suppose \(T\) takes finitely many values and \(P(T = t \mid Z) > 0\) with probability one. Show that \[ \mathbb E[Y \mid do(T = t)] = \mathbb E_Z\Big[\,\mathbb E[Y \mid T = t, Z]\,\Big], \] and say where each assumption is used: the independence of the noises, \(Z\) not being a consequence of \(T\), and the condition on \(P(T = t \mid Z)\).

When \(Z\) takes finitely many values, the formula is commonly written as a sum. The reading writes it this way, for probabilities of \(Y\) rather than its mean, and calls it the adjustment formula. It is the same formula: \[ \mathbb E[Y \mid do(T = t)] = \sum_z \mathbb E[Y \mid T = t, Z = z]\, P(Z = z). \]

For a continuous \(T\) (as in problem 2 and world D), \(T = t\) exactly has probability \(0\), and the theorem needs two more conditions: \((T, Z)\) has a continuous joint density that is positive everywhere, and \(m_t(z) = \mathbb E[g(t, z, \varepsilon_Y)]\) is continuous in \(t\) and \(z\). Then \(\mathbb E[Y \mid T = t, Z = z]\) is the limit of the average \(Y\) over units with \((T, Z)\) near \((t, z)\), that limit is \(m_t(z)\), and the same proof goes through. The linear normal models in this notebook satisfy both conditions.

Start from the intervened world. What is \(Y\) under \(do(T = t)\), and which variables keep their distributions? Give the average over the noise with \(Z\) frozen a name: \(m_t(z) = \mathbb E\big[g(t, z, \varepsilon_Y)\big]\).

Use the second ST102 fact, with \(V = \varepsilon_Y\) and \(Z\) in the role of \(X\), to show \(\mathbb E[Y \mid do(T = t)] = \mathbb E_Z\big[m_t(Z)\big]\).

Now the observed world: on \(\{T = t, Z = z\}\), \(Y = g(t, z, \varepsilon_Y)\). Why is \(\varepsilon_Y\) independent of \((T, Z)\), and what does that give for \(\mathbb E[Y \mid T = t, Z = z]\)?

Why is the distribution of \(Z\) the same in the observed and the intervened worlds?

Exercise 4 (inverse propensity weighting theorem). The lecture proved this one. Close the slides and prove it again. The structural model of exercise 3, with \(T \in \{0, 1\}\) and the propensity score \(\pi(z) = P(T = 1 \mid Z = z)\), where \(0 < \pi(Z) < 1\) with probability one.

  1. Show that \(\mathbb E[Y \mid do(T = 1)] = \mathbb E\big[T\,Y/\pi(Z)\big]\), and say where you used the independence of the noises, \(Z\) not being a consequence of \(T\), and \(\pi(Z) > 0\).

  2. Show that \(\mathbb E[Y \mid do(T = 0)] = \mathbb E\big[(1 - T)\,Y/(1 - \pi(Z))\big]\).

  3. Without the structural model, show by conditioning on \(Z\) that \(\mathbb E\big[T\,Y/\pi(Z)\big] = \mathbb E_Z\big[\,\mathbb E[Y \mid T = 1, Z]\,\big]\). What does that say about the two theorems?

  1. Start from what you want. Under \(do(T = 1)\), what is \(Y\), and what is random? To bring \(T\) in, multiply and divide by \(\pi(Z)\).
  1. Replace the \(\pi(Z)\) on top by \(\mathbb E_{\varepsilon_T}[T \mid Z, \varepsilon_Y]\). Why are they equal? (The second fact, with \(V = \varepsilon_T\).) Then take the fixed factor inside and use the tower property.
  1. At the end, check that \(T\,g(1, Z, \varepsilon_Y) = T\,Y\) when \(T = 0\) and when \(T = 1\).
  1. \(T\,Y = Y\,\mathbb 1[T = 1]\), so \(\mathbb E[T\,Y \mid Z] = P(T = 1 \mid Z)\, \mathbb E[Y \mid T = 1, Z]\).

Problems

Problem 1 (costs 5 and 1). A conservation survey classifies penguins from a photograph. Labeling a Chinstrap as an Adelie (a miss) wastes a return visit, cost 5; labeling an Adelie as a Chinstrap (a false alarm) wastes a phone call, cost 1. Write \(p(x) = P(\text{Chinstrap} \mid x)\).

  1. Find the threshold on \(p(x)\) above which predicting Chinstrap has the smaller expected cost. What happens at the threshold itself?

  2. A fitted model gives \(\hat p = 0.25\) for one bird. What does the rule predict, and what would have to be true of \(\hat p\) for that prediction to be the best one at that \(x\)?

  3. The survey hands the same \(\hat p\) to two field teams. For team A a miss costs 5 and a false alarm 1; for team B both cost 1. For which values of \(\hat p\) do the two teams make different calls?

  4. Say in one sentence what the survey decided when it chose 5 and 1.

Problem 2 (the lecture’s DAG, with numbers). Let \(Z\), \(\varepsilon_T\) and \(\varepsilon_Y\) be independent standard normal variables, \(T = Z + \varepsilon_T\) and \(Y = T + 2Z + \varepsilon_Y\).

  1. What is the effect on \(\mathbb E[Y]\) of raising \(t\) by one under \(do(T = t)\)?

  2. Compute \(\mathrm{Var}(T)\), \(\mathrm{Cov}(T, Y)\) and the population slope of the regression of \(Y\) on \(T\) alone, and compare the slope with (a).

  3. Recover the answer to (a) from the distribution of the observed \((T, Z, Y)\) alone. Which result are you using, and where do you need the noises to be independent?

Problem 3 (inverse propensity weighting by hand). \(Z\) is \(0\) or \(1\), each with probability \(\tfrac12\). Treatment depends on \(Z\): \(T = h(Z, \varepsilon_T)\) with \(\pi(0) = 0.2\) and \(\pi(1) = 0.8\). The outcome is \(Y = 2T + 4Z + \varepsilon_Y\), with \(\varepsilon_Y\) of mean \(0\), and \(Z\), \(\varepsilon_T\), \(\varepsilon_Y\) independent.

  1. What is \(\mathbb E[Y \mid do(T = 1)] - \mathbb E[Y \mid do(T = 0)]\)?

  2. Picture 100 units in exactly these proportions. How many are treated with \(Z = 0\), treated with \(Z = 1\), untreated with \(Z = 0\), untreated with \(Z = 1\)? Compute \(\mathbb E[Y \mid T = 1] - \mathbb E[Y \mid T = 0]\) and compare it with (a).

  3. Give each unit the weight that inverse propensity weighting gives it, and multiply each count from (b) by its weight. What do the weighted treated units look like? Compute the inverse propensity weighting estimates of \(\mathbb E[Y \mid do(T = 1)]\) and \(\mathbb E[Y \mid do(T = 0)]\) and compare with (a).

  4. Which units get the largest weights? What would happen to them if \(\pi(0)\) were \(0.01\)?

Coding

1. Penguins: a probability and a threshold

Adelie against Chinstrap, from the bill. Keep the two species and the birds with complete measurements.

pg <- palmerpenguins::penguins |>   # modeldata also has a penguins; be explicit
  filter(species != "Gentoo") |>
  drop_na(bill_length_mm, bill_depth_mm, flipper_length_mm) |>
  mutate(chinstrap = as.numeric(species == "Chinstrap"))
ggplot(pg, aes(bill_length_mm, bill_depth_mm, colour = species)) +
  geom_point() +
  scale_colour_viridis_d(end = 0.7)

Step 1. How many birds, and what proportion are Chinstraps? What accuracy would a classifier that always predicts the more common species achieve?

# ...

Step 2. Fit the logistic regression of chinstrap on the two bill measurements. Look at its coefficients as they are and exponentiated (tidy() has an exponentiate argument).

# fit_bill <- ____
# ...

Checkpoint (run it; if it stops, fix step 2 before going on):

stopifnot(inherits(fit_bill, "glm"), length(coef(fit_bill)) == 3,
          all(fitted(fit_bill) >= 0 & fitted(fit_bill) <= 1))

You may have seen a warning that fitted probabilities are numerically 0 or 1: the two species are nearly separable on these two measurements (step 3 counts the errors). This is the lecture’s “Challenges” slide in real life; next week says what it does to the fitting.

Question 1. Read the exponentiated coefficient on bill_length_mm. What happens to the odds of Chinstrap when bill length increases by one millimeter and bill depth is held fixed? Is that the change in the probability? Is it what would happen to a bird whose bill grew by a millimeter?

Step 3. The confusion matrix at threshold \(\tfrac12\), with table(). From it, compute the accuracy, the true positive rate and the false positive rate.

# p_hat <- ____
# tab <- table(truth = ____, predicted = ____)
# accuracy <- ____; tpr <- ____; fpr <- ____
# ...

Checkpoint:

stopifnot(sum(tab) == nrow(pg), tpr >= 0, tpr <= 1, fpr >= 0, fpr <= 1)

Step 4. The bill model is nearly perfect, so it cannot show the trade-off. Fit a weaker model on flipper_length_mm alone, then compute the ROC curve yourself: at every threshold where it can change (0, then each distinct fitted probability), the true positive rate and the false positive rate. Plot one against the other. The thresholds are set up; you write the rest.

# fit_flipper <- ____
# pg <- pg |> mutate(p_flipper = ____)
# thresholds <- c(0, sort(unique(pg$p_flipper)))   # 0, then every distinct fitted probability
# roc_by_hand <- map_dfr(thresholds, function(t) {
#   ____
#   tibble(t = t, tpr = ____, fpr = ____)
# })
# ggplot(roc_by_hand, aes(____, ____)) + geom_path() + coord_equal()
# ...

Checkpoint:

stopifnot(nrow(roc_by_hand) == length(thresholds),
          all(roc_by_hand$tpr >= 0 & roc_by_hand$tpr <= 1),
          roc_by_hand$tpr[1] == 1, roc_by_hand$fpr[1] == 1)

Question 2. The area under your curve equals the probability that a randomly chosen Chinstrap gets a higher fitted probability than a randomly chosen Adelie, with a tie counting half. Why does accumulating one point per threshold produce that probability? Why are ties certain to occur with this model?

Step 5. Costs. A miss costs 5, a false alarm 1. For each threshold on a grid from 0 to 1 in steps of 0.01, compute the average cost per bird on the 219 birds, find the threshold that minimizes it, and compare with problem 1’s threshold.

grid <- seq(0, 1, by = 0.01)
# costs <- map_dfr(grid, function(t) {
#   ____
#   tibble(t = t, cost = ____)
# })
# ...

Checkpoint:

stopifnot(nrow(costs) == length(grid), all(costs$cost >= 0))

Question 3. The cost-minimizing threshold is well below \(\tfrac12\), and the classifier now raises many false alarms. Who decided that, and in what sense is this classifier “better” than the one at \(\tfrac12\)? Is the number you found the population optimum?

Optional, can be done at home: the same with yardstick

yardstick computes the confusion matrix, the ROC curve and the AUC. It wants factors, and the event_level = "second" argument tells it that 1 is the class of interest.

# pg_y <- pg |> mutate(truth = factor(chinstrap),
#                      predicted = factor(as.numeric(p_hat > 0.5), levels = c(0, 1)))
# conf_mat(pg_y, truth = truth, estimate = predicted)
# roc_curve(pg_y, truth, p_flipper, event_level = "second") |> autoplot()
# roc_auc(pg_y, truth, ____, event_level = "second")
# ...

2. Simulated worlds: adjustment and inverse propensity weighting

Suppose we find a statistically significant correlation between \(X\) and \(Y\). What does it mean?

  1. A false positive (a spurious correlation)
  2. \(X\) causes \(Y\)
  3. \(Y\) causes \(X\)
  4. Both have a common cause \(U\), possibly unobserved

Without experimental data, explanations 2 to 4 can be statistically indistinguishable, and they have importantly different consequences. Worlds A, C and B below are explanations 2, 3 and 4, built so that the slope of \(Y\) on \(X\) is about the same in all three. In the lecture the target variable was called \(T\); here it is X, because T is TRUE in R. The function draws the noise once and keeps it, so that step 2 can reuse it. World D is the lecture’s DAG with the confounder \(Z\) observed: \(Z \to X\), \(Z \to Y\), \(X \to Y\). Step 4 adds a fifth world, E, with the same DAG and a binary treatment.

simulate_world <- function(world, n = 500) {
  e_u <- rnorm(n)
  e_x <- rnorm(n)
  e_y <- rnorm(n)
  if (world == "A") {           # X causes Y
    X <- e_x
    Y <- 2 * X + e_y
    Z <- NA
  } else if (world == "B") {    # U causes both; U is not observed
    U <- e_u
    X <- U + 0.5 * e_x
    Y <- 2.5 * U + e_y
    Z <- NA
  } else if (world == "C") {    # Y causes X
    Y <- 2 * e_y
    X <- 0.4 * Y + 0.4 * e_x
    Z <- NA
  } else {                      # D: Z causes both, and X causes Y; Z is observed
    Z <- e_u
    X <- 5 - 2 * Z + e_x
    Y <- 4 + 2 * Z - 3 * X + e_y
  }
  tibble(world = world, X = X, Y = Y, Z = Z, e_u = e_u, e_x = e_x, e_y = e_y)
}
worlds <- bind_rows(simulate_world("A"), simulate_world("B"), simulate_world("C"), simulate_world("D"))

Step 1. In each world, fit lm(Y ~ X) and report the slope.

# worlds |>
#   group_by(world) |>
#   summarise(slope = coef(lm(____))[2])
# ...

Question 4. Before you intervene: for each world, write down what you expect to happen to the average of \(Y\) when \(X\) is set to \(X + 1\) for every unit. Use each world’s equations, and exercise 2 for world D.

Step 2. Now intervene: set \(X\) to \(X + 1\) for every unit and recompute \(Y\) with the same noise, following each world’s equations. Write the four Y_do lines and compute the average change in \(Y\) in each world.

Note: recomputing each unit’s \(Y\) with its own noise is different from the intervention calculations in the lecture, which use only the distribution of the noise. We will learn more about it later in the course. Here it makes the average change exact, with no simulation noise.

intervene <- function(d) {
  d |> mutate(
    X_do = X + 1,
    # Y_do = case_when(
    #   world == "A" ~ ____,
    #   world == "B" ~ ____,
    #   world == "C" ~ ____,
    #   world == "D" ~ ____
    # )
    # ...
  )
}
# worlds_do <- intervene(worlds)
# worlds_do |> group_by(world) |> summarise(average_change = ____)
# ...

Checkpoint:

stopifnot(all(c("X_do", "Y_do") %in% names(worlds_do)), !anyNA(worlds_do$Y_do))

Compare the average changes with your predictions in question 4.

Step 3. World D, the adjustment theorem in code. Regress Y on X and Z and read off the coefficient on X. Then use the same regression as your estimate of \(\mathbb E[Y \mid X = x, Z = z]\) and do what the adjustment theorem says, to estimate the average change in \(Y\) when every unit’s \(X\) is raised by one.

d <- worlds |> filter(world == "D")
# fit_d <- ____
# ...

Question 5. Compare the coefficient on X with the averaged difference, and both with the effect you found in step 2. Which result says the coefficient estimates the effect, which says the averaged difference does, and what about world D made both work? What would go wrong in world B, where the confounder e_u is in the table but would not be observed in practice?

Step 4. World E: the same DAG, with a binary treatment. Each unit is treated (X = 1) with probability \(\pi(Z)\), where \(\pi(z) = 1/(1 + e^{-(-0.5 + 1.5 z)})\), and \(Y = 1 + 2X + 3Z + \varepsilon_Y\). Compute three estimates of the effect of treatment: the naive difference in means between treated and untreated, the inverse propensity weighting estimate, and regression adjustment. Then look at the five units with the largest weights.

simulate_world_e <- function(n = 2000) {
  Z <- rnorm(n)
  e_y <- rnorm(n)
  X <- rbinom(n, 1, plogis(-0.5 + 1.5 * Z))   # treated with probability pi(Z)
  Y <- 1 + 2 * X + 3 * Z + e_y
  tibble(X = X, Y = Y, Z = Z)
}
d_e <- simulate_world_e()
# fit_pi <- ____
# d_e <- d_e |> mutate(pi_hat = ____, w = ____)
# c(naive = ____, weighting = ____, regression = ____)
# d_e |> slice_max(w, n = 5)
# ...

Checkpoint, then a diagnostic. As in problem 3(c), the weighted treated units and the weighted untreated units should each count as about the whole sample, so both numbers printed should be near 1; a few large weights can move them.

stopifnot(all(d_e$pi_hat > 0 & d_e$pi_hat < 1), all(d_e$w >= 1))
c(treated = sum(d_e$w[d_e$X == 1]), untreated = sum(d_e$w[d_e$X == 0])) / nrow(d_e)

Question 6. What is the effect of treatment in world E, which of the three estimates are near it, and which result explains each one? Which units have the largest weights, and what is unusual about them? If the unit with the largest weight had not been in the sample, and the other weights stayed as they are, by how much would the weighting estimate move?

Optional, can be done at home: another dataset

Attrition (modeldata::attrition; fictional, made up by IBM data scientists, used only as an exercise): whether an employee leaves (Attrition), with a threshold-and-costs story (a missed leaver costs a replacement, a false alarm costs a conversation). Leavers are 16 percent, so the accuracy-against-cost point is sharper: at threshold \(\tfrac12\) almost nobody is predicted to leave.

library(modeldata)
data(attrition)
att <- attrition |> mutate(left = as.numeric(Attrition == "Yes"))
# fit_att <- ____   # predictors: OverTime, MonthlyIncome, Age, JobSatisfaction
# then steps 3, 4 and 5, with costs 10 (a missed leaver) and 1 (a false alarm)
# ...