Propensity Score-Based Analysis Using R
Lab for Assignment 2
Xu Qin
A comprehensive list of propensity score software can be found here: https://www.biostat.jhsph.edu/~estuart/propensityscoresoftware.html
We will use Lalonde’s data on the evaluation of the National Supported Work program. The statistical quantity of interest is the causal effect of the treatment (treat) on 1978 earnings (re78). The other variables are pre-treatment covariates. See ?lalonde for more information on this dataset.
install.packages("MatchIt")# If you have already installed it, then you do not need to reinstall it.library("MatchIt")
data(lalonde)
head(lalonde)## treat age educ race married nodegree re74 re75 re78
## NSW1 1 37 11 black 1 1 0 0 9930.0460
## NSW2 1 22 9 hispan 0 1 0 0 3595.8940
## NSW3 1 30 12 black 0 0 0 0 24909.4500
## NSW4 1 27 11 black 0 1 0 0 7506.1460
## NSW5 1 33 8 black 0 1 0 0 289.7899
## NSW6 1 22 9 black 0 1 0 0 4056.49401. Logistic regression in R
To increase the efficiency of a propensity score-based method, it is suggested to include covariates in the propensity score model that are related to the outcome variable, no matter if they are significantly related to the treatment or not. Theories show that age, number of years of schooling, race, married status, high school degree, and income before the treatment are all related to the the income after the treatment. We further check the association between each covariate with the outcome based on the sample data. Below we use age as an illustration.
summary(lm(re78 ~ age, data = lalonde))##
## Call:
## lm(formula = re78 ~ age, data = lalonde)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9013 -6041 -1826 4024 53464
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4594.75 884.01 5.198 2.75e-07 ***
## age 80.33 30.39 2.643 0.00842 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7435 on 612 degrees of freedom
## Multiple R-squared: 0.01129, Adjusted R-squared: 0.009673
## F-statistic: 6.988 on 1 and 612 DF, p-value: 0.008418Similarly, we found that in addition to age, educ, race, married, nodegree, re74, and re75 are all significantly associated with the outcome re78. Therefore, we include all of them in the following propensity score model.
l = glm(treat ~ age + educ + race + married + nodegree + re74 + re75, data = lalonde, family = binomial(link = "logit")) # run the logistic regression
summary(l) # extract the output of the regression##
## Call:
## glm(formula = treat ~ age + educ + race + married + nodegree +
## re74 + re75, family = binomial(link = "logit"), data = lalonde)
##
## 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 ***
## married -8.321e-01 2.903e-01 -2.866 0.00415 **
## nodegree 7.073e-01 3.377e-01 2.095 0.03620 *
## 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: 52. Overlapping in logit of propensity score
2.1 Calculate the logit of propensity score
logit.p = predict(l, type = "link") # save the logit propensity score2.2 Compare the distribution of the logit score between the treatment group and the control group
install.packages("ggplot2") # If you have already installed it, then you do not need to reinstall it.library(ggplot2)
lalonde = cbind(lalonde, logit.p)
ggplot(lalonde, aes(x = logit.p, fill = as.factor(treat))) +
geom_density(alpha = .5) +
labs(x = "Logit Propensity Scores", y = "Density", fill = "") +
theme_bw()
The plot indicates that the distribution of the logit propensity score in the treatment group well overlaps with that in the control group.
A good overlap indicates that for many levels of propensity scores, you have both treated and untreated subjects. This means that for an individual in one group, there are observations with the same baseline characteristics (captured by propensity scores) in the other group to provide counterfactual information for the individual.
3. Matching
A step-by-step guidance on the implementation of MatchIt: https://cran.r-project.org/web/packages/MatchIt/vignettes/MatchIt.html
A more comprehensive introduction of the package: https://imai.fas.harvard.edu/research/files/matchit.pdf
library(MatchIt)
m.out1 <- matchit(treat ~ age + educ + race + married + nodegree + re74 + re75, data = lalonde, method = "nearest", distance = "glm", caliper = 0.25, mahvars = c("age", "educ", "race", "married", "nodegree", "re74", "re75"))This code matches the specific one-to-one matching described in Lecture 4.
distance = “glm” indicates that logit propensity score is used as distance measure in matching.
caliper = 0.25 indicates that control units are drawn within a quarter of standard deviations of the distance measure.
mahvars specifies variables on which to perform Mahalanobis-metric matching within each caliper.
3.1 Balance checking
summary(m.out1) # Checking balance after matching##
## Call:
## matchit(formula = treat ~ age + educ + race + married + nodegree +
## re74 + re75, data = lalonde, method = "nearest", distance = "glm",
## mahvars = c("age", "educ", "race", "married", "nodegree",
## "re74", "re75"), caliper = 0.25)
##
## Summary of Balance for All Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance 0.5774 0.1822 1.7941 0.9211 0.3774
## age 25.8162 28.0303 -0.3094 0.4400 0.0813
## educ 10.3459 10.2354 0.0550 0.4959 0.0347
## raceblack 0.8432 0.2028 1.7615 . 0.6404
## racehispan 0.0595 0.1422 -0.3498 . 0.0827
## racewhite 0.0973 0.6550 -1.8819 . 0.5577
## married 0.1892 0.5128 -0.8263 . 0.3236
## nodegree 0.7081 0.5967 0.2450 . 0.1114
## re74 2095.5737 5619.2365 -0.7211 0.5181 0.2248
## re75 1532.0553 2466.4844 -0.2903 0.9563 0.1342
## eCDF Max
## distance 0.6444
## age 0.1577
## educ 0.1114
## raceblack 0.6404
## racehispan 0.0827
## racewhite 0.5577
## married 0.3236
## nodegree 0.1114
## re74 0.4470
## re75 0.2876
##
## Summary of Balance for Matched Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance 0.5321 0.4892 0.1948 1.2219 0.0497
## age 26.8879 25.4310 0.2036 0.5432 0.0765
## educ 10.5517 10.2931 0.1286 0.6676 0.0272
## raceblack 0.7500 0.7414 0.0237 . 0.0086
## racehispan 0.0948 0.1034 -0.0365 . 0.0086
## racewhite 0.1552 0.1552 0.0000 . 0.0000
## married 0.1983 0.2586 -0.1541 . 0.0603
## nodegree 0.7155 0.6379 0.1707 . 0.0776
## re74 2710.0046 2777.4127 -0.0138 1.4745 0.0489
## re75 1944.4850 1700.8527 0.0757 1.7308 0.0234
## eCDF Max Std. Pair Dist.
## distance 0.2672 0.2206
## age 0.2845 1.0976
## educ 0.1293 0.7375
## raceblack 0.0086 0.0237
## racehispan 0.0086 0.0365
## racewhite 0.0000 0.0000
## married 0.0603 0.4182
## nodegree 0.0776 0.3982
## re74 0.2845 0.5385
## re75 0.1034 0.5979
##
## Sample Sizes:
## Control Treated
## All 429 185
## Matched 116 116
## Unmatched 313 69
## Discarded 0 0The table of “Summary of Balance for All Data” is the balance checking results before matching.
The table of “Summary of Balance for Matched Data” is the balance checking results after matching.
The first row “distance” in each summary table indicates the chosen distance measure, which is logit propensity score here.
The third column in each summary table is the standardized mean difference. A covariate is considered to be balanced on average if the standardized mean difference is less than 0.25 and preferably less than 0.10 in magnitude.
The last table of “Sample Sizes” indicates that 116 treated units are matched with 116 control units, leaving 69 treated units and 313 control units unmatched.
3.2 Visualize the absolute standardized mean difference
plot(summary(m.out1))
If the dots are outside the range of x axis, you can adjust the range by modifying the above coding to e.g., plot(summary(m.out1), xlim = c(0, 2.5)). The first value of xlim is the lower bound, while the second is the upper bound. They can be changed based on your need.
Each circle represents the absolute standardized mean difference in the corresponding covariate before matching.
Each dot represents the absolute standardized mean difference in the corresponding covariate after matching.
The right solid vertical line indicates an absolute standardized mean difference of 0.1.
The dashed vertical line indicates an absolute standardized mean difference of 0.05.
The plot shows that all the covariates have standardized mean difference smaller than 0.25 after matching, indicating good balance is achieved after matching.
3.3 Obtain matched sample
After matching, we can save the matched sample and run a regression of the outcome on the treatment and covariates in the matched sample.
m.data <- match.data(m.out1)
head(m.data)## treat age educ race married nodegree re74 re75 re78 logit.p
## NSW1 1 37 11 black 1 1 0 0 9930.046 0.5700293
## NSW2 1 22 9 hispan 0 1 0 0 3595.894 -1.2388614
## NSW4 1 27 11 black 0 1 0 0 7506.146 1.2443718
## NSW6 1 22 9 black 0 1 0 0 4056.494 0.8428727
## NSW8 1 32 11 black 0 1 0 0 8472.158 1.3232572
## NSW9 1 22 16 black 0 0 0 0 2164.022 1.2647240
## distance weights subclass
## NSW1 0.6387699 1 1
## NSW2 0.2246342 1 2
## NSW4 0.7763241 1 3
## NSW6 0.6990699 1 4
## NSW8 0.7897231 1 5
## NSW9 0.7798382 1 6For one-to-one matching, weights (the computed matching weights) are all 1.
subclass indicates matching pair membership.
nrow(m.data) # This is for checking the sample size of the matched sample.## [1] 232Note that this matched sample only has 232 individuals.
3.4 Estimate causal effect using matching
We can then model the outcome in this matched sample using the standard regression functions in R, like lm().
Including the covariates used for matching in the regression can provide additional robustness to slight imbalances remaining after matching and can improve precision.
Including strong predictors of the outcome in the outcome model (even though they do not predict the treatment) can improve the estimation efficiency.
install.packages("marginaleffects")library(marginaleffects)
fit <- lm(re78 ~ treat * (age + educ + race + married + nodegree + re74 + re75), data = m.data, weights = weights)Finally, we use avg_comparisons() to perform g-computation to estimate the treatment effect. We recommend using cluster-robust standard errors for most analyses, with pair membership as the clustering variable. avg_comparisons() makes this straightforward.
avg_comparisons(fit,
variables = "treat",
vcov = ~subclass,
newdata = subset(m.data, treat == 1),
wts = "weights")##
## Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## 1190 1032 1.15 0.249 2.0 -833 3213
##
## Term: treat
## Type: response
## Comparison: 1 - 0Calculate the effect size by dividing the original estimate by the standard deviation of the control group’s outcome.
1190/sd(lalonde$re78[lalonde$treat == 0])## [1] 0.1631442The treatment effect is estimated to be 1190, with a standard error of 1029. The result indicates that one’s potential 1978 earnings if assigned to the National Supported Work program is $1190 higher than if assigned to the control group on average. It accounts for 16% of the standard deviation of the outcome in the control group. The confidence interval is [-826, 3206], indicating that the treatment effect is statistically insignificant.
3.5 Alternative matching method
The one-to-one matching method has 313 untreated units and 69 treated units unmatched. We may try a different matching method to keep as much as information as possible. Here, let’s try full matching. An introduction to full matching can be found in Lecture 4.
install.packages("optmatch")m.out2 <- matchit(treat ~ age + educ + race + married + nodegree + re74 + re75, data = lalonde, method = "full", distance = "glm")
summary(m.out2)##
## Call:
## matchit(formula = treat ~ age + educ + race + married + nodegree +
## re74 + re75, data = lalonde, method = "full", distance = "glm")
##
## Summary of Balance for All Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance 0.5774 0.1822 1.7941 0.9211 0.3774
## age 25.8162 28.0303 -0.3094 0.4400 0.0813
## educ 10.3459 10.2354 0.0550 0.4959 0.0347
## raceblack 0.8432 0.2028 1.7615 . 0.6404
## racehispan 0.0595 0.1422 -0.3498 . 0.0827
## racewhite 0.0973 0.6550 -1.8819 . 0.5577
## married 0.1892 0.5128 -0.8263 . 0.3236
## nodegree 0.7081 0.5967 0.2450 . 0.1114
## re74 2095.5737 5619.2365 -0.7211 0.5181 0.2248
## re75 1532.0553 2466.4844 -0.2903 0.9563 0.1342
## eCDF Max
## distance 0.6444
## age 0.1577
## educ 0.1114
## raceblack 0.6404
## racehispan 0.0827
## racewhite 0.5577
## married 0.3236
## nodegree 0.1114
## re74 0.4470
## re75 0.2876
##
## Summary of Balance for Matched Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance 0.5774 0.5762 0.0054 0.9930 0.0041
## age 25.8162 24.8095 0.1407 0.4976 0.0795
## educ 10.3459 10.3452 0.0004 0.5830 0.0206
## raceblack 0.8432 0.8347 0.0236 . 0.0086
## racehispan 0.0595 0.0657 -0.0266 . 0.0063
## racewhite 0.0973 0.0996 -0.0078 . 0.0023
## married 0.1892 0.1368 0.1338 . 0.0524
## nodegree 0.7081 0.7056 0.0056 . 0.0025
## re74 2095.5737 2363.4473 -0.0548 1.1080 0.0424
## re75 1532.0553 1632.4020 -0.0312 1.8588 0.0704
## eCDF Max Std. Pair Dist.
## distance 0.0486 0.0192
## age 0.3131 1.3111
## educ 0.0548 1.2390
## raceblack 0.0086 0.0324
## racehispan 0.0063 0.5400
## racewhite 0.0023 0.3911
## married 0.0524 0.4715
## nodegree 0.0025 0.9593
## re74 0.2492 0.8654
## re75 0.2366 0.8099
##
## Sample Sizes:
## Control Treated
## All 429. 185
## Matched (ESS) 52.11 185
## Matched 429. 185
## Unmatched 0. 0
## Discarded 0. 0plot(summary(m.out2))
The results indicate that all the observations in the sample are matched, and a better balance is achieved after matching.
4. Stratification
See page 6 of the paper: https://imai.fas.harvard.edu/research/files/matchit.pdf
m.out3 <- matchit(treat ~ age + educ + race + married + nodegree + re74 + re75, data = lalonde, method = "subclass")All the other steps stay the same as those for matching.
5. Weighting
5.1 Calculate the weight
\(\hat{w}_{i} = \frac{Pr(T=1)}{\hat{\theta}}\) for the treated
\(\hat{w}_{i} = \frac{1 - Pr(T=1)}{1-\hat{\theta}}\) for the untreated
lalonde$p = predict(l, type = "response") # save the propensity score
lalonde$ipw[lalonde$treat == 1] = mean(lalonde$treat == 1)/lalonde$p[lalonde$treat == 1] # calculate ipw for treated units
lalonde$ipw[lalonde$treat == 0] = (1 - mean(lalonde$treat == 1))/ (1 - lalonde$p[lalonde$treat == 0]) # calculate ipw for untreated units5.2 Balance checking
For each covariate \(X\), we calculate the standardized weighted mean difference as
\(\displaystyle \frac{\frac{\sum_{i=1}^{n} \hat{w}_{i}T_{i}X_i}{\sum_{i=1}^{n} \hat{w}_{i}T_{i}} - \frac{\sum_{i=1}^{n} \hat{w}_{i}(1-T_{i})X_i}{\sum_{i=1}^{n} \hat{w}_{i}(1-T_{i})}}{sd(X)}\)
Please run the following balance checking function without changing anything.
balance = function(x, data, treat){
mean.1 = sum(x * data$ipw * data[, treat])/sum(data$ipw * data[, treat])
mean.0 = sum(x * data$ipw * (1 - data[, treat]))/sum(data$ipw * (1 - data[, treat]))
mean.diff = mean.1 - mean.0
sd.mean.diff = mean.diff/sd(x)
return(sd.mean.diff)
}Check the balance of the logit of propensity score. The first element in the function is the name of the logit of propensity score, the second element is the name of the data, and the third is the name of the treatment variable.
balance(logit.p, lalonde, "treat") ## [1] 0.1738566Check the balance of each observed covariate. Below I use age, edu, and married as illustration. The first element in the function is the name of the covariate, and the other two as the same as above.
balance(lalonde$age, lalonde, "treat") ## [1] -0.1552148balance(lalonde$edu, lalonde, "treat") ## [1] 0.1217614balance(lalonde$married, lalonde, "treat") ## [1] -0.191339Good balance is achieved if the absolute value of the standardized mean difference is smaller than 0.25.
5.3 Estimate causal effect using weighting
The treatment effect is estimated as
\(\displaystyle \hat{\delta} =\frac{\sum_{i=1}^{n} \hat{w}_{i}T_{i}Y_i}{\sum_{i=1}^{n} \hat{w}_{i}T_{i}} - \frac{\sum_{i=1}^{n} \hat{w}_{i}(1-T_{i})Y_i}{\sum_{i=1}^{n} \hat{w}_{i}(1-T_{i})}\)
Because the propensity score is estimated, the standard error of the casual effect estimate contains uncertainty in the estimated propensity score.
Bootstrapping can be used to account for the uncertainty in the estimated propensity score.
When running the following estimation function, please change the regression l = glm() only.
install.packages("boot")library(boot)
ATE = function(y, treat, data, indices) { # ATE stands for average treatment effect
data.boot = data[indices,] #allows boot to select sample
l = glm(treat ~ age + educ + race + married + nodegree + re74 + re75, data = data.boot, family =
binomial(link = "logit")) # rerun the logistic regression for each bootstrapped sample
data.boot$p = predict(l, type = "response") # save the propensity score
data.boot$ipw[data.boot[, treat] == 1] = 1/ data.boot$p[data.boot[, treat] == 1]
data.boot$ipw[data.boot[, treat] == 0] = 1/ (1 - data.boot$p[data.boot[, treat] == 0])
y.1 = sum(data.boot[,y] * data.boot$ipw * data.boot[, treat])/sum(data.boot$ipw * data.boot[, treat])
y.0 = sum(data.boot[,y] * data.boot$ipw * (1 - data.boot[, treat]))/sum(data.boot$ipw * (1 - data.boot[, treat]))
return(y.1 - y.0)
}Estimate the ATE. The first element specifies the name of the outcome variable. The second specifies the name of the treatment variable. The third specifies the name of the data. The fourth is the name of the function, ATE, as defined above. The fifth element is the number of bootstrapped samples, which is usually set at 1000.
test = boot(y = "re78", treat = "treat", data = lalonde, statistic = ATE, R = 1000)
test$t0 # original point estimate## [1] 224.6763test$t0/sd(lalonde$re78[lalonde$treat == 0]) # effect size## [1] 0.03080221quantile(test$t, prob = c(0.025, 0.975)) # confidence interval## 2.5% 97.5%
## -1647.521 2176.2466. Sensitivity analysis in R
R package: causalsens
https://cran.r-project.org/web/packages/causalsens/causalsens.pdf
R package: EValue
https://cran.r-project.org/web/packages/EValue/readme/README.html
EValue works best for a binary outcome although it can be extended to other non-negative outcomes (Ding and VanderWeele, 2016).