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())Week 3 seminar (complete)
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\).
At a fixed \(x\), write the expected cost of predicting \(1\) and the expected cost of predicting \(0\).
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\)?
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\).
With equal costs, what is \(c^\star\), and what is the probability of a wrong label at \(x\) under \(\delta^\star\)?
NoteHint 1
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?
NoteHint 2
- Compare the two expected costs and solve the inequality for \(p(x)\). “Not both \(0\)” is what lets you divide.
NoteHint 3
- 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\)?
TipReveal solution
Predicting \(1\) is a mistake when \(Y = 0\), which has probability \(1 - p(x)\), and then costs \(c_{FP}\); predicting \(0\) is a mistake when \(Y = 1\), probability \(p(x)\), cost \(c_{FN}\): \[ \mathbb E[\text{cost} \mid \text{predict } 1,\ X = x] = c_{FP}\,(1 - p(x)), \qquad \mathbb E[\text{cost} \mid \text{predict } 0,\ X = x] = c_{FN}\, p(x). \]
Predicting \(1\) is cheaper when \[ c_{FN}\,p > c_{FP}\,(1 - p) \iff p\,(c_{FP} + c_{FN}) > c_{FP} \iff p > \frac{c_{FP}}{c_{FP} + c_{FN}}, \] dividing by \(c_{FP} + c_{FN} > 0\). At \(p(x) = c^\star\) the two expected costs are equal, so either prediction is optimal.
For any rule \(\delta\), by the tower property, \[ \mathbb E\big[\text{cost}(\delta(X), Y)\big] = \mathbb E_X\Big[\,\mathbb E\big[\text{cost}(\delta(X), Y) \mid X\big]\Big]. \] At \(X = x\) the inner expectation is the expected cost of the prediction \(\delta(x)\) from (a), and by (b) it is at least the expected cost of \(\delta^\star(x)\). An inequality that holds at every \(x\) survives averaging over \(X\). \(\blacksquare\) This is week 2’s step: minimize the expected loss at each \(x\), then average.
\(c^\star = \tfrac12\): predict the more probable class, the Bayes classifier. The probability of a wrong label at \(x\) is \(1 - p(x)\) if you predict \(1\) and \(p(x)\) if you predict \(0\), so under \(\delta^\star\) it is \(\min(p(x), 1 - p(x))\). Equal costs make the expected cost a multiple of the probability of error, so by (c) no rule has a smaller overall probability of error.
Common errors: the ratio upside down, \(c_{FN}/(c_{FP} + c_{FN})\) (check an extreme: if a miss were infinitely costly the threshold must go to \(0\), and only \(c_{FP}/(c_{FP} + c_{FN})\) does that); choosing one threshold to minimize the overall cost directly, without the argument at each \(x\); reading the theorem as a statement about \(\hat p\). It is about the true \(p(x)\); with \(\hat p\) in its place the rule is a plug-in rule, only as good as the estimate.
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. \]
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\).
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?
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\)?
When is the gap zero?
NoteHint 1
- Delete \(T\)’s equation and set \(T = t\). Is \(Z\) upstream or downstream of \(T\)?
NoteHint 2
- 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\)?
NoteHint 3
- Expand \(\mathrm{Cov}(T, Y)\) using \(Y\)’s equation. \(\mathrm{Cov}(T, \varepsilon_Y) = 0\) for the reason you found in (b).
TipReveal solution
\(Z\) is upstream of \(T\) and the noises are their own, so \(Z\) and \(\varepsilon_Y\) keep their distributions: \(Y = \tau t + \beta Z + \varepsilon_Y\), and \(\mathbb E[Y \mid do(T = t)] = \tau t + \beta\,\mathbb E Z + 0\). The difference between \(t = 1\) and \(t = 0\) is \(\tau\), the coefficient on \(T\) in \(Y\)’s own equation.
\(\mathbb E[Y \mid T = t] = \tau t + \beta\,\mathbb E[Z \mid T = t] + \mathbb E[\varepsilon_Y \mid T = t]\). \(T = \alpha Z + \varepsilon_T\) is a function of \((Z, \varepsilon_T)\), and \(\varepsilon_Y\) is independent of \((Z, \varepsilon_T)\), so \(\varepsilon_Y\) is independent of \(T\) and \(\mathbb E[\varepsilon_Y \mid T = t] = \mathbb E\,\varepsilon_Y = 0\). That is where independence is used. Subtracting the expressions at \(t = 1\) and \(t = 0\) gives the gap. \(\blacksquare\)
\(\mathrm{Cov}(T, Y) = \tau\,\mathrm{Var}(T) + \beta\,\mathrm{Cov}(T, Z) + \mathrm{Cov}(T, \varepsilon_Y)\), and the last term is \(0\). Dividing by \(\mathrm{Var}(T)\) gives \(\tau + \beta\,\delta\). Here \(\delta\), the population slope of \(Z\) on \(T\), plays the role of \(\hat\delta\) (the slope of the omitted \(\mathbf x_2\) on the kept \(\mathbf x_1\)), \(\tau\) plays the role of \(\beta_1\) and \(\beta\) that of \(\beta_2\): omitted-variable bias, read causally. With these equations \(\delta = \alpha\,\mathrm{Var}(Z)/(\alpha^2\,\mathrm{Var}(Z) + \mathrm{Var}(\varepsilon_T))\).
When \(\beta = 0\) (\(Z\) does not affect \(Y\)), or when the mean of \(Z\) does not move with \(T\); for example \(\alpha = 0\), when \(T\) does not depend on \(Z\) at all, which is what randomizing \(T\) arranges. (The gap in (b) compares two conditional means; the slope in (c) is a least-squares slope. They agree, gap \(= \beta\delta\), when \(\mathbb E[Z \mid T]\) is linear in \(T\), as it is when everything is normal.)
Common errors: letting \(Z\) change under the intervention; writing \(\mathbb E[\varepsilon_Y \mid T = t] = 0\) without saying why (it needs the independence, which an unobserved common cause of \(T\) and \(Y\) would break); writing the bias as \(\beta\alpha\) instead of \(\beta\delta\).
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.
NoteHint 1
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]\).
NoteHint 2
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]\).
NoteHint 3
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]\)?
NoteHint 4
Why is the distribution of \(Z\) the same in the observed and the intervened worlds?
TipReveal solution
Hint 1: under \(do(T = t)\), \(Y = g(t, Z, \varepsilon_Y)\), with \((Z, \varepsilon_Y)\) distributed as before (\(Z\) is upstream of \(T\)). Hint 2: \(\varepsilon_Y\) is independent of \(Z\), so \(\mathbb E[g(t, Z, \varepsilon_Y) \mid Z = z] = m_t(z)\), and by the tower property \[ \mathbb E[Y \mid do(T = t)] = \mathbb E\big[g(t, Z, \varepsilon_Y)\big] = \mathbb E_Z\big[m_t(Z)\big]. \] Hint 3: \(T = h(Z, \varepsilon_T)\) is a function of \((Z, \varepsilon_T)\) and \(\varepsilon_Y\) is independent of \((Z, \varepsilon_T)\), so \(\varepsilon_Y\) is independent of \((T, Z)\). Wherever \(P(T = t \mid Z = z) > 0\), the event \(\{T = t, Z = z\}\) can be conditioned on, and on it \(Y = g(t, z, \varepsilon_Y)\), so \[ \mathbb E[Y \mid T = t, Z = z] = \mathbb E\big[g(t, z, \varepsilon_Y) \mid T = t, Z = z\big] = m_t(z). \] Hint 4: \(Z\) is upstream of \(T\), so the intervention does not change it, and its distribution is the same in the two worlds. Averaging the last display over it, \[ \mathbb E[Y \mid do(T = t)] = \mathbb E_Z\big[m_t(Z)\big] = \mathbb E_Z\Big[\,\mathbb E[Y \mid T = t, Z]\,\Big]. \qquad \blacksquare \] Where each assumption is used: the independence of the noises in the observed-world step (a common cause of \(T\) and \(Y\) left out of \(Z\) would sit in both \(\varepsilon_T\) and \(\varepsilon_Y\), making them dependent); \(Z\) not a consequence of \(T\) in the intervened-world step (a \(Z\) downstream of \(T\) would change when \(T\) is set, and its distribution would differ between the worlds); \(P(T = t \mid Z) > 0\) with probability one, so that for the values \(z\) the average runs over there are units with \(T = t\) and \(Z = z\) to average (if a set of \(z\) with positive probability had none, the right-hand side would not be defined there). For a continuous \(T\), the note under the exercise replaces this step. If \(Z\) is empty, \(T = h(\varepsilon_T)\) is independent of \(\varepsilon_Y\) and the same lines give \(\mathbb E[Y \mid do(T = t)] = \mathbb E[Y \mid T = t]\): randomization.
Common errors: letting \(Z\) respond to the intervention; asserting \(\mathbb E[Y \mid T = t, Z = z] = m_t(z)\) without the independence of \(\varepsilon_Y\) from \((T, Z)\); reading the theorem as “regress \(Y\) on \(T\) and \(Z\) and take the coefficient”, which is one way to estimate the inside when the conditional mean is linear, not what the theorem says.
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.
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\).
Show that \(\mathbb E[Y \mid do(T = 0)] = \mathbb E\big[(1 - T)\,Y/(1 - \pi(Z))\big]\).
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?
NoteHint 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)\).
NoteHint 2
- 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.
NoteHint 3
- At the end, check that \(T\,g(1, Z, \varepsilon_Y) = T\,Y\) when \(T = 0\) and when \(T = 1\).
NoteHint 4
- \(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]\).
TipReveal solution
- The subscript on each \(\mathbb E\) names what is averaged over. Start from what we want and change one term at a time: \[ \begin{aligned} \mathbb E[Y \mid do(T = 1)] &= \mathbb E_{Z, \varepsilon_Y}\big[g(1, Z, \varepsilon_Y)\big] & \text{(a)} \\ &= \mathbb E_{Z, \varepsilon_Y}\Big[\frac{\pi(Z)}{\pi(Z)}\, g(1, Z, \varepsilon_Y)\Big] & \text{(b)} \\ &= \mathbb E_{Z, \varepsilon_Y}\Big[\frac{\mathbb E_{\varepsilon_T}[T \mid Z, \varepsilon_Y]}{\pi(Z)}\, g(1, Z, \varepsilon_Y)\Big] & \text{(c)} \\ &= \mathbb E_{Z, \varepsilon_Y}\Big[\,\mathbb E_{\varepsilon_T}\Big[\frac{T\, g(1, Z, \varepsilon_Y)}{\pi(Z)} \;\Big|\; Z, \varepsilon_Y\Big]\Big] & \text{(d)} \\ &= \mathbb E_{Z, \varepsilon_T, \varepsilon_Y}\Big[\frac{T\, g(1, Z, \varepsilon_Y)}{\pi(Z)}\Big] & \text{(e)} \\ &= \mathbb E_{T, Z, Y}\Big[\frac{T\,Y}{\pi(Z)}\Big]. & \text{(f)} \end{aligned} \]
- Under \(do(T = 1)\), \(Y = g(1, Z, \varepsilon_Y)\), with \((Z, \varepsilon_Y)\) distributed as before. (b) \(\pi(Z)/\pi(Z) = 1\) because \(\pi(Z) > 0\). (c) \(T = h(Z, \varepsilon_T)\) and \(\varepsilon_T\) is independent of \((Z, \varepsilon_Y)\), so by the second fact \(\mathbb E[T \mid Z = z, \varepsilon_Y = e] = P\big(h(z, \varepsilon_T) = 1\big)\); the same fact with \(Z\) alone gives \(\pi(z) = P\big(h(z, \varepsilon_T) = 1\big)\). (d) \(g(1, Z, \varepsilon_Y)/\pi(Z)\) is fixed once \((Z, \varepsilon_Y)\) is given: taking out what is known, in reverse. (e) The tower property. (f) If \(T = 0\), both \(T\,g(1, Z, \varepsilon_Y)\) and \(T\,Y\) are \(0\); if \(T = 1\), \(Y = g(1, Z, \varepsilon_Y)\). The last line is the same average, written over the observed \((T, Z, Y)\), which are functions of \((Z, \varepsilon_T, \varepsilon_Y)\). \(\blacksquare\)
Where each assumption is used: the independence of \(\varepsilon_T\) from \((Z, \varepsilon_Y)\) in (c) (a common cause of \(T\) and \(Y\) left out of \(Z\) would sit in both noises, and then \(\varepsilon_Y\) would carry information about \(T\) beyond \(Z\)); \(\pi(Z) > 0\) in (b); \(Z\) not a consequence of \(T\) in (a) (a \(Z\) downstream of \(T\) would change when \(T\) is set).
The same steps with \(1 - T\) and \(1 - \pi(Z)\): in (b), multiply and divide by \(1 - \pi(Z)\), which needs \(\pi(Z) < 1\); in (c), \(\mathbb E_{\varepsilon_T}[1 - T \mid Z, \varepsilon_Y] = 1 - \pi(Z)\); in (f), \((1 - T)\,g(0, Z, \varepsilon_Y) = (1 - T)\,Y\).
\(\mathbb E[T\,Y \mid Z] = \mathbb E\big[Y\,\mathbb 1[T = 1] \mid Z\big] = \pi(Z)\,\mathbb E[Y \mid T = 1, Z]\). Since \(1/\pi(Z)\) is a function of \(Z\), the tower property gives \[ \mathbb E\Big[\frac{T\,Y}{\pi(Z)}\Big] = \mathbb E\Big[\frac{\mathbb E[T\,Y \mid Z]}{\pi(Z)}\Big] = \mathbb E_Z\Big[\,\mathbb E[Y \mid T = 1, Z]\,\Big]. \] This holds for any joint distribution with \(\pi(Z) > 0\): the two theorems have the same right-hand side, written two ways, and the causal assumptions are what make it equal to \(\mathbb E[Y \mid do(T = 1)]\). In practice they differ in what you estimate: the adjustment theorem needs a regression model for \(\mathbb E[Y \mid T, Z]\), inverse propensity weighting a regression model for \(\pi(Z)\).
Common errors: writing \(\mathbb E[T \mid Z, \varepsilon_Y] = \pi(Z)\) without saying why; weighting by \(\pi(Z)\) instead of \(1/\pi(Z)\); dropping the \(T\) in the numerator, which weights the untreated units too.
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)\).
Find the threshold on \(p(x)\) above which predicting Chinstrap has the smaller expected cost. What happens at the threshold itself?
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\)?
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?
Say in one sentence what the survey decided when it chose 5 and 1.
NoteHint 1
- Which kind of mistake is the false positive here, and which the false negative? Then exercise 1(b).
NoteHint 2
- Each team has its own threshold. Between the two thresholds, what does each team predict?
TipReveal answer
Predicting Chinstrap is a false alarm when the bird is an Adelie: \(c_{FP} = 1\), \(c_{FN} = 5\). Expected costs \(1 \cdot (1 - p)\) and \(5p\), so Chinstrap is cheaper when \(p > 1/6 \approx 0.167\). At \(p = 1/6\) the two costs are equal and either prediction is optimal.
\(\hat p = 0.25 > 1/6\), so the rule says Chinstrap. That is the best decision at \(x\) if the true \(p(x)\) is at least \(1/6\) too: \(\hat p\) need not be exact, only on the same side of the threshold as \(p(x)\). With an estimate in its place it is a plug-in rule: if the model is badly calibrated there (say the true \(p(x)\) is \(0.1\)) the call is wrong.
Team A’s threshold is \(1/6\), team B’s is \(1/2\). For \(1/6 < \hat p < 1/2\), A calls Chinstrap and B calls Adelie; at the two end points one team is indifferent. Same model, same bird, different decisions: the difference comes from the costs alone.
That five false alarms are worth one miss: a statement about the survey’s resources, not about penguins.
Common errors: swapping the costs (threshold \(5/6\)); “the model predicts Chinstrap”, when the model gives a probability and the rule makes the call.
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\).
What is the effect on \(\mathbb E[Y]\) of raising \(t\) by one under \(do(T = t)\)?
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).
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?
NoteHint 1
- and (b) are exercise 2 with \(\alpha = \tau = 1\) and \(\beta = 2\). For the covariance, \(\mathrm{Cov}(T, Z) = \mathrm{Var}(Z)\) here.
NoteHint 2
- Take \(\mathbb E[\,\cdot \mid T = t, Z = z]\) of \(Y\)’s equation, then exercise 3.
TipReveal answer
Under \(do(T = t)\), \(Y = t + 2Z + \varepsilon_Y\), so \(\mathbb E[Y \mid do(T = t)] = t\) and raising \(t\) by one raises it by \(1\).
\(\mathrm{Var}(T) = 1 + 1 = 2\). \(\mathrm{Cov}(T, Y) = \mathrm{Var}(T) + 2\,\mathrm{Cov}(T, Z) + \mathrm{Cov}(T, \varepsilon_Y) = 2 + 2 + 0 = 4\). Slope \(4/2 = 2 = 1 + 2 \cdot \tfrac12\): the effect \(\tau = 1\) plus the gap \(\beta\delta = 2 \cdot \tfrac12\), with \(\delta = \mathrm{Cov}(Z, T)/\mathrm{Var}(T) = \tfrac12\). The naive slope credits \(T\) with what \(Z\) contributes through selection, and here the gap is as large as the effect.
The adjustment theorem. \(\mathbb E[Y \mid T = t, Z = z] = t + 2z + \mathbb E[\varepsilon_Y \mid T = t, Z = z] = t + 2z\), because \(\varepsilon_Y\) is independent of \((Z, \varepsilon_T)\) and so of \((T, Z)\). Averaging over \(Z\): \(\mathbb E_Z[t + 2Z] = t\), which is (a). The independence is what made \(\mathbb E[\varepsilon_Y \mid T = t, Z = z]\) zero. Normality is what makes this work at every \(t\): given \(Z\), \(T\) has a positive density everywhere, so the adjustment theorem’s condition holds.
Common errors: in (a), letting \(Z\) change when \(T\) is set; in (b), \(\mathrm{Cov}(T, Z) = 0\).
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.
What is \(\mathbb E[Y \mid do(T = 1)] - \mathbb E[Y \mid do(T = 0)]\)?
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).
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).
Which units get the largest weights? What would happen to them if \(\pi(0)\) were \(0.01\)?
NoteHint 1
- Under \(do(T = t)\), \(Z\) keeps its distribution. (b) 50 units have \(Z = 0\), and 20 percent of them are treated. Among the treated, what fraction has \(Z = 1\)?
NoteHint 2
- In the 100 units, \(\mathbb E\big[T\,Y/\pi(Z)\big]\) is \(\tfrac{1}{100}\) times the sum of \(Y/\pi(Z)\) over the treated units. Group them by \(Z\): within a group, \(Y\) averages \(\mathbb E[Y \mid T = 1, Z = z] = 2 + 4z\).
TipReveal answer
Under \(do(T = t)\), \(Y = 2t + 4Z + \varepsilon_Y\), so \(\mathbb E[Y \mid do(T = t)] = 2t + 4 \cdot \tfrac12 = 2t + 2\): \(4\) and \(2\), effect \(2\).
Of the 50 units with \(Z = 0\), 10 are treated and 40 are not; of the 50 with \(Z = 1\), 40 are treated and 10 are not. So \(\mathbb E[Z \mid T = 1] = 40/50 = 0.8\) and \(\mathbb E[Z \mid T = 0] = 0.2\), and \(\mathbb E[Y \mid T = 1] = 2 + 4 \cdot 0.8 = 5.2\), \(\mathbb E[Y \mid T = 0] = 4 \cdot 0.2 = 0.8\). The naive difference is \(4.4 = 2 + 4\,(0.8 - 0.2)\): the effect \(2\) plus a gap of \(2.4\). This is the lecture’s confounding gap; its derivation (exercise 2(b)) used only \(Y\)’s equation and the independence of \(\varepsilon_Y\), not the linear equation for \(T\).
Treated: 10 units with \(Z = 0\) at weight \(1/0.2 = 5\) count as 50; 40 units with \(Z = 1\) at weight \(1/0.8 = 1.25\) count as 50. The weighted treated units have the whole population’s \(Z\): half \(0\), half \(1\). The untreated are the same: \(40 \cdot 1.25 = 50\) and \(10 \cdot 5 = 50\). Then \[ \mathbb E\Big[\frac{T\,Y}{\pi(Z)}\Big] = \frac{10 \cdot 5 \cdot 2 + 40 \cdot 1.25 \cdot 6}{100} = 4, \qquad \mathbb E\Big[\frac{(1 - T)\,Y}{1 - \pi(Z)}\Big] = \frac{40 \cdot 1.25 \cdot 0 + 10 \cdot 5 \cdot 4}{100} = 2, \] which are \(\mathbb E[Y \mid do(T = 1)]\) and \(\mathbb E[Y \mid do(T = 0)]\) from (a). Difference \(2\).
The treated units with \(Z = 0\) and the untreated units with \(Z = 1\): the units whose treatment was unlikely given their \(Z\). Each stands in for the many units like it that went the other way. With \(\pi(0) = 0.01\), about one in a hundred units with \(Z = 0\) is treated, and each gets weight \(100\). The condition \(\pi(Z) > 0\) still holds, but in a sample of a few hundred units one or two treated units with \(Z = 0\) would carry the estimate for half the population, so it would be unstable.
Common errors: weighting by \(\pi(Z)\) instead of \(1/\pi(Z)\); weighting the untreated by \(1/\pi(Z)\) instead of \(1/(1 - \pi(Z))\); dividing by the 50 treated units instead of all 100.
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?
# ...
TipReveal solution
nrow(pg)[1] 219mean(pg$chinstrap)[1] 0.3105023max(mean(pg$chinstrap), 1 - mean(pg$chinstrap))[1] 0.6894977Step 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 <- ____
# ...
TipReveal solution
fit_bill <- glm(chinstrap ~ bill_length_mm + bill_depth_mm, family = binomial, data = pg)
tidy(fit_bill)# A tibble: 3 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) -24.1 13.5 -1.78 0.0744
2 bill_length_mm 2.20 0.687 3.20 0.00136
3 bill_depth_mm -3.98 1.49 -2.67 0.00759tidy(fit_bill, exponentiate = TRUE)# A tibble: 3 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 3.43e-11 13.5 -1.78 0.0744
2 bill_length_mm 9.03e+ 0 0.687 3.20 0.00136
3 bill_depth_mm 1.87e- 2 1.49 -2.67 0.00759Checkpoint (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?
TipReveal answer
The odds of Chinstrap are multiplied by about \(e^{2.20} \approx 9.0\) per millimeter, with bill depth held fixed. That is a multiplier on the odds \(p/(1 - p)\), not a change in \(p\): how much \(p\) moves depends on where it starts (near \(0\) or \(1\) a ninefold change in the odds moves \(p\) very little; near \(\tfrac12\) it moves it a lot). The standard error is large (\(0.69\) on the log-odds scale, so the interval for the multiplier runs from about \(2\) to about \(35\)), because the classes are nearly separable. And “holding depth fixed” is the lecture’s ceteris paribus: a comparison between birds in the fitted model, not what would happen to one bird if its bill grew. Nothing here is causal.
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 <- ____
# ...
TipReveal solution
p_hat <- fitted(fit_bill)
tab <- table(truth = pg$chinstrap, predicted = as.numeric(p_hat > 0.5))
tab predicted
truth 0 1
0 148 3
1 3 65accuracy <- (tab[1, 1] + tab[2, 2]) / sum(tab)
tpr <- tab[2, 2] / sum(tab[2, ])
fpr <- tab[1, 2] / sum(tab[1, ])
c(accuracy = accuracy, tpr = tpr, fpr = fpr) accuracy tpr fpr
0.97260274 0.95588235 0.01986755 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()
# ...
TipReveal solution
fit_flipper <- glm(chinstrap ~ flipper_length_mm, family = binomial, data = pg)
pg <- pg |> mutate(p_flipper = fitted(fit_flipper))
thresholds <- c(0, sort(unique(pg$p_flipper))) # 0, then every distinct fitted probability
roc_by_hand <- map_dfr(thresholds, function(t) {
yhat <- as.numeric(pg$p_flipper > t)
tibble(t = t,
tpr = sum(yhat == 1 & pg$chinstrap == 1) / sum(pg$chinstrap == 1),
fpr = sum(yhat == 1 & pg$chinstrap == 0) / sum(pg$chinstrap == 0))
})
ggplot(roc_by_hand, aes(fpr, tpr)) +
geom_path() +
geom_abline(slope = 1, intercept = 0, linetype = "dotted") +
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?
TipReveal answer
Moving the threshold down past one bird’s score adds that bird to the predicted Chinstraps: a Chinstrap raises the TPR by \(1/n_1\), an Adelie raises the FPR by \(1/n_0\). The area under the path therefore adds up, for each Adelie step to the right, the height of the curve there, which is the fraction of Chinstraps already passed, those scoring above that Adelie. Summing over Adelies and dividing by \(n_0\) gives the fraction of (Chinstrap, Adelie) pairs in which the Chinstrap scores higher; a tie is a diagonal step and contributes half. Ties are certain because flipper_length_mm is recorded in whole millimeters and \(\hat p\) is a monotone function of it, so birds with the same flipper length have the same score; about 4 percent of the Chinstrap and Adelie pairs tie. By hand:
p1 <- pg$p_flipper[pg$chinstrap == 1]
p0 <- pg$p_flipper[pg$chinstrap == 0]
mean(outer(p1, p0, ">") + 0.5 * outer(p1, p0, "=="))[1] 0.7337846The area under your curve, by the trapezoid rule (a tie is a diagonal step, so it counts half), is the same number:
r <- arrange(roc_by_hand, fpr, tpr)
sum(diff(r$fpr) * (head(r$tpr, -1) + tail(r$tpr, -1)) / 2)[1] 0.7337846Step 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 = ____)
# })
# ...
TipReveal solution
grid <- seq(0, 1, by = 0.01)
costs <- map_dfr(grid, function(t) {
yhat <- as.numeric(pg$p_flipper > t)
tibble(t = t,
cost = (5 * sum(yhat == 0 & pg$chinstrap == 1) + 1 * sum(yhat == 1 & pg$chinstrap == 0)) / nrow(pg))
})
ggplot(costs, aes(t, cost)) + geom_line()
costs |> slice_min(cost, n = 3)# A tibble: 3 × 2
t cost
<dbl> <dbl>
1 0.16 0.580
2 0.17 0.580
3 0.15 0.612c(formula = 1 / (1 + 5)) formula
0.1666667 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?
TipReveal answer
The survey did, by choosing costs 5 and 1; the model contributed only \(\hat p\). On these 219 birds it has a lower cost per bird under those costs, and a lower accuracy: “better” has no meaning until a loss is chosen. The formula gives the population optimum for the true \(p(x)\); the curve is an empirical cost on the birds the model was fitted to, with \(\hat p\) plugged in, and its minimum is flat (several thresholds near \(1/6\) tie). So the agreement is between an estimate and the formula, not a proof that the formula is right for new birds. That evaluating on the training rows flatters the model is the validation week’s subject.
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")
# ...
TipReveal solution
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) Truth
Prediction 0 1
0 148 3
1 3 65roc_curve(pg_y, truth, p_flipper, event_level = "second") |> autoplot()
roc_auc(pg_y, truth, p_flipper, event_level = "second")# A tibble: 1 × 3
.metric .estimator .estimate
<chr> <chr> <dbl>
1 roc_auc binary 0.7342. Simulated worlds: adjustment and inverse propensity weighting
Suppose we find a statistically significant correlation between \(X\) and \(Y\). What does it mean?
- A false positive (a spurious correlation)
- \(X\) causes \(Y\)
- \(Y\) causes \(X\)
- 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])
# ...
TipReveal solution
worlds |>
group_by(world) |>
summarise(slope = coef(lm(Y ~ X))[2])# A tibble: 4 × 2
world slope
<chr> <dbl>
1 A 1.98
2 B 1.99
3 C 2.03
4 D -3.79Question 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.
TipReveal answer
A: the average rises by \(2\) (\(Y\)’s equation has \(2X\)). B: no change; \(Y\) depends on \(U\) and its own noise only. C: no change; \(Y\) is upstream of \(X\). D: the average falls by \(3\), the coefficient on \(X\) in \(Y\)’s equation, while the slope in step 1 is about \(-3.8\): \(\mathrm{Cov}(X, Y)/\mathrm{Var}(X) = (2 \cdot (-2) - 3 \cdot 5)/5 = -19/5\), the effect \(-3\) plus the confounding gap \(2 \cdot (-2/5)\).
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 = ____)
# ...
TipReveal solution
intervene <- function(d) {
d |> mutate(
X_do = X + 1,
Y_do = case_when(
world == "A" ~ 2 * X_do + e_y,
world == "B" ~ 2.5 * e_u + e_y,
world == "C" ~ Y,
world == "D" ~ 4 + 2 * Z - 3 * X_do + e_y
)
)
}
worlds_do <- intervene(worlds)
worlds_do |>
group_by(world) |>
summarise(slope = coef(lm(Y ~ X))[2],
effect_of_intervention = mean(Y_do - Y))# A tibble: 4 × 3
world slope effect_of_intervention
<chr> <dbl> <dbl>
1 A 1.98 2
2 B 1.99 0
3 C 2.03 0
4 D -3.79 -3The average change is the average effect of shifting every unit’s \(X\) up by one, computed on these units: exactly \(2\), \(0\), \(0\) and \(-3\). Each world’s equations give every unit the same change, and because the noise is reused nothing else moves.
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 <- ____
# ...
TipReveal solution
d <- worlds |> filter(world == "D")
fit_d <- lm(Y ~ X + Z, data = d)
tidy(fit_d)# A tibble: 3 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 4.27 0.221 19.3 1.56e- 62
2 X -3.04 0.0432 -70.4 9.11e-261
3 Z 1.90 0.0981 19.3 2.02e- 62mean(predict(fit_d, newdata = mutate(d, X = X + 1))) - mean(predict(fit_d, newdata = d))[1] -3.043221Question 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?
TipReveal answer
The two numbers are equal, and both are near \(-3\). The coefficient: the partialling-out theorem (week 2) says the coefficient on X in the long regression is the slope of Y on what is left of X after Z; in world D what is left is e_x, which is independent of everything else in Y, so the population value is \(-3\) and the sample coefficient estimates it. The averaged difference: the adjustment theorem, \(\mathbb E[Y \mid do(X = x)] = \mathbb E_Z[\mathbb E[Y \mid X = x, Z]]\), with the regression as the estimate of the inside. (Step 3 shifts each unit’s own X rather than setting every unit to one value \(x\); in world D the two agree because the effect is the same at every \(x\).) The two numbers coincide exactly because the fitted model is linear in X, so every unit’s difference is the coefficient. Both work because Z holds the only common cause of X and Y, is not a consequence of X, and the noises are independent. In world B the common cause is U, which you would not have: regress Y on X alone and you get the slope of about \(2\) from step 1, while the effect is \(0\). With e_u in the regression the coefficient on X is about \(0\): the theorem working, on a variable you have only because you simulated the world.
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)
# ...
TipReveal solution
fit_pi <- glm(X ~ Z, family = binomial, data = d_e)
d_e <- d_e |> mutate(pi_hat = fitted(fit_pi),
w = ifelse(X == 1, 1 / pi_hat, 1 / (1 - pi_hat)))
c(naive = mean(d_e$Y[d_e$X == 1]) - mean(d_e$Y[d_e$X == 0]),
weighting = mean(d_e$X * d_e$Y / d_e$pi_hat) - mean((1 - d_e$X) * d_e$Y / (1 - d_e$pi_hat)),
regression = coef(lm(Y ~ X + Z, data = d_e))[["X"]]) naive weighting regression
5.185507 1.990335 1.996637 d_e |> slice_max(w, n = 5)# A tibble: 5 × 5
X Y Z pi_hat w
<int> <dbl> <dbl> <dbl> <dbl>
1 1 -3.47 -2.54 0.0161 62.0
2 1 -4.58 -2.01 0.0346 28.9
3 1 -3.61 -1.95 0.0377 26.5
4 1 -1.89 -1.93 0.0387 25.8
5 1 0.245 -1.62 0.0603 16.6Checkpoint, 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) treated untreated
1.0145105 0.9923636 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?
TipReveal answer
The effect is \(2\): \(Y\)’s equation has \(2X\), and nothing else in it depends on \(X\). The naive difference is about \(5.2\): the treated have higher \(Z\), and the confounding gap is \(3\,(\mathbb E[Z \mid X = 1] - \mathbb E[Z \mid X = 0]) \approx 3 \cdot 1.07\). The weighting estimate is centered on \(2\) by the inverse propensity weighting theorem, with \(\hat\pi\) in place of \(\pi\); the logistic model for \(\pi\) is the right one here, because X was generated from a logistic function of Z. The regression coefficient is near \(2\) by the adjustment theorem (or the partialling-out theorem), because \(\mathbb E[Y \mid X, Z]\) is linear in X and Z. Weighting needs a correct regression model for the treatment; regression adjustment needs a correct regression model for the outcome; both need Z to hold the common causes. Of the three, the weighting estimate varies most from run to run: in about one run in ten it falls outside \(1.7\) to \(2.3\), so compare yours with your neighbors’.
The largest weights belong to treated units with low Z or untreated units with high Z: units whose treatment was unlikely given their Z. Their \(\hat\pi\) (or \(1 - \hat\pi\)) is a few percent, so each stands in for tens of units. The largest weight is usually between about 15 and 100, and it changes a lot from run to run. With the other weights held fixed (not refitting \(\hat\pi\)), write the estimate as an average, \(\hat\tau = \frac1n \sum_i s_i\) with \(s_i = (2X_i - 1)\, w_i Y_i\): a treated unit adds \(w_i Y_i\), an untreated one subtracts it. Leaving unit \(i\) out gives \((n\hat\tau - s_i)/(n - 1)\), so the estimate moves by \((\hat\tau - s_i)/(n - 1)\), in either arm. For the unit with the largest weight in your sample:
s <- (2 * d_e$X - 1) * d_e$w * d_e$Y
i <- which.max(d_e$w)
c(estimate = mean(s), largest_weight = d_e$w[i], change_if_left_out = (mean(s) - s[i]) / (nrow(d_e) - 1)) estimate largest_weight change_if_left_out
1.9903346 61.9683026 0.1086913 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)
# ...