9  Using the propensity score

Work-in-progress 🚧

You are reading the work-in-progress first edition of Causal Inference in R. This chapter has its foundations written but is still undergoing changes.

The propensity score is a balancing tool – we use it to help us make our exposure groups exchangeable. There are many ways to incorporate the propensity score into an analysis. Commonly used techniques include stratification (estimating the causal effect within propensity score stratum), matching, weighting, and direct covariate adjustment. In this section, we will focus on matching and weighting; other techniques will be discussed once we introduce the outcome model. Recall at this point in the book we are still in the design phase. We have not yet incorporated the outcome into our analysis at all.

9.1 Matching

Ultimately, we want the exposed and unexposed observations to be exchangeable with respect to the confounders we have proposed in our DAG (so we can use the observed effect for one to estimate the counterfactual for the other). One way to do this is to ensure that each observation in our analysis sample has at least one observation of the opposite exposure that has matching values for each of these confounders. If we had a small number of binary confounders, for example, we might be able to construct an exact match for observations (and only include those for whom such a match exists), but as the number and continuity of confounders increases, exact matching becomes less feasible. This is where the propensity score, a summary measure of all of the confounders, comes in to play.

Let’s setup the data as we did in Chapter 8.

library(broom)
library(touringplans)

seven_dwarfs_9 <- seven_dwarfs_train_2018 |> filter(wait_hour == 9)

We can re-fit the propensity score using the MatchIt package, as below. Notice here the matchit function fit a logistic regression model for our propensity score, as we had in Chapter 8. There were 60 days in 2018 where the Magic Kingdom had extra magic morning hours. For each of these 60 exposed days, matchit found a comparable unexposed day, by implementing a nearest-neighbor match using the constructed propensity score. Examining the output, we also see that the target estimand is an “ATT” (do not worry about this yet, we will discuss this and several other estimands in Chapter 11).

library(MatchIt)
m <- matchit(
  park_extra_magic_morning ~ park_ticket_season + park_close + park_temperature_high,
  data = seven_dwarfs_9
)
m
A matchit object
 - method: 1:1 nearest neighbor matching without replacement
 - distance: Propensity score
             - estimated with logistic regression
 - number of obs.: 354 (original), 120 (matched)
 - target estimand: ATT
 - covariates: park_ticket_season, park_close, park_temperature_high

We can use the get_matches function to create a data frame with the original variables that only consists of those who were matched. Notice here our sample size has been reduced from the original 354 days to 120.

matched_data <- get_matches(m)
glimpse(matched_data)
Rows: 120
Columns: 18
$ id                       <chr> "5", "340", "12", "1…
$ subclass                 <fct> 1, 1, 2, 2, 3, 3, 4,…
$ weights                  <dbl> 1, 1, 1, 1, 1, 1, 1,…
$ park_date                <date> 2018-01-05, 2018-12…
$ wait_hour                <int> 9, 9, 9, 9, 9, 9, 9,…
$ attraction_name          <chr> "Seven Dwarfs Mine T…
$ wait_minutes_actual_avg  <dbl> 33.0, 8.0, 114.0, 32…
$ wait_minutes_posted_avg  <dbl> 70.56, 80.62, 79.29,…
$ attraction_park          <chr> "Magic Kingdom", "Ma…
$ attraction_land          <chr> "Fantasyland", "Fant…
$ park_open                <time> 09:00:00, 09:00:00,…
$ park_close               <time> 24:00:00, 23:00:00,…
$ park_extra_magic_morning <dbl> 1, 0, 1, 0, 1, 0, 1,…
$ park_extra_magic_evening <dbl> 0, 0, 0, 0, 0, 0, 0,…
$ park_ticket_season       <chr> "regular", "regular"…
$ park_temperature_average <dbl> 43.56, 57.61, 70.91,…
$ park_temperature_high    <dbl> 54.29, 65.44, 78.26,…
$ distance                 <dbl> 0.18410, 0.18381, 0.…

9.2 Weighting

One way to think about matching is as a crude “weight” where everyone who was matched gets a weight of 1 and everyone who was not matched gets a weight of 0 in the final sample. Another option is to allow this weight to be smooth, applying a weight to allow, on average, the covariates of interest to be balanced in the weighted population. To do this, we will construct a weight using the propensity score. There are many different weights that can be applied, depending on your target estimand of interest (see Chapter 11 for details). For this section, we will focus on the “Average Treatment Effect” weights, commonly referred to as an “inverse probability weight”. The weight is constructed as follows, where each observation is weighted by the inverse of the probability of receiving the exposure they received.

\[w_{ATE} = \frac{X}{p} + \frac{(1 - X)}{1 - p}\]

For example, if observation 1 had a very high likelihood of being exposed given their pre-exposure covariates (\(p = 0.9\)), but they in fact were not exposed, their weight would be 10 (\(w_1 = 1 / (1 - 0.9)\)). Likewise, if observation 2 had a very high likelihood of being exposed given their pre-exposure covariates (\(p = 0.9\)), and they were exposed, their weight would be 1.1 (\(w_2 = 1 / 0.9\)). Intuitively, we give more weight to observations who, based on their measured confounders, appear to have useful information for constructing a counterfactual – we would have predicted that they were exposed and but by chance they were not, or vice-versa. The propensity package is useful for implementing propensity score weighting.

library(propensity)

seven_dwarfs_9_with_ps <-
  glm(
    park_extra_magic_morning ~ park_ticket_season + park_close + park_temperature_high,
    data = seven_dwarfs_9,
    family = binomial()
  ) |>
  augment(type.predict = "response", data = seven_dwarfs_9)
seven_dwarfs_9_with_wt <- seven_dwarfs_9_with_ps |>
  mutate(w_ate = wt_ate(.fitted, park_extra_magic_morning))

Table 9.1 shows the weights in the first column.

seven_dwarfs_9_with_wt |>
  select(
    w_ate,
    .fitted,
    park_date,
    park_extra_magic_morning,
    park_ticket_season,
    park_close,
    park_temperature_high
  ) |>
  head() |>
  knitr::kable()
w_ate .fitted park_date park_extra_magic_morning park_ticket_season park_close park_temperature_high
1.433 0.3019 2018-01-01 0 peak 23:00:00 58.63
1.392 0.2815 2018-01-02 0 peak 24:00:00 53.65
1.409 0.2900 2018-01-03 0 peak 24:00:00 51.11
1.232 0.1881 2018-01-04 0 regular 24:00:00 52.66
5.432 0.1841 2018-01-05 1 regular 24:00:00 54.29
1.262 0.2074 2018-01-06 0 regular 23:00:00 56.25
Table 9.1: The first six observations in the seven_dwarfs_9_with_wt dataset, including their propensity scores in the .fitted column and weight in the w_ate column.