Definition:

The propensity score is typically defined as the probability of receiving some program (or “treatment”) versus a comparison condition, given a set of observed baseline covariates. (Rosenbaum and Rubin, 1983).

We use \(e_i\) to denote propensity score, the definition can be written as : \(e_i = Pr(Z_i|X_i)\), where \(X_i = (X_1,X_2,…,X_n)\) is the set of subject \(i\)’s observed baseline covariates, and \(Z_i\) is an indicator variable denoting subject \(i\)’s treatment received (\(Z_i=1\) for active treatment and \(Z_i=0\) for control treatment).

Motivation:

Because the proposition of propropensity score is to perfect the Rubin Causal Model (Rubin, 1974) which is described as potential outcome framework, I begin this part by describing the conceptual framework. In this process, I will exhibit the advantages on the propensity score.

Here I refer to Dr. Peng Ding’s blog on “Capital of Statistics”. In this blog, Dr.Ding gives a detailed explanation on the propensity score.

Let’s denote \(\{Y_i(1),Y_i(0)\}\) as subject \(i\)’s potential outcome under the active treatment and control treatment. Then we can use \(Y_i(1)-Y_i(0)\) to denote subject \(i\)’s causal effect (treatment effect). Unfortunately, for individual patient,the treatment received is either active or control. The observed outcome only can be \(Y_i = Z_i*Y_i(1)+(1-Z_i)*Y_i(0)\). Individual patient’s causal effect can’t be observed. However, if \(Z_i\) is completely randomized, we can observe the whole sample subjects’ average causal effect (ACE):

\[ACE(Z\to Y) = E\{Y_i(1)-Y_i(0)\}\] Since \[
\begin{aligned}
ACE(Z\to Y) &= E\{Y_i(1)\} – E\{Y_i(0)\} \\
&= E\{Y_i(1)|Z_i=1\} – E\{Y_i(0)|Z_i=0\}\\
&= E\{Y_i|Z_i=1\} – E\{Y_i|Z_i=0\}
\end{aligned}
\] based on the assumption: \[Z\perp\{Y_i(1),Y_i(0)\}\] Nonetheless, in some studies, we prior know some covariates have strong relationships with the result. We control these covariates to reduce the variance in experiments. Then we find \(Z\perp\{Y_i(1),Y_i(0)\}\) can’t be satisfied. More generally, \[Z\perp\{Y_i(1),Y_i(0)|X\}\] Then, \[
\begin{aligned}
ACE(Z\to Y) &= E\{Y(1)\} – E\{Y(0)\} \\
&= E\{Y(1)|X\} – E\{Y(0)|X\} \\
&= E\{Y(1)|X,Z=1\} – E\{Y(0)|X,Z=0\}\\
&= E\{Y|X,Z=1\} – E\{Y|X,Z=0\}
\end{aligned}
\] From above formular, to get the average causal effect (ACE), we need to estimate \(E\{Y|X,Z=z\}(z=0,1)\). It becames a regression problem. One assumption is to use simple linear model: \(E\{Y|X,Z\}=\beta_0+\beta_xX+\beta_zZ\) to fit the data, and here \(\beta_z\) is the average causal effect (ACE). Although the linear model is easy to be realized, it has obvious defects.

The core idea of Rubin Causal model is “matching”. More specifically, we want to find “similar” subjects in different groups by “matching”. By comparing the responses of these “similar” subjects, we estimate the average causal effect (ACE). However, as dimensions of covariates increasing, using linear model to realize “matching” becomes hard. To perfect Rubin Causal model, in 1983, Paul Rosenbaum and Donald Rubin defined the propensity score \[e(X) = Pr(Z=1|X)\] The propensity score is a balancing score: conditional on the propensity score, the distribution of measured baseline covariates is similar between treated and untreated subjects. Thus, in a set of subjects all of whom have the same propensity score, the distribution of observed baseline covariates will be the same between the treated and untreated subjects.

It satisfies \(Z\perp X|e(X)\) (also called sufficient dimension reduction), and \(Z\perp \{Y_i(1),Y_i(0)\}|e(X)\). Hence, we can rewrite the formular as: \[
\begin{aligned}
ACE &= E\{Y|e(X),Z=1\} – E\{Y|e(X),Z=0\}\\
&= E\{Y(1)\}-E\{Y(0)\}\\
&= E\{\frac{ZY}{e(X)}\}-E\{\frac{(1-Z)Y}{1-e(X)}\}
\end{aligned}
\].

How to use R to estimate propensity score:

Here we ues logistics regression to estimate the propensity score

data("lalonde", package = "cobalt")   # use dataset "lalonde" in package"cobalt"
head(lalonde)     # first 6 rows of dataset "lalonde"

##   treat age educ   race married nodegree re74 re75       re78
## 1     1  37   11  black       1        1    0    0  9930.0460
## 2     1  22    9 hispan       0        1    0    0  3595.8940
## 3     1  30   12  black       0        0    0    0 24909.4500
## 4     1  27   11  black       0        1    0    0  7506.1460
## 5     1  33    8  black       0        1    0    0   289.7899
## 6     1  22    9  black       0        1    0    0  4056.4940

m_ps <- glm(treat ~ age + treat + educ + race + nodegree + married +  re74 + re75, family = binomial(), data = lalonde) # use logistics regression model for propensity score estimation

## Warning in model.matrix.default(mt, mf, contrasts): the response appeared
## on the right-hand side and was dropped

## Warning in model.matrix.default(mt, mf, contrasts): problem with term 2 in
## model.matrix: no columns are assigned

summary(m_ps)

## 
## Call:
## glm(formula = treat ~ age + treat + educ + race + nodegree + 
##     married + re74 + re75, family = binomial(), data = lalonde)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7645  -0.4736  -0.2862   0.7508   2.7169  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.663e+00  9.709e-01  -1.713  0.08668 .  
## age          1.578e-02  1.358e-02   1.162  0.24521    
## educ         1.613e-01  6.513e-02   2.477  0.01325 *  
## racehispan  -2.082e+00  3.672e-01  -5.669 1.44e-08 ***
## racewhite   -3.065e+00  2.865e-01 -10.699  < 2e-16 ***
## nodegree     7.073e-01  3.377e-01   2.095  0.03620 *  
## married     -8.321e-01  2.903e-01  -2.866  0.00415 ** 
## re74        -7.178e-05  2.875e-05  -2.497  0.01253 *  
## re75         5.345e-05  4.635e-05   1.153  0.24884    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 751.49  on 613  degrees of freedom
## Residual deviance: 487.84  on 605  degrees of freedom
## AIC: 505.84
## 
## Number of Fisher Scoring iterations: 5

prs_df <- data.frame(pr_score = predict(m_ps, type = "response"),
                     treat = m_ps$model$treat)
head(prs_df)

##    pr_score treat
## 1 0.6387699     1
## 2 0.2246342     1
## 3 0.6782439     1
## 4 0.7763241     1
## 5 0.7016387     1
## 6 0.6990699     1

Propensity score matching (PSM) method

Propensity score matching (PSM) method is one-to-one or pair matching, in which pairs of treated and untreated subjects are formed, such that matched subjects have similar values of the propensity score.

Here we use “MatchIt” package for PSM analysis

library(MatchIt)

## 
## Attaching package: 'MatchIt'

## The following object is masked _by_ '.GlobalEnv':
## 
##     lalonde

match_model <- matchit(treat ~ age + educ + race + nodegree + married + re74 + re75, data = lalonde, method = "nearest")
# summary(match_model)
match_data <- match.data(match_model)
library(dplyr) 

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

match_data <- match_data %>% rename(pr_score = distance)
head(match_data %>%
       select(pr_score,treat))

##    pr_score treat
## 1 0.6387699     1
## 2 0.2246342     1
## 3 0.6782439     1
## 4 0.7763241     1
## 5 0.7016387     1
## 6 0.6990699     1

Now we estimate the average causal effect (ACE) after we doing PSM analysis

match_data %>%
  group_by(treat) %>%
  summarise(Income1978 = mean(re78),
            n = n())

## # A tibble: 2 x 3
##   treat Income1978     n
##   <int>      <dbl> <int>
## 1     0      5441.   185
## 2     1      6349.   185

Propensity score weighting (IPTW) method

Since the PSM method is one-to-one or pair matching, the control group would be shrunk down to the same size as the treatment group. A possibly better use of the propensity scores is to keep all observations in play but weight them according to the propensity score.

Weights can be defined as \(w_i = \frac{Z_i}{e(X_i)}-\frac{(1-Z_i)}{1-e(X_i)}\). A subject’s weight is equal to the inverse of the probability of receiving the treatment that the subject actually received. Inverse probability of treatment weighting was first proposed by Rosenbaum (1987a) as a form of model-based direct standardization. \[
\begin{aligned}
ACE &= E\{\frac{ZY}{e(X)}\}-E\{\frac{(1-Z)Y}{1-e(X)}\}\\
&= E\{w_iY_i(1)\}-E\{w_iY_i(0)\}
\end{aligned}
\].

Here we use “WeightIt” package for IPTW analysis

library("WeightIt")
library("cobalt")

##  cobalt (Version 3.9.0, Build Date: 2019-10-06 19:10:02 UTC)
##    Please read the documentation at ?bal.tab to understand the default outputs.
##    Submit bug reports and feature requests to https://github.com/ngreifer/cobalt/issues
##    Install the development version (not guaranteed to be stable) with:
##      devtools::install_github("ngreifer/cobalt")
##    Thank you for using cobalt!

## 
## Attaching package: 'cobalt'

## The following object is masked from 'package:MatchIt':
## 
##     lalonde

W.out <- weightit(treat ~ age + educ + race + married + nodegree + re74 + re75,
        data = lalonde, estimand = "ATT", method = "ps")

lalonde %>%
  mutate(weights = get.w(W.out)) %>%
  group_by(treat) %>%
  summarise(Income1978_weighted = weighted.mean(re78, weights),
            n = n())

## # A tibble: 2 x 3
##   treat Income1978_weighted     n
##   <int>               <dbl> <int>
## 1     0               5135.   429
## 2     1               6349.   185