Observational Causal Inference

The Essential Approach

Learn to predict outcomes for all combinations of treatment and covariates.
Compare predictions for different treatments to estimate the effect of treatment.

Average these comparisons to estimate a summary of the effect.

Step 1: Learning to Predict

formula_options = new Map(
  [["lines","y ~ w*(x.1 + x.2 + x.3 + x.4 + x.5 + x.6)"],
   ["parallel lines","y ~ w + x.1 + x.2 + x.3 + x.4 + x.5 + x.6"]])
viewof example_formula = Inputs.radio(formula_options, 
  {value: formula_options.get("lines"), label: "model"})

formula_options = Map(2) {"lines" => "y ~ w*(x.1 + x.2 + x.3 + x.4 + x.5 + x.6)", "parallel lines" => "y ~ w + x.1 + x.2 + x.3 + x.4 + x.5 + x.6"}

model

linesparallel lines

example_formula = "y ~ w*(x.1 + x.2 + x.3 + x.4 + x.5 + x.6)"

R Code

Start Over

Run Code

›

linear.predictor = lm(example_formula, data=sam)

muhat.WX = predict(linear.predictor, newdata=sam)

_webr_editor_16 = Object {code: "linear.predictor = lm(example_formula, data=sam)\nmuhat.WX = predict(linear.predictor, newdata=sam) ", options: Object, indicator: Ke}

The Predictions
The Data

R Code

Start Over

Run Code

›

ggplot() + scatter(sam |> mutate(muhat=muhat.WX), xx, 'muhat',

shape=square.shape, jitter=no.jitter)

_webr_editor_17 = Object {code: "ggplot() + scatter(sam |> mutate(muhat=muhat.WX), … shape=square.shape, jitter=no.jitter)", options: Object, indicator: Ke}

R Code

Start Over

Run Code

›

ggplot() + scatter(sam, xx, yy)

_webr_editor_18 = Object {code: "\nggplot() + scatter(sam, xx, yy) ", options: Object, indicator: Ke}

You’re looking at predictions based a higher-dimensional version of the lines model.⁴
- The little squares you see are the predictions we get from the model.
- To compare the predictions to the actual data, flip over the the ‘Data’ tab.
Everything is plotted against the covariate you chose on Slide 0.6.⁵

R Code

Start Over

Run Code

›

example_formula

_webr_editor_19 = Object {code: "example_formula", options: Object, indicator: Ke}

Step 2: Comparing Predictions

R Code

Start Over

Run Code

›

muhat.1 = predict(linear.predictor, newdata=sam |> mutate(w=1))

muhat.0 = predict(linear.predictor, newdata=sam |> mutate(w=0))

tau = muhat.1 - muhat.0

tau.plot = ggplot() + scatter(sam |> mutate(tau=tau), xx, 'tau',

shape=triangle.shape, jitter=no.jitter)

tau.plot

_webr_editor_20 = Object {code: "muhat.1 = predict(linear.predictor, newdata=sam |>… shape=triangle.shape, jitter=no.jitter)\ntau.plot", options: Object, indicator: Ke}

Here I’ve predicted the treatment effect for each person in the sample.
And I’ve plotted them, as little triangles, against your chosen covariate.

Step 3: Averaging the Comparisons

R Code

Start Over

Run Code

›

Deltahat.all = mean(tau)

Deltahat.1 = mean(tau[sam$w==1])

tau.plot + geom_hline(yintercept=Deltahat.all, color='purple') +

geom_hline(yintercept=Deltahat.1, color=groupb, linetype='dashed')

_webr_editor_21 = Object {code: "Deltahat.all = mean(tau)\nDeltahat.1 = mean(tau[s…pt=Deltahat.1, color=groupb, linetype='dashed')", options: Object, indicator: Ke}

Here I’ve calculated two averages of this predicted effect.
- is the average over all people in the sample.
- is the average over the treated people in the sample.
I’ve shown them as horizontal lines on the plot of above.
Here are the numbers.

R Code

Start Over

Run Code

›

data.frame(Deltahat.all, Deltahat.1)

_webr_editor_22 = Object {code: "data.frame(Deltahat.all, Deltahat.1)", options: Object, indicator: Ke}

Figure 1

Model Choice

If you haven’t touched anything, what you’re seeing is based on the lines model.
Go back to Step 1 and click the button for the parallel lines model.
- Everything else will update automatically.
- And the results will be different. What’s going on?

Questions
Answers

What’s different in Step 2?
What does the mean for Step 3?

In Step 2, will be the same for everybody.
- For every function in the parallel lines model, is the same.
In Step 3, as a result, and will be the same.
- If you average a constant function over two different samples, you get the same thing for both.
- You get that constant.

Shortcuts

Letting R Do the Work

With certain models, we can skip steps.
- e.g., if we predict a constant effect of treatment, there’s no need to average.
- This happens when we use certain regression models.
If is in , then

And if we parameterize it the right way, is just a coefficient.
And our prediction for the treatment effect is just a coefficient.
This means that if we get R to tell us the coefficient , we’re done.

R Code

Start Over

Run Code

›

linear.predictor.parallel = lm(y ~ w + x.1 + x.2 + x.3 + x.4 + x.5 + x.6, data=sam)

summary(linear.predictor.parallel)

_webr_editor_23 = Object {code: "linear.predictor.parallel = lm(y ~ w + x.1 + x.2 +…x.6, data=sam)\nsummary(linear.predictor.parallel)", options: Object, indicator: Ke}

Take a look at the coefficient on in this summary output.
Look back at Figure 1 to compare it to the estimates and we calculated in the previous section.
If you used the parallel lines model there, this coefficient should be the same as both of those.

A Problem: We Predict Frequencies Outside

R Code

Start Over

Run Code

›

muhat.WX = predict(linear.predictor.parallel, newdata=sam)

ggplot() + scatter(sam |> mutate(muhat=muhat.WX), xx, 'muhat',

shape=square.shape, jitter=no.jitter)

_webr_editor_24 = Object {code: "muhat.WX = predict(linear.predictor.parallel, newd… shape=square.shape, jitter=no.jitter)", options: Object, indicator: Ke}

The little squares you see are the predictions .
Some of them are outside the interval . In particular, a couple are bigger than .
This makes people uncomfortable. To fix this, they often use what are called generalized linear models.
Here is a function that maps the real line into the interval : .

R Code

Start Over

Run Code

›

invlogit = function(x) { 1 / (1 + exp(-x)) }

xs = seq(-5, 5, length=1000)

ggplot() + geom_line(aes(x=xs, y=pnorm(xs)), color='red') +

geom_line(aes(x=xs, y=invlogit(xs)), color='blue') +

labs(y='inverse logit(x)')

_webr_editor_25 = Object {code: "\ninvlogit = function(x) { 1 / (1 + exp(-x)) }\nxs =…(xs)), color='blue') +\nlabs(y='inverse logit(x)')", options: Object, indicator: Ke}

Sometimes people use other functions that do the same thing.
- e.g. cumulative distribution function of the standard normal. They call that a ‘probit’ model.
- These two curves have about the same shape, so it usually doesn’t matter that much which which you use.

Problem Solved?

Step 1: Learning to Predict
Step 2: Comparing Predictions
Step 3: Averaging the Comparisons

R Code

Start Over

Run Code

›

logit.predictor = glm(y ~ w + x.1 + x.2 + x.3 + x.4 + x.5 + x.6,

family='binomial', data=sam)

muhat = predict(logit.predictor, newdata=sam, type='response')

_webr_editor_26 = Object {code: "logit.predictor = glm(y ~ w + x.1 + x.2 + x.3 + x.…t(logit.predictor, newdata=sam, type='response') ", options: Object, indicator: Ke}

R Code

Start Over

Run Code

›

ggplot() + scatter(sam |> mutate(muhat=muhat), xx, 'muhat',

shape=square.shape, jitter=no.jitter)

_webr_editor_27 = Object {code: "ggplot() + scatter(sam |> mutate(muhat=muhat), xx,… shape=square.shape, jitter=no.jitter)", options: Object, indicator: Ke}

Here are the predictions we get using the logistic regression version of the parallel lines model.
- Now all the predictions are in the interval . That’s how logistic regression works.
- No matter what you put into , you get something in .

We can use our predictor to estimate the ATE and ATT just like we did before.

R Code

Start Over

Run Code

›

muhat1 = predict(logit.predictor, newdata=sam |> mutate(w=1), type='response')

muhat0 = predict(logit.predictor, newdata=sam |> mutate(w=0), type='response')

tauX = muhat1 - muhat0

tau.plot = ggplot() + scatter(sam |> mutate(tau=tauX), xx, 'tau',

shape=triangle.shape, jitter=no.jitter)

tau.plot

_webr_editor_28 = Object {code: "muhat1 = predict(logit.predictor, newdata=sam |> m… shape=triangle.shape, jitter=no.jitter)\ntau.plot", options: Object, indicator: Ke}

Here I’ve predicted the treatment effect for each person in the sample.
And I’ve plotted them, as little triangles, against your chosen covariate.

R Code

Start Over

Run Code

›

Deltahat.all = mean(tauX)

Deltahat.1 = mean(tauX[sam$w==1])

tau.plot +

geom_hline(yintercept=Deltahat.all, color='purple') +

geom_hline(yintercept=Deltahat.1, color=groupb, linetype='dashed')

_webr_editor_29 = Object {code: "Deltahat.all = mean(tauX)\nDeltahat.1 = mean(tauX…pt=Deltahat.1, color=groupb, linetype='dashed')", options: Object, indicator: Ke}

Here I’ve calculated two averages of this predicted effect.
- is the average over all people in the sample.
- is the average over the treated people in the sample.
I’ve shown them as horizontal lines on the plot of above.
Here are the numbers. Go ahead and compare them to the ones in Figure 1, when we used a linear model.

R Code

Start Over

Run Code

›

data.frame(Deltahat.all, Deltahat.1)

_webr_editor_30 = Object {code: "data.frame(Deltahat.all, Deltahat.1)", options: Object, indicator: Ke}

The Fix and the Shortcut

R Code

Start Over

Run Code

›

Deltahat.all = mean(tauX)

Deltahat.1 = mean(tauX[sam$w==1])

tau.plot +

geom_hline(yintercept=Deltahat.all, color='purple') +

geom_hline(yintercept=Deltahat.1, color=groupb, linetype='dashed')

_webr_editor_31 = Object {code: "Deltahat.all = mean(tauX)\nDeltahat.1 = mean(tauX…pt=Deltahat.1, color=groupb, linetype='dashed')", options: Object, indicator: Ke}

Now what happens when we use our shortcut?
- One thing you’ll notice is that we’re not predicting the same treatment effect for everyone.
- So the ATT and ATE are different. And the coefficient can’t be both.
What is the coefficient in our predictor ?

Here’s what R says.

R Code

Start Over

Run Code

›

summary(logit.predictor)

_webr_editor_32 = Object {code: "summary(logit.predictor)", options: Object, indicator: Ke}

This — or something like it — is what people usually report.
- Why? Because it’s what their software (usually R) gives them.
- It’s totally different from our estimates of the ATE and ATT.
- It’s an estimate of the log odds-ratio of getting PCI with those two treatments.

It’s a bit of a calculation to show this, but all you’ve got to do is plug in your formula for and do some algebra. We did it in class.

Interpreting Odds Ratios

We should know a bit about how to interpret odds ratios because people report them.
- Shah et al. (2007) reports an odds ratio point estimate and a 95% confidence interval .⁶
What the odds ratio that’s meant to be in this interval is a bit unclear.
- We’d expect it to vary with so there isn’t any one odds ratio.
But let’s suppose the odds ratio were actually constant and try to work out what that’d mean.
A starting point is to observe that odds are increasing as a function of frequency.

R Code

Start Over

Run Code

›

p = seq(.01, .99, by=.01)

ggplot() + geom_line(aes(x=p, y=p/(1-p)), color='blue') +

labs(x='p', y='p/(1-p)')

_webr_editor_33 = Object {code: "\np = seq(.01, .99, by=.01)\nggplot() + geom_line(ae…-p)), color='blue') + \n labs(x='p', y='p/(1-p)')", options: Object, indicator: Ke}

That means it’s equivalent to say that …
- the frequency is higher than
- the odds are higher than .
- and the odds ratio is bigger than 1.
This is basically all people tend to know about odds ratios.
- If , we’re predicting that PCI happens more with treatment than without treatment.
- If , we’re predicting that PCI happens less with treatment than without treatment.
Beyond that, it’s hard to interpret. Let’s look at a picture.
Here’s the odds ratio plotted against two frequencies .

The Estimates in Shah et al. (2007)

In the data analyzed in Shah et al. (2007), for-profit hospitals do PCI less often than non-profit hospitals.
- 46.5% of the time vs. 47.1% of the time in non-profit hospitals.
- If we really want this as an odds ratio, we can calculate it. It’s a little less than 1.
But that’s a raw difference. We know those are impacted by covariate shift. e.g …
- differences in the patients who end up in for-profit vs. non-profit hospitals.
- hospital characteristics other than for/non-profit status: size, region, etc.
They report an adjusted odds ratio estimate of 1.14. It’s bigger than 1. That tells the opposite story.
- It says similar patients at similar for and non-profit hospitals get PCI more often at the for-profit ones.
- It’s Simpson’s paradox, just like in our Berkeley admissions example.
But 1.14 is a point estimate. Their 95% confidence interval includes numbers less than 1, too.
- So, in essence, their study is inconclusive about PCI and for-profit hospitals. That’s what they say.
- But inconclusiveness is a matter of degree. If you want to get in some practice thinking about inconclusiveness, odds ratios, and how calibration works, take a look here in the Appendix.

What About Misspecification?

If we’re trying to think of all this as meaningful, we’re making some pretty strong assumptions.
To start with, we haven’t been clear about what odds ratio we’re trying to estimate, exactly.
Even if the odds ratio were constant — i.e. if the ‘unless’ part was true — misspecification is likely.
- We get constancy if the log odds are additive in and .
- Our model includes a small subset of functions like that: the ones where
  the log odds are linear combinations of .
We’re doing something to adjust for covariate shift, we may not be adjusting for it successfully.

Misspecification and Simulated Data

Let’s think through what happens with the simulated data we’ve been using.
- In that data, the odds ratio is constant. It’s 1. Treatment does nothing at all.
- for all . See Equation 1.
But the estimate we get using our logistic regression model isn’t all that close to one.
- We can base a normal-approximation interval for its log on the summary output from Slide 2.4.
- Here’s our interval. We’re confidently wrong (Remember that ).

R Code

Start Over

Run Code

›

point.estimate = coef(logit.predictor)['w']

its.standard.deviation = sqrt(vcov(logit.predictor)['w','w'])

point.estimate + c(-1,1) * 1.96 * its.standard.deviation

_webr_editor_34 = Object {code: "point.estimate = coef(logit.predictor)['w']\nits.st…timate + c(-1,1) * 1.96 * its.standard.deviation ", options: Object, indicator: Ke}

What We Learn from This

We’re doing essentially the same thing they did in Shah et al. (2007). It’s not working.
- Our log-odds-ratio target is zero but our interval estimate nowhere near there.
- That’s worrisome. Perhaps we shouldn’t trust the concusions of Shah et al. (2007).
We see the same problem in the summary output for our ‘parallel lines model’ on Slide 2.1.

R Code

Start Over

Run Code

›

point.estimate = coef(linear.predictor.parallel)['w']

its.standard.deviation = sqrt(vcov(linear.predictor.parallel)['w','w'])

point.estimate + c(-1,1) * 1.96 * its.standard.deviation

_webr_editor_35 = Object {code: "point.estimate = coef(linear.predictor.parallel)['…timate + c(-1,1) * 1.96 * its.standard.deviation ", options: Object, indicator: Ke}

Our target in that context—whether we’re thinking of the ATE or ATT—is zero.
- Our interval is tiny but everything in it is way off. Way bigger than zero. This is bad.
- It’s bad enough that our estimate is far off target.
- What makes it worse is that we’re saying that it’s very close (because our interval is so small).
We get similar point estimates when we calculate the ATE and ATT ‘the long way’
using the other models we’ve looked at.
- The ‘lines’ model (Slide 1.4).
- The ‘logit version’ of the ‘parallel lines model’ (Slide 2.3).
All of these estimators are badly biased.
- We’ll take a look at their sampling distributions when we discuss Kim et al. (2023).
- It seems like we should be doing something else. Let’s get started on that.

Lecture 15

Summary

The Heart-Attack Context

Assumptions: Randomization and Random Sampling

About the Data

Fake Data from Kim et al. (2023)

What the Data Looks Like

Observational Causal Inference

The Essential Approach

Step 1: Learning to Predict

Step 2: Comparing Predictions

Step 3: Averaging the Comparisons

Model Choice

Shortcuts

Letting R Do the Work

A Problem: We Predict Frequencies Outside

Problem Solved?

The Fix and the Shortcut

Interpreting Odds Ratios

The Estimates in Shah et al. (2007)

What About Misspecification?

Misspecification and Simulated Data

What We Learn from This

Using Weighting to Reduce Bias

Balanced Covariate Distributions

Stable Balancing Weights

What to Read

Shah et al. (2007)

Kim et al. (2023)

Appendix

About Inconclusiveness

Footnotes