In jbryer/psa: Applied Propensity Score Analysis with R
# options(width=70)
# options(digits=2)
# options(continue=" ")
# options(warn=-1)
knitr::opts_chunk$set(digits = 3, width = 120,
warning = FALSE, message = FALSE, error = FALSE,
fig.width = 12, fig.height=6, fig.align = 'center',
digits = 3)
# The following is to fix a DT::datatable issue with Xaringan
# https://github.com/yihui/xaringan/issues/293
options(htmltools.dir.version = FALSE, htmltools.preserve.raw = FALSE)
palette2 <- c('#fc8d62', '#66c2a5')
palette3 <- c('#fc8d62', '#66c2a5', '#8da0cb')
palette3_darker <- c('#d95f02', '#1b9e77', '#7570b3')
palette4 <- c('#1f78b4', '#33a02c', '#a6cee3', '#b2df8a')
# This style was adapted from Max Kuhn: https://github.com/rstudio-conf-2020/applied-ml
# And Rstudio::conf 2020: https://github.com/rstudio-conf-2020/slide-templates/tree/master/xaringan
# This slide deck shows a lot of the features of Xaringan: https://www.kirenz.com/slides/xaringan-demo-slides.html
# To use, add this to the slide title: `r I(hexes(c("DATA606")))`
# It will use images in the images/hex_stickers directory (i.e. the filename is the parameter)
hexes <- function(x) {
x <- rev(sort(x))
markup <- function(pkg) glue::glue('<img src="images/hex/{pkg}.png" class="title-hex">')
res <- purrr::map_chr(x, markup)
paste0(res, collapse = "")
}
printLaTeXFormula <- function(fit, digits=2) {
vars <- all.vars(fit$terms)
result <- paste0('\\hat{', vars[1], '} = ', prettyNum(fit$coefficients[[1]], digits=2))
for(i in 2:length(vars)) {
val <- fit$coefficients[[i]]
result <- paste0(result, ifelse(val < 0, ' - ', ' + '),
prettyNum(abs(val), digits=digits),
' ', names(fit$coefficients)[i])
}
return(result)
}
library(multilevelPSA)
library(Matching)
library(MatchIt)
library(multilevelPSA)
library(party)
library(PSAgraphics)
library(granovaGG)
library(rbounds)
library(rpart)
library(TriMatch)
library(psa)
library(gridExtra)
library(psych)
library(tidyverse)
library(plyr)
library(knitr)
library(rmarkdown)
theme_set(theme_bw())
data(pisana)
data(tutoring)
class: center, middle, inverse, title-slide
knitr::include_graphics('images/hex/psa.png')
r metadata$title
r metadata$subtitle
r metadata$author
Last updated: r format(Sys.Date(), '%B %d, %Y')
Agenda
- Overview of causal inference
- Counterfactuals
- Randomized control trials
- Rubin's Causal Model
- Propensity Score Analysis in Three Phases
- Phase I: Estimate propensity scores and check balance
- Phase II: Estimate causal effects
- Matching
- Stratification
- Weighting
- Phase III: Check sensitivity to unobserved confounders
- Sensitivity analysis
- Bootstrapping PSA
- Advanced topics in PSA
- Matching for Non-Binary Treatments
- Multilevel PSA
Getting Started in R
Installing the psa package from Github with dependencies = 'Enhances' should install all the packages we will use in this workshop.
install.packages('devtools')
remotes::install_github('jbryer/psa', build_vignettes = TRUE, dependencies = 'Enhances')
You can download an R script to work through the commands in this slide deck at:
https://github.com/jbryer/psa/blob/master/Slides/Intro_PSA.R
Other resources are available on Github here: https://github.com/jbryer/psa
class: inverse, middle, center
Overview of Causal Inference
Popularity of Propensity Score Analysis
data(psa_citations)
ggplot(psa_citations, aes(x = Year, y = Citations, color = Search_Term)) +
geom_path() + geom_point() +
scale_color_brewer('Search Term', palette = 'Set1') +
ylab('Number of Citations') +
theme_bw() +
ggtitle('Number of Citations for Propensity Score Analysis',
subtitle = 'Source: Web of Science and Google Scholar')
Counterfactuals
knitr::include_graphics('images/Causality.png')
The Randomized Experiment
Considered to be the gold standard for estimating causal effects.
-
Effects can be estimated using simple means between groups, or blocks in randomized block design.
-
Randomization presumes unbiasedness and balance between groups.
However, randomization is often not feasible for many reasons, especially in educational contexts.
The strong ignorability assumption states that:
$$({ Y }{ i }(1),{ Y }{ i }(0)) \; \unicode{x2AEB} \; { T }{ i }|{ X }{ i }=x$$
for all ${X}_{i}$.
class: font80
RCT Example
set.seed(2112)
pop.mean <- 100
pop.sd <- 15
pop.es <- .3
n <- 30
thedata <- data.frame(
id = 1:30,
center = rnorm(n, mean = pop.mean, sd = pop.sd),
stringsAsFactors = FALSE
)
val <- pop.sd * pop.es / 2
thedata$placebo <- thedata$center - val
thedata$treatment <- thedata$center + val
thedata$diff <- thedata$treatment - thedata$placebo
thedata$RCT_Assignment <- sample(c('placebo', 'treatment'), n, replace = TRUE)
thedata$RCT_Value <- as.numeric(apply(thedata, 1,
FUN = function(x) { return(x[x['RCT_Assignment']]) }))
head(thedata, n = 3)
tab.out <- describeBy(thedata$RCT_Value, group = thedata$RCT_Assignment, mat = TRUE, skew = FALSE)
True Counterfactual
p1 <- ggplot(thedata) +
geom_segment(aes(x = placebo, xend = treatment, y = id, yend = id)) +
geom_point(aes(x = placebo, y = id), color = 'blue') +
geom_point(aes(x = treatment, y = id), color = 'red') +
ylab('') + xlab('Outcome') +
xlim(pop.mean - 3 * pop.sd, pop.mean + 3 * pop.sd) +
ggtitle(paste0('True Counterfactual Difference = ', mean(thedata$diff)))
p1b <- p1 +
geom_vline(xintercept = mean(thedata$treatment), color = 'red') +
geom_vline(xintercept = mean(thedata$placebo), color = 'blue')
p2 <- ggplot(thedata, aes(x = RCT_Value, color = RCT_Assignment, y = id)) +
geom_point() +
scale_color_manual(values = c('placebo' = 'blue', 'treatment' = 'red')) +
theme(legend.position = 'none') +
ylab('') + xlab('Outcome') +
xlim(pop.mean - 3 * pop.sd, pop.mean + 3 * pop.sd) +
ggtitle('Observed values in an RCT')
p2b <- p2 +
geom_vline(data = tab.out, aes(xintercept = mean, color = group1)) +
ggtitle(paste0('RCT Difference = ', round(diff(tab.out$mean), digits = 2)))
.pull-left[
p1
]
True Counterfactual (left) vs. One RCT (right)
.pull-left[
p1
]
.pull-right[
p2
]
True Counterfactual (left) vs. One RCT (right)
.pull-left[
p1b
]
.pull-right[
p2b
]
Distribution of Differences from 1,000 RCTs
sim.diff <- numeric(1000)
for(i in seq_along(sim.diff)) {
treats <- sample(c(T,F), n, replace = TRUE)
sim.diff[i] <- mean(thedata[treats,]$treatment) - mean(thedata[!treats,]$placebo)
}
ggplot(data.frame(x = sim.diff), aes(x = x)) +
geom_histogram(alpha = 0.5, bins = 20) +
geom_vline(xintercept = mean(thedata$diff), color = 'red') +
geom_vline(xintercept = mean(sim.diff)) +
xlab('RCT Different') + ylab('Count')
Rubin's Causal Model
- The causal effect of a treatment is the difference in an individual's outcome under the situation they were given the treatment and not (referred to as a counterfactual).
$${\delta}{i} ={ Y }{ i1 }-{ Y }_{ i0 }$$
-
However, it is impossible to directly observe ${\delta}_{i}$ (referred to as The Fundamental Problem of Causal Inference, Holland 1986).
-
Rubin frames this problem as a "missing data problem" .font70[(see Rubin, 1974, 1977, 1978, 1980, and Holland, 1986)].
Propensity Score Analysis
The propensity score is the "conditional probability of assignment to a particular treatment given a vector of observed covariates" (Rosenbaum & Rubin, 1983, p. 41). The probability of being in the treatment:
$$\pi ({ X }{ i }) \; \equiv \; Pr({ T }{ i }=1|{ X }_{ i })$$
The balancing property under exogeneity:
$${ T }{ i } \; \unicode{x2AEB} { X }{ i } \;| \; \pi ({ X }_{ i })$$
We can then restate the ignorability assumption with the propensity score:
$$({ Y }{ i }(1),{ Y }{ i }(0)) \; \unicode{x2AEB} \; { T }{ i } \; | \; \pi({ X }{ i })$$
Treatment Effects
The average treatment effect (ATE) is defined as:
$$E({ r }{ 1 })-E({ r }{ 0 })$$
where $E(.)$ is the expectation in the population. For a set of covariates, $X$, and outcomes $Y$ where 0 denotes control and 1 treatment, we define ATE as:
$$ATE=E(Y_{1}-Y_{0}|X)=E(Y_{1}|X)-E(Y_{0}|X)$$
As we will see later there are alternative treatment effects (estimands) we can estimate instead of ATE.
What Rosenbaum and Rubin (1983) proved in their seminal paper is that the propensity score is a univariate representation of the multivariate matrix. As we will see later, two observations with very similar propensity scores will look similar across all the observed covariates.
class: inverse, middle, center
Propensity Score Analysis in Three Phases
class: font90
Simulated Example
We will simulate a dataset with three covariates, x1 and x2 which are continuous and x3 which is categorical. The assumed treatment effect is 1.5.
cols <- c('#fc8d62', '#66c2a5')
set.seed(2112)
.pull-left[
n <- 500
treatment_effect <- 1.5
X <- mvtnorm::rmvnorm(
n,
mean = c(0.5, 1, 0),
sigma = matrix(c(2, 1, 1,
1, 1, 1,
1, 1, 1),
ncol = 3) )
dat <- tibble(
x1 = X[, 1],
x2 = X[, 2],
x3 = X[, 3] > 0,
treatment = as.numeric(- 0.5 +
0.25 * x1 +
0.75 * x2 +
0.05 * x3 +
rnorm(n, 0, 1) > 0),
outcome = treatment_effect * treatment +
rnorm(n, 0, 1)
)
]
.pull-right[
head(dat, n = 6)
]
Scatterplot
ggplot(dat, aes(x = x1, y = x2, shape = x3, color = factor(treatment))) +
geom_point() + scale_color_manual('Treatment', values = cols)
Steps for Implementing Propensity Score Analysis
knitr::include_graphics('images/PSA_Flow.png')
Propensity score methods
There are three major approaches for conducting PSA:
-
Stratification Treatment and comparison units are divided into strata (or subclasses) so that treated and comparison units are similar within each strata. Cochran (1968) observed that creating five subclassifications (stratum) removes at least 90% of the bias in the estimated treatment effect.
-
Matching - Each treatment unit is paired with a comparison unit based upon the pre-treatment covariates.
-
Weighting Each observation is weighted by the inverse of the probability of being in that group.
Stratification
Stratification involves dividing (or stratifying) the observations into subgroups based upon the propensity score. Here, we used quintiles on the propensity scores where were estimated using logistic regression. For classification trees the stratum is determined by the leaf nodes.
lr.out <- glm(treatment ~ x1 + x2 + x3, data = dat, family = binomial(link='logit'))
dat$ps <- fitted(lr.out) # Get the propensity scores
# Stratification
breaks5 <- psa::get_strata_breaks(dat$ps)
dat$strata5 <- cut(x = dat$ps,
breaks = breaks5$breaks,
include.lowest = TRUE,
labels = breaks5$labels$strata)
ggplot(dat, aes(x = ps, color = as.logical(treatment))) +
geom_density(aes(fill = as.logical(treatment)), alpha = 0.2) +
geom_vline(xintercept = breaks5$breaks, alpha = 0.5) +
geom_text(data = breaks5$labels,
aes(x = xmid, y = 0, label = strata),
color = 'black', vjust = 0.8, size = 8) +
scale_fill_manual('Treatment', values = palette2) +
scale_color_manual('Treatment', values = palette2) +
xlab('Propensity Score') + ylab('Density') +
xlim(c(0, 1)) +
ggtitle('Density distribution of propensity scores by treatment',
subtitle = 'Five strata represented by vertical lines')
Stratification (cont.)
Independent sample tests (e.g. t-tests) are conducted within each stratum and pooled to provide an overall estimate.
psa::stratification_plot(ps = dat$ps,
treatment = dat$treatment,
outcome = dat$outcome)
Matching
Dependent sample tests (e.g. t-tests) are conducted using match pairs to provide a treatment.
match_out <- Matching::Match(Y = dat$outcome,
Tr = dat$treatment,
X = dat$ps,
caliper = 0.1,
estimand = 'ATE')
dat_match <- data.frame(treat_ps = dat[match_out$index.treated,]$ps,
treat_outcome = dat[match_out$index.treated,]$outcome,
control_ps = dat[match_out$index.control,]$ps,
control_outcome = dat[match_out$index.control,]$outcome)
psa::matching_plot(ps = dat$ps,
treatment = dat$treatment,
outcome = dat$outcome,
index_treated = match_out$index.treated,
index_control = match_out$index.control)
Matching Methods
There are many choices and approaches to matching, including:
- Propensity score matching.
- Limited exact matching.
- Full matching.
- Nearest neighbor matching.
- Optimal/Genetic matching.
- Mahalanobis distance matching (for quantitative covariates only).
- Matching with and without replacement.
- One-to-one or one-to-many matching.
Which method should you use?
Whichever one gives the best balance!
Weighting
Propensity score weights can be used as regression weights, the specific weights depend on the desired estimand and will be provided in later slides.
dat <- dat |> mutate(
ate_weight = psa::calculate_ps_weights(treatment, ps, estimand = 'ATE'),
att_weight = psa::calculate_ps_weights(treatment, ps, estimand = 'ATT'),
atc_weight = psa::calculate_ps_weights(treatment, ps, estimand = 'ATC'),
atm_weight = psa::calculate_ps_weights(treatment, ps, estimand = 'ATM')
)
psa::weighting_plot(ps = dat$ps,
treatment = dat$treatment,
outcome = dat$outcome)
Shiny Application
We can explore how these three plots change as the treatment effects change using the psa::psa_simulation_shiny() application.
psa::psa_simulation_shiny()
knitr::include_graphics('images/psa_simulation_screenshot.png')
Phase I: Estimate Propensity Scores
.pull-left[
In this example we will use logistic regression to estimate the propensity scores.
lr.out <- glm(
treatment ~ x1 + x2 + x3,
data = dat,
family = binomial(link='logit'))
dat$ps <- fitted(lr.out) # Propensity scores
For stratification we will use quintiles to split the observations into five equal groups.
breaks5 <- psa::get_strata_breaks(dat$ps)
dat$strata5 <- cut(
x = dat$ps,
breaks = breaks5$breaks,
include.lowest = TRUE,
labels = breaks5$labels$strata)
]
.pull-right[.font70[
summary(lr.out)
]]
Distribution of Propensity Scores
ggplot(dat) +
geom_histogram(data = dat[dat$treatment == 1,], aes(x = ps, y = after_stat(count)), bins = 50, fill = cols[2]) +
geom_histogram(data = dat[dat$treatment == 0,], aes(x = ps, y = -after_stat(count)), bins = 50, fill = cols[1]) +
geom_hline(yintercept = 0, lwd = 0.5) + scale_y_continuous(label = abs)
Check Balance: Multiple Covariates
PSAgraphics::cv.bal.psa(dat[,1:3], dat$treatment, dat$ps, strata = 5)
Check Balance: Single Covariate
.pull-left[
PSAgraphics::box.psa(dat$x1,
dat$treatment,
dat$strata5)
]
.pull-right[
PSAgraphics::cat.psa(dat$x3,
dat$treatment,
dat$strata5)
]
PS Weights for Understanding Treatment Effects
Given that the distribution of treatment and control observations across the propensity score range are not the same, there are a number of alternative estimates of treatment effect. We will explore three additional esimates in addition to the classic average treatment effect.
dat <- dat |> mutate(
ate_weight = psa::calculate_ps_weights(treatment, ps, estimand = 'ATE'),
att_weight = psa::calculate_ps_weights(treatment, ps, estimand = 'ATT'),
atc_weight = psa::calculate_ps_weights(treatment, ps, estimand = 'ATC'),
atm_weight = psa::calculate_ps_weights(treatment, ps, estimand = 'ATM')
)
dat |> head(n = 4)
Average Treatment Effect (ATE)
$$ATE = E(Y_1 - Y_0 | X) = E(Y_1|X) - E(Y_0|X)$$
ggplot() +
geom_histogram(data = dat[dat$treatment == 1,],
aes(x = ps, y = after_stat(count)),
bins = 50, alpha = 0.5) +
geom_histogram(data = dat[dat$treatment == 1,],
aes(x = ps, weight = ate_weight, y = after_stat(count)),
bins = 50,
fill = cols[2], alpha = 0.5) +
geom_histogram(data = dat[dat$treatment == 0,],
aes(x = ps, y = -after_stat(count)),
bins = 50, alpha = 0.5) +
geom_histogram(data = dat[dat$treatment == 0,],
aes(x = ps, weight = ate_weight, y = -after_stat(count)),
bins = 50,
fill = cols[1], alpha = 0.5) +
ggtitle('Average Treatment Effect (ATE)')
Average Treatment Effect Among the Treated (ATT)
$$ATT=E(Y_{1}-Y_{0}|X,C=1)=E(Y_{1}|X,C=1)-E(Y_{0}|X,C=1)$$
ggplot() +
geom_histogram(data = dat[dat$treatment == 1,],
aes(x = ps, y = after_stat(count)),
bins = 50, alpha = 0.5) +
geom_histogram(data = dat[dat$treatment == 1,],
aes(x = ps, weight = att_weight, y = after_stat(count)),
bins = 50,
fill = cols[2], alpha = 0.5) +
geom_histogram(data = dat[dat$treatment == 0,],
aes(x = ps, y = -after_stat(count)),
bins = 50, alpha = 0.5) +
geom_histogram(data = dat[dat$treatment == 0,],
aes(x = ps, weight = att_weight, y = -after_stat(count)),
bins = 50,
fill = cols[1], alpha = 0.5) +
ggtitle('Average Treatment Effect Among the Treated (ATT)')
Average Treatment Effect Among the Control (ATC)
$$ATC = E(Y_1 - Y_0 | X = 0) = E(Y_1 | X = 0) - E(Y_0 | X = 0)$$
ggplot() +
geom_histogram(data = dat[dat$treatment == 1,],
aes(x = ps, y = after_stat(count)),
bins = 50, alpha = 0.5) +
geom_histogram(data = dat[dat$treatment == 1,],
aes(x = ps, weight = atc_weight, y = after_stat(count)),
bins = 50,
fill = cols[2], alpha = 0.5) +
geom_histogram(data = dat[dat$treatment == 0,],
aes(x = ps, y = -after_stat(count)),
bins = 50, alpha = 0.5) +
geom_histogram(data = dat[dat$treatment == 0,],
aes(x = ps, weight = atc_weight, y = -after_stat(count)),
bins = 50,
fill = cols[1], alpha = 0.5) +
ggtitle('Average Treatment Effect Among the Control (ATC)')
Average Treatment Effect Among the Evenly Matched
$$ATM_d = E(Y_1 - Y_0 | M_d = 1)$$
ggplot() +
geom_histogram(data = dat[dat$treatment == 1,],
aes(x = ps, y = after_stat(count)),
bins = 50, alpha = 0.5) +
geom_histogram(data = dat[dat$treatment == 1,],
aes(x = ps, weight = atm_weight, y = after_stat(count)),
bins = 50,
fill = cols[2], alpha = 0.5) +
geom_histogram(data = dat[dat$treatment == 0,],
aes(x = ps, y = -after_stat(count)),
bins = 50, alpha = 0.5) +
geom_histogram(data = dat[dat$treatment == 0,],
aes(x = ps, weight = atm_weight, y = -after_stat(count)),
bins = 50,
fill = cols[1], alpha = 0.5) +
ggtitle('Average Treatment Effect Among the Evenly Matched (ACM)')
Treatment Effects for Weighting
$$Treatment\ Effect = \frac{\sum Y_{i}Z_{i}w_{i}}{\sum Z_{i} w_{i}} - \frac{\sum Y_{i}(1 - Z_{i}) w_{i}}{\sum (1 - Z_{i}) w_{i} }$$
Where $w$ is the weight (as defined in the following sections), $Z_i$ is the treatment assignment such that $Z = 1$ is treatment and $Z = 0$ is control, and $Y_i$ is the outcome
.pull-left[
$$w_{ATE} = \frac{Z_i}{\pi_i} + \frac{1 - Z_i}{1 - \pi_i}$$
$$w_{ATT} = \frac{\pi_i Z_i}{\pi_i} + \frac{\pi_i (1 - Z_i)}{1 - \pi_i}$$
]
.pull-right[
$$w_{ATC} = \frac{(1 - \pi_i) Z_i}{\pi_i} + \frac{(1 - e_i)(1 - Z_i)}{1 - \pi_i}$$
$$w_{ATM} = \frac{min{\pi_i, 1 - \pi_i}}{Z_i \pi_i (1 - Z_i)(1 - \pi_i)}$$
]
class: font90
Treatment Effects
.pull-left[
.font80[Average Treatment Effect]
psa::treatment_effect(
treatment = dat$treatment,
outcome = dat$outcome,
weights = dat$ate_weight)
lm(outcome ~ treatment,
data = dat,
weights = dat$ate_weight)
]
.pull-right[
.font80[Average Treatment Effect Among the Treated]
psa::treatment_effect(
treatment = dat$treatment,
outcome = dat$outcome,
weights = dat$att_weight)
lm(outcome ~ treatment,
data = dat,
weights = dat$att_weight)
]
class: font90
Treatment Effects (cont.)
.pull-left[
.font80[Average Treatment Effect Among the Control]
psa::treatment_effect(
treatment = dat$treatment,
outcome = dat$outcome,
weights = dat$atc_weight)
lm(outcome ~ treatment,
data = dat,
weights = dat$atc_weight)
]
.pull-right[
.font80[Average Treatment Effect Among the Evenly Matched]
psa::treatment_effect(
treatment = dat$treatment,
outcome = dat$outcome,
weights = dat$atm_weight)
lm(outcome ~ treatment,
data = dat,
weights = dat$atm_weight)
]
class: inverse, middle, center
Example: National Supported Work Demonstration
National Supported Work
The National Supported Work (NSW) Demonstration was a federally and privately funded randomized experiment done in the 1970s to estimate the effects of a job training program for disadvantaged workers.
- Participants were randomly selected to participate in the training program.
- Both groups were followed up to determine the effect of the training on wages.
- Analysis of the mean differences (unbiased given randomization), was approximately $800.
Lalonde (1986) used data from the Panel Survey of Income Dynamics (PSID) and the Current Population Survey (CPS) to investigate whether non-experimental methods would result in similar results to the randomized experiment. He found results ranging from $700 to $16,000.
National Supported Work (cont.)
Dehejia and Wahba (1999) later used propensity score matching to analyze the data. The found that,
- Comparison groups selected by Lalonde were very dissimilar to the treated group.
- By restricting the comparison group to those that were similar to the treated group, they could replicate the original NSW results.
- Using the CPS data, the range of treatment effect was between $1,559 to $1,681. The experimental results for the sample sample was approximately $1,800.
The covariates available include: age, education level, high school degree, marital status, race, ethnicity, and earning sin 1974 and 1975.
Outcome of interest is earnings in 1978.
data(lalonde, package='Matching')
class: font90
Estimating Propensity Scores
.pull-left[
Estimate propensity scores using logistic regression.
lalonde.formu <- treat ~ age + educ + black + hisp +
married + nodegr + re74 + re75
glm1 <- glm(lalonde.formu,
data = lalonde,
family = binomial(link = 'logit'))
Get the propensity scores:
lalonde$ps <- fitted(glm1)
Define the stratification:
strata5 <- cut(lalonde$ps,
quantile(lalonde$ps, seq(0, 1, 1/5)),
include.lowest = TRUE,
labels = letters[1:5])
]
.pull-right[
.font70[
summary(glm1)
]
]
Checking Balance: Covariate Balance Plot
covars <- all.vars(lalonde.formu)
covars <- lalonde[,covars[2:length(covars)]]
cv.bal.psa(covars, lalonde$treat, lalonde$ps, strata = 5)
Checking Balance: Continuous Covariates
.pull-left[
box.psa(lalonde$age, lalonde$treat, strata5)
]
.pull-right[
box.psa(lalonde$re74, lalonde$treat, strata5)
]
Checking Balance: Continuous Covariates (cont.)
.pull-left[
box.psa(lalonde$educ, lalonde$treat, strata5)
]
.pull-right[
box.psa(lalonde$re75, lalonde$treat, strata5)
]
Checking Balance: Categorical Covariates
.pull-left[
cat.psa(lalonde$married, lalonde$treat, strata5)
]
.pull-right[
cat.psa(lalonde$hisp, lalonde$treat, strata5)
]
Checking Balance: Categorical Covariates (cont.)
.pull-left[
cat.psa(lalonde$black, lalonde$treat, strata5)
]
.pull-right[
cat.psa(lalonde$nodegr, lalonde$treat, strata5)
]
Loess Regression
psadf <- data.frame(ps = lalonde$ps, Y = lalonde$re78, Tr = lalonde$treat)
loess_plot(ps = psadf[psadf$Y < 30000,]$ps,
outcome = psadf[psadf$Y < 30000,]$Y,
treatment = as.logical(psadf[psadf$Y < 30000,]$Tr))
Stratification
.pull-left[
psa::stratification_plot(ps = psadf$ps,
treatment = psadf$Tr,
outcome = psadf$Y,
n_strata = 5)
]
.pull-right[
psa::stratification_plot(ps = psadf$ps,
treatment = psadf$Tr,
outcome = psadf$Y,
n_strata = 10)
]
Stratification (cont.)
.pull-left[
strata5 <- cut(lalonde$ps,
quantile(lalonde$ps, seq(0, 1, 1/5)),
include.lowest = TRUE,
labels = letters[1:5])
circ.psa(lalonde$re78, lalonde$treat, strata5)
]
.pull-right[
strata10 <- cut(lalonde$ps,
quantile(lalonde$ps, seq(0, 1, 1/10)),
include.lowest = TRUE,
labels = letters[1:10])
circ.psa(lalonde$re78, lalonde$treat, strata10)
]
class: font60
Stratification (cont.)
.pull-left[
circ.psa(lalonde$re78, lalonde$treat, strata5)
]
.pull-right[
circ.psa(lalonde$re78, lalonde$treat, strata10)
]
Matching
rr <- Match(Y = lalonde$re78,
Tr = lalonde$treat,
X = lalonde$ps,
M = 1,
estimand = 'ATT',
ties = FALSE)
summary(rr)
Visualizing Matching Results
matches <- data.frame(Treat = lalonde[rr$index.treated,'re78'],
Control = lalonde[rr$index.control,'re78'])
granovagg.ds(matches[,c('Control','Treat')], xlab = 'Treat', ylab = 'Control')
Balance for Matching
psa::MatchBalance(df = lalonde, formu = lalonde.formu,
formu.Y = update.formula(lalonde.formu, re78 ~ .),
M = 1, estimand = 'ATT', ties = FALSE) |> plot()
Balance for Matching (cont.)
psa::MatchBalance(df = lalonde, formu = lalonde.formu,
formu.Y = update.formula(lalonde.formu, re78 ~ .),
exact.covs = c('nodegr'), #<<
M = 1, estimand = 'ATT', ties = FALSE) |> plot()
class: inverse, middle, center
Sensitivity Analysis
Sensitivity Analysis
- An observational study is free of hidden bias if the propensity scores for each subject depend only on the observed covariates.
- That is, the p-value is valid if there are no unobserved confounders.
- However, there are very likely covariates that would better model treatment. These introduce hidden bias.
- Hidden bias exists if two subjects have the same covariates, but different propensity scores.
$X_a = X_b$ but ${ \pi }{ a }\neq { \pi }{ b }$ for some a and b.
Sensitivity Analysis
Each person in the treatment is matched to exactly one person in the control. The odds of being in the treatment for persons a and b are:
$O_a = \frac{ \pi_a }{ 1 - \pi_a }$ and $O_b = \frac{ \pi_b }{ 1 - \pi_b }$
The ratio of these odds, $\Gamma$, measures the bias after matching.
$$\Gamma =\frac { { O }{ a } }{ { O }{ b } } =\frac { { { \pi }{ a } / ( }{ 1-{ \pi }{ a }) } }{ { { \pi }{ b } / (1-{ \pi }{ b }) } }$$
This is the ratio of the odds the treated unit being in the treatment group to the matched control unit being in the treatment group.
Sensitivity Analysis
Sensitivity analysis tests whether the results hold for various ranges of $\Gamma$. That is, we test how large the differences in $\pi$ (i.e. propensity scores) would have to be to change our basic inference. Let $p_a$ and $p_b$ be the probability of each unit of the matched pair being treated, conditional on exactly one being treated. For example:
- If $\Gamma = 1$, the treatment and control unit within each pair has the same value of treatment assignment ( $p_a = 0.5$ and $p_b = 0.5$).
- If $\frac{1}{2} \le \Gamma \le 2$, no unit can be more than twice as likely as its match to get treated ( $0.33 \le p_a$, $p_b \le 0.66$).
- If $\frac{1}{3} \le \Gamma \le 3$, no unit can be more than three times as likely as its match to get treated ( $0.25 \le p_a$, $p_b \le 0.75$)
To get the bounds:
$$ \frac{1}{\Gamma +1 } \le p_a, p_b \le \frac{\Gamma}{\Gamma +1} $$
Wilcoxon Signed Rank Test
-
Drop pairs where the matches have the same outcome.
-
Calculate the difference in outcomes within each pair.
-
Rank the pairs from smallest absolute difference to largest absolute difference (i.e. the smallest = 1).
-
Take the sum of the ranks where the treated unit had the higher outcome.
$$W=\left| \sum { 1 }^{ { N }{ r } }{ sgn({ x }{ T,i }-{ x }{ C,i })\cdot { R }{ i } } \right|$$
Where $N$ is the number of ranked pairs; $R_i$ is the rank for pair r; $x{T,i}$ and $x_{C,i}$ are the outcomes for the $i^{th}$ treated and control pair, respectively.
Sensitivity Analysis
The process for sensitivity analysis:
-
Select a series of values for $\Gamma$. For social science research, values between 1 and 2 is an appropriate start.
-
For each $\Gamma$, estimate the p-values to see how the p-values increase for larger values of $\Gamma$.
-
For binary outcomes, use McNemar's test, for all others use Wilcoxon sign rank test and the Hodges-Lehmann point estimate. See Keele (2010) for more information.
Sensitivity Analysis
Children of parents who had worked in a factory where lead was used in making batteries were matched by age, exposure to traffic, and neighborhood with children whose parents did not work in lead-related industries. Whole blood was assessed for lead content yielding measurements in mg/dl
require(rbounds)
psens(lalonde$re78[rr$index.treated],
lalonde$re78[rr$index.control],
Gamma = 2, GammaInc = 0.1)
class: inverse, middle, center
Bootstrapping Propensity Score Analysis
Boostrapping
-
Bootstrapping was first introduced by Efron (1979) in Bootstrap Methods: Another Look at the Jackknife.
-
Estimates confidence of statistics by resampling with replacement.
-
The bootstrap sample provides an estimate of the sampling distribution.
For PSA, sensitivity analysis is only well defined for matched samples.
Rosenbaum (2012) suggested that one way to test for sensitivity of model selection is to test the null hypothesis twice.
Bootstrapping Propensity Score Analysis
The PSAboot implements bootstrapping for propensity score analysis.
-
A stratified bootstrap sample is drawn to ensure the ratio of treatment-to-control observations is the same (i.e. sampling with replacement is done for the treatment and control observations is done separately). Note that the control.ratio and treat.ratio parameters allow for under sampling in the case of imbalanced data.
-
For each bootstrap sample balance statistics and treatment effects are estimated using each method (five by default).
-
Overall treatment effect with confidence interval is estimated from the bootstrap samples.
PSA Bootstrapping Example
library(PSAboot)
psaboot <- PSAboot(Tr = lalonde$treat,
Y = lalonde$re78,
X = lalonde,
formu = lalonde.formu)
summary(psaboot)
Checking Balance
psaboot_bal <- balance(psaboot)
plot(psaboot_bal)
Plotting all Boostrap Samples
plot(psaboot)
Boxplot of Boostrapping PSA
boxplot(psaboot)
Comparing Across PSA Methods
matrixplot(psaboot)
class: inverse, middle, center
Matching of Non-Binary Treatments
Matching of Non-Binary Treatments
require(TriMatch)
data(tutoring)
formu <- ~ Gender + Ethnicity + Military + ESL + EdMother + EdFather + Age +
Employment + Income + Transfer + GPA
tutoring.tpsa <- trips(tutoring, tutoring$treat, formu)
tutoring.matched.n <- trimatch(tutoring.tpsa, method=OneToN, M1=5, M2=3)
- The
TriMatch package provides functions for finding matched triplets.
- Estimates propensity scores for three separate logistic regression models (one for each pair of groups, that is, treat1-to-control, treat2-to-control, and treat1-to-treat2).
- Finds matched triplets that minimize the total distance (i.e. sum of the standardized distance between propensity scores within the three models). within a caliper.
- Provides multiple methods for determining which matched triplets are retained:
- Optimal which attempts to retain all treatment units.
- Full which retains all matched triplets within the specified caliper (.25 by default as suggested by Rosenbaum).
- Analog of the one-to-many for matched triplets. Specify how many times each treat1 and treat2 unit can be matched.
- Unique which allows each unit to be matched once, and only once.
- Functions for conducting repeated measures ANOVA and Freidman Ranksum Tests are provided.
Example: Tutoring
Students can opt to utilize tutoring services to supplement math courses. Of those who used tutoring services, approximately 58% of students used the tutoring service once, whereas the remaining 42% used it more than once. Outcome of interest is course grade.
- Military Active military status.
- Income Income level.
- Employment Employment level.
- NativeEnglish Is English their native language
- EdLevelMother Education level of their mother.
- EdLevelFather Education level of their father.
- Ethnicity American Indian or Alaska Native, Asian, Black or African American, Hispanic, Native Hawaiian or Other Pacific Islander, Two or more races, Unknown, White
- Gender Male, Female
- Age Age at course start.
- GPA Student GPA at the beginning of the course.
New Student Outreach: Covariates
Newly enrolled students received outreach contacts until they registered for a course or six months have passed, whichever came first. Outreach was conducted by two academic advisors and a comparison group was drawn from students who enrolled prior to the start of the outreach program. Outcome of interest is number of credits attempted within the first seven months of enrollment.
PSA for Non-Binary Treatments
The TriMatch algorithm works as follows:
- Estimate three separate propensity score models for each pair of groups (i.e. Control-to-Treat1, Control-to-Treat2, Treat1-to-Treat2).
- Determine the matching order. The default is to start with the largest of two treatments, then the other treatment, followed by the control.
- For each unit in group 1, find all units from group 2 within a certain threshold (i.e. difference between PSs is within a specified caliper).
- For each unit in group 2, find all units from group 3 within a certain threshold.
- Calculate the distance (difference) between each unit 3 found and the original unit 1. Eliminate candidates that exceed the caliper.
- Calculate a total distance (sum of the three distances) and retain the smallest unique M group 1 units (by default M=2)
Matching Triplets
plot(tutoring.matched.n, rows=c(50), draw.segments=TRUE)
Checking Balance
multibalance.plot(tutoring.tpsa, grid=TRUE)
Results
boxdiff.plot(tutoring.matched.n, tutoring$Grade)
class: inverse, middle, center
Multilevel Propensity Score Analysis
Multilevel PSA
The use of PSA for clustered, or multilevel data, has been limited (Thoemmes \& Felix, 2011). Bryer and Pruzek (2012, 2013) have introduced an approach to analyzing multilevel or clustered data using stratification methods and implemented in the multilevelPSA R package.
- Exact and partially exact matching methods implicitly adjust for clustering. That is, the covariates chosen to exactly match are, in essence, clustering variables.
- Exact matching only applies to phase I of PSA. How are the clusters related to outcome of interest.
The multilevelPSA uses stratification methods (e.g. quintiles, classification trees) by:
- Estimate separate propensity scores for each cluster.
- Identify strata within each cluster (e.g. leaves of classification trees, quintiles).
- Estimate ATE (or ATT) within each cluster.
- Aggregate estimated ATE to provide an overall ATE estimate.
- Several functions to summarize and visualize results and check balance.
The Programme of International Student Assessment
-
International assessment conducted by the Organization for Economic Co-operation and Development (OECD).
-
Assesses students towards the end of secondary school (approximately 15-year-old children) in math, reading, and science.
-
Collects a robust set of background information from students, parents, teachers, and schools.
-
Assess both private and public school students in many countries.
-
We will use PISA to estimate the effects of private school attendance on PISA outcomes.
Phase I of Multilevel PSA
The multilevelPSA provides two functions, mlpsa.ctree and mlpsa.logistic, that will estimate propensity scores using classification trees and logistic regression, respectively. Since logistic regression requires a complete dataset (i.e. no missing values), we will use classification trees in this example.
data(pisana)
data(pisa.colnames)
data(pisa.psa.cols)
student = pisana
mlctree = mlpsa.ctree(student[,c('CNT','PUBPRIV',pisa.psa.cols)],
formula=PUBPRIV ~ ., level2='CNT')
student.party = getStrata(mlctree, student, level2='CNT')
student.party$mathscore = apply(
student.party[,paste0('PV', 1:5, 'MATH')], 1, sum) / 5
To assess what covariates were used in each tree model, as well as the relative importance, we can create a heat map of covariate usage by level.
Covariate Heat Map
tree.plot(mlctree, level2Col=student$CNT, colLabels=pisa.colnames[,c('Variable','ShortDesc')])
Phase II of Multilevel PSA
The mlpsa function will compare the outcome of interest.
results.psa.math = mlpsa(response=student.party$mathscore,
treatment=student.party$PUBPRIV, strata=student.party$strata,
level2=student.party$CNT, minN=5)
results.psa.math$overall.wtd
results.psa.math$overall.ci
results.psa.math$level2.summary[,c('level2','Private','Private.n',
'Public','Public.n','diffwtd','ci.min','ci.max')]
Multilevel PSA Assessment Plot
The multilevel PSA assessment plot is an extension of the circ.psa plot in PSAgraphics introduced by Helmreich and Pruzek (2009).
plot(results.psa.math)
Multilevel PSA Difference Plot
mlpsa.difference.plot(results.psa.math, sd=mean(student.party$mathscore, na.rm=TRUE))
class: inverse, middle, center
Shiny Application
Shiny Application
psa::psa_shiny()
knitr::include_graphics('images/shiny_screenshot.png', dpi=350)
class: inverse, right, middle, hide-logo
Thank You!
r icons::fontawesome("paper-plane") jason.bryer@cuny.edu
r icons::fontawesome("github") @jbryer
r icons::fontawesome('mastodon') @jbryer@vis.social
r icons::fontawesome("link") psa.bryer.org
r icons::fontawesome("link") github.com/jbryer/psa
jbryer/psa documentation built on Nov. 17, 2023, 8:21 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# options(width=70) # options(digits=2) # options(continue=" ") # options(warn=-1) knitr::opts_chunk$set(digits = 3, width = 120, warning = FALSE, message = FALSE, error = FALSE, fig.width = 12, fig.height=6, fig.align = 'center', digits = 3) # The following is to fix a DT::datatable issue with Xaringan # https://github.com/yihui/xaringan/issues/293 options(htmltools.dir.version = FALSE, htmltools.preserve.raw = FALSE) palette2 <- c('#fc8d62', '#66c2a5') palette3 <- c('#fc8d62', '#66c2a5', '#8da0cb') palette3_darker <- c('#d95f02', '#1b9e77', '#7570b3') palette4 <- c('#1f78b4', '#33a02c', '#a6cee3', '#b2df8a') # This style was adapted from Max Kuhn: https://github.com/rstudio-conf-2020/applied-ml # And Rstudio::conf 2020: https://github.com/rstudio-conf-2020/slide-templates/tree/master/xaringan # This slide deck shows a lot of the features of Xaringan: https://www.kirenz.com/slides/xaringan-demo-slides.html # To use, add this to the slide title: `r I(hexes(c("DATA606")))` # It will use images in the images/hex_stickers directory (i.e. the filename is the parameter) hexes <- function(x) { x <- rev(sort(x)) markup <- function(pkg) glue::glue('<img src="images/hex/{pkg}.png" class="title-hex">') res <- purrr::map_chr(x, markup) paste0(res, collapse = "") } printLaTeXFormula <- function(fit, digits=2) { vars <- all.vars(fit$terms) result <- paste0('\\hat{', vars[1], '} = ', prettyNum(fit$coefficients[[1]], digits=2)) for(i in 2:length(vars)) { val <- fit$coefficients[[i]] result <- paste0(result, ifelse(val < 0, ' - ', ' + '), prettyNum(abs(val), digits=digits), ' ', names(fit$coefficients)[i]) } return(result) } library(multilevelPSA) library(Matching) library(MatchIt) library(multilevelPSA) library(party) library(PSAgraphics) library(granovaGG) library(rbounds) library(rpart) library(TriMatch) library(psa) library(gridExtra) library(psych) library(tidyverse) library(plyr) library(knitr) library(rmarkdown) theme_set(theme_bw()) data(pisana) data(tutoring)
class: center, middle, inverse, title-slide
knitr::include_graphics('images/hex/psa.png')
r metadata$title
r metadata$subtitle
r metadata$author
Last updated: r format(Sys.Date(), '%B %d, %Y')
Agenda
- Overview of causal inference
- Counterfactuals
- Randomized control trials
- Rubin's Causal Model
- Propensity Score Analysis in Three Phases
- Phase I: Estimate propensity scores and check balance
- Phase II: Estimate causal effects
- Matching
- Stratification
- Weighting
- Phase III: Check sensitivity to unobserved confounders
- Sensitivity analysis
- Bootstrapping PSA
- Advanced topics in PSA
- Matching for Non-Binary Treatments
- Multilevel PSA
Getting Started in R
Installing the psa package from Github with dependencies = 'Enhances' should install all the packages we will use in this workshop.
install.packages('devtools') remotes::install_github('jbryer/psa', build_vignettes = TRUE, dependencies = 'Enhances')
You can download an R script to work through the commands in this slide deck at:
https://github.com/jbryer/psa/blob/master/Slides/Intro_PSA.R
Other resources are available on Github here: https://github.com/jbryer/psa
class: inverse, middle, center
Overview of Causal Inference
Popularity of Propensity Score Analysis
data(psa_citations) ggplot(psa_citations, aes(x = Year, y = Citations, color = Search_Term)) + geom_path() + geom_point() + scale_color_brewer('Search Term', palette = 'Set1') + ylab('Number of Citations') + theme_bw() + ggtitle('Number of Citations for Propensity Score Analysis', subtitle = 'Source: Web of Science and Google Scholar')
Counterfactuals
knitr::include_graphics('images/Causality.png')
The Randomized Experiment
Considered to be the gold standard for estimating causal effects.
-
Effects can be estimated using simple means between groups, or blocks in randomized block design.
-
Randomization presumes unbiasedness and balance between groups.
However, randomization is often not feasible for many reasons, especially in educational contexts.
The strong ignorability assumption states that:
$$({ Y }{ i }(1),{ Y }{ i }(0)) \; \unicode{x2AEB} \; { T }{ i }|{ X }{ i }=x$$
for all ${X}_{i}$.
class: font80
RCT Example
set.seed(2112) pop.mean <- 100 pop.sd <- 15 pop.es <- .3 n <- 30 thedata <- data.frame( id = 1:30, center = rnorm(n, mean = pop.mean, sd = pop.sd), stringsAsFactors = FALSE ) val <- pop.sd * pop.es / 2 thedata$placebo <- thedata$center - val thedata$treatment <- thedata$center + val thedata$diff <- thedata$treatment - thedata$placebo thedata$RCT_Assignment <- sample(c('placebo', 'treatment'), n, replace = TRUE) thedata$RCT_Value <- as.numeric(apply(thedata, 1, FUN = function(x) { return(x[x['RCT_Assignment']]) })) head(thedata, n = 3) tab.out <- describeBy(thedata$RCT_Value, group = thedata$RCT_Assignment, mat = TRUE, skew = FALSE)
True Counterfactual
p1 <- ggplot(thedata) + geom_segment(aes(x = placebo, xend = treatment, y = id, yend = id)) + geom_point(aes(x = placebo, y = id), color = 'blue') + geom_point(aes(x = treatment, y = id), color = 'red') + ylab('') + xlab('Outcome') + xlim(pop.mean - 3 * pop.sd, pop.mean + 3 * pop.sd) + ggtitle(paste0('True Counterfactual Difference = ', mean(thedata$diff))) p1b <- p1 + geom_vline(xintercept = mean(thedata$treatment), color = 'red') + geom_vline(xintercept = mean(thedata$placebo), color = 'blue') p2 <- ggplot(thedata, aes(x = RCT_Value, color = RCT_Assignment, y = id)) + geom_point() + scale_color_manual(values = c('placebo' = 'blue', 'treatment' = 'red')) + theme(legend.position = 'none') + ylab('') + xlab('Outcome') + xlim(pop.mean - 3 * pop.sd, pop.mean + 3 * pop.sd) + ggtitle('Observed values in an RCT') p2b <- p2 + geom_vline(data = tab.out, aes(xintercept = mean, color = group1)) + ggtitle(paste0('RCT Difference = ', round(diff(tab.out$mean), digits = 2)))
.pull-left[
p1
]
True Counterfactual (left) vs. One RCT (right)
.pull-left[
p1
] .pull-right[
p2
]
True Counterfactual (left) vs. One RCT (right)
.pull-left[
p1b
] .pull-right[
p2b
]
Distribution of Differences from 1,000 RCTs
sim.diff <- numeric(1000) for(i in seq_along(sim.diff)) { treats <- sample(c(T,F), n, replace = TRUE) sim.diff[i] <- mean(thedata[treats,]$treatment) - mean(thedata[!treats,]$placebo) } ggplot(data.frame(x = sim.diff), aes(x = x)) + geom_histogram(alpha = 0.5, bins = 20) + geom_vline(xintercept = mean(thedata$diff), color = 'red') + geom_vline(xintercept = mean(sim.diff)) + xlab('RCT Different') + ylab('Count')
Rubin's Causal Model
- The causal effect of a treatment is the difference in an individual's outcome under the situation they were given the treatment and not (referred to as a counterfactual).
$${\delta}{i} ={ Y }{ i1 }-{ Y }_{ i0 }$$
-
However, it is impossible to directly observe ${\delta}_{i}$ (referred to as The Fundamental Problem of Causal Inference, Holland 1986).
-
Rubin frames this problem as a "missing data problem" .font70[(see Rubin, 1974, 1977, 1978, 1980, and Holland, 1986)].
Propensity Score Analysis
The propensity score is the "conditional probability of assignment to a particular treatment given a vector of observed covariates" (Rosenbaum & Rubin, 1983, p. 41). The probability of being in the treatment: $$\pi ({ X }{ i }) \; \equiv \; Pr({ T }{ i }=1|{ X }_{ i })$$
The balancing property under exogeneity:
$${ T }{ i } \; \unicode{x2AEB} { X }{ i } \;| \; \pi ({ X }_{ i })$$
We can then restate the ignorability assumption with the propensity score:
$$({ Y }{ i }(1),{ Y }{ i }(0)) \; \unicode{x2AEB} \; { T }{ i } \; | \; \pi({ X }{ i })$$
Treatment Effects
The average treatment effect (ATE) is defined as:
$$E({ r }{ 1 })-E({ r }{ 0 })$$
where $E(.)$ is the expectation in the population. For a set of covariates, $X$, and outcomes $Y$ where 0 denotes control and 1 treatment, we define ATE as:
$$ATE=E(Y_{1}-Y_{0}|X)=E(Y_{1}|X)-E(Y_{0}|X)$$
As we will see later there are alternative treatment effects (estimands) we can estimate instead of ATE.
What Rosenbaum and Rubin (1983) proved in their seminal paper is that the propensity score is a univariate representation of the multivariate matrix. As we will see later, two observations with very similar propensity scores will look similar across all the observed covariates.
class: inverse, middle, center
Propensity Score Analysis in Three Phases
class: font90
Simulated Example
We will simulate a dataset with three covariates, x1 and x2 which are continuous and x3 which is categorical. The assumed treatment effect is 1.5.
cols <- c('#fc8d62', '#66c2a5') set.seed(2112)
.pull-left[
n <- 500 treatment_effect <- 1.5 X <- mvtnorm::rmvnorm( n, mean = c(0.5, 1, 0), sigma = matrix(c(2, 1, 1, 1, 1, 1, 1, 1, 1), ncol = 3) ) dat <- tibble( x1 = X[, 1], x2 = X[, 2], x3 = X[, 3] > 0, treatment = as.numeric(- 0.5 + 0.25 * x1 + 0.75 * x2 + 0.05 * x3 + rnorm(n, 0, 1) > 0), outcome = treatment_effect * treatment + rnorm(n, 0, 1) )
] .pull-right[
head(dat, n = 6)
]
Scatterplot
ggplot(dat, aes(x = x1, y = x2, shape = x3, color = factor(treatment))) + geom_point() + scale_color_manual('Treatment', values = cols)
Steps for Implementing Propensity Score Analysis
knitr::include_graphics('images/PSA_Flow.png')
Propensity score methods
There are three major approaches for conducting PSA:
-
Stratification Treatment and comparison units are divided into strata (or subclasses) so that treated and comparison units are similar within each strata. Cochran (1968) observed that creating five subclassifications (stratum) removes at least 90% of the bias in the estimated treatment effect.
-
Matching - Each treatment unit is paired with a comparison unit based upon the pre-treatment covariates.
-
Weighting Each observation is weighted by the inverse of the probability of being in that group.
Stratification
Stratification involves dividing (or stratifying) the observations into subgroups based upon the propensity score. Here, we used quintiles on the propensity scores where were estimated using logistic regression. For classification trees the stratum is determined by the leaf nodes.
lr.out <- glm(treatment ~ x1 + x2 + x3, data = dat, family = binomial(link='logit')) dat$ps <- fitted(lr.out) # Get the propensity scores # Stratification breaks5 <- psa::get_strata_breaks(dat$ps) dat$strata5 <- cut(x = dat$ps, breaks = breaks5$breaks, include.lowest = TRUE, labels = breaks5$labels$strata)
ggplot(dat, aes(x = ps, color = as.logical(treatment))) + geom_density(aes(fill = as.logical(treatment)), alpha = 0.2) + geom_vline(xintercept = breaks5$breaks, alpha = 0.5) + geom_text(data = breaks5$labels, aes(x = xmid, y = 0, label = strata), color = 'black', vjust = 0.8, size = 8) + scale_fill_manual('Treatment', values = palette2) + scale_color_manual('Treatment', values = palette2) + xlab('Propensity Score') + ylab('Density') + xlim(c(0, 1)) + ggtitle('Density distribution of propensity scores by treatment', subtitle = 'Five strata represented by vertical lines')
Stratification (cont.)
Independent sample tests (e.g. t-tests) are conducted within each stratum and pooled to provide an overall estimate.
psa::stratification_plot(ps = dat$ps, treatment = dat$treatment, outcome = dat$outcome)
Matching
Dependent sample tests (e.g. t-tests) are conducted using match pairs to provide a treatment.
match_out <- Matching::Match(Y = dat$outcome, Tr = dat$treatment, X = dat$ps, caliper = 0.1, estimand = 'ATE') dat_match <- data.frame(treat_ps = dat[match_out$index.treated,]$ps, treat_outcome = dat[match_out$index.treated,]$outcome, control_ps = dat[match_out$index.control,]$ps, control_outcome = dat[match_out$index.control,]$outcome) psa::matching_plot(ps = dat$ps, treatment = dat$treatment, outcome = dat$outcome, index_treated = match_out$index.treated, index_control = match_out$index.control)
Matching Methods
There are many choices and approaches to matching, including:
- Propensity score matching.
- Limited exact matching.
- Full matching.
- Nearest neighbor matching.
- Optimal/Genetic matching.
- Mahalanobis distance matching (for quantitative covariates only).
- Matching with and without replacement.
- One-to-one or one-to-many matching.
Which method should you use?
Whichever one gives the best balance!
Weighting
Propensity score weights can be used as regression weights, the specific weights depend on the desired estimand and will be provided in later slides.
dat <- dat |> mutate( ate_weight = psa::calculate_ps_weights(treatment, ps, estimand = 'ATE'), att_weight = psa::calculate_ps_weights(treatment, ps, estimand = 'ATT'), atc_weight = psa::calculate_ps_weights(treatment, ps, estimand = 'ATC'), atm_weight = psa::calculate_ps_weights(treatment, ps, estimand = 'ATM') ) psa::weighting_plot(ps = dat$ps, treatment = dat$treatment, outcome = dat$outcome)
Shiny Application
We can explore how these three plots change as the treatment effects change using the psa::psa_simulation_shiny() application.
psa::psa_simulation_shiny()
knitr::include_graphics('images/psa_simulation_screenshot.png')
Phase I: Estimate Propensity Scores
.pull-left[ In this example we will use logistic regression to estimate the propensity scores.
lr.out <- glm( treatment ~ x1 + x2 + x3, data = dat, family = binomial(link='logit')) dat$ps <- fitted(lr.out) # Propensity scores
For stratification we will use quintiles to split the observations into five equal groups.
breaks5 <- psa::get_strata_breaks(dat$ps) dat$strata5 <- cut( x = dat$ps, breaks = breaks5$breaks, include.lowest = TRUE, labels = breaks5$labels$strata)
] .pull-right[.font70[
summary(lr.out)
]]
Distribution of Propensity Scores
ggplot(dat) + geom_histogram(data = dat[dat$treatment == 1,], aes(x = ps, y = after_stat(count)), bins = 50, fill = cols[2]) + geom_histogram(data = dat[dat$treatment == 0,], aes(x = ps, y = -after_stat(count)), bins = 50, fill = cols[1]) + geom_hline(yintercept = 0, lwd = 0.5) + scale_y_continuous(label = abs)
Check Balance: Multiple Covariates
PSAgraphics::cv.bal.psa(dat[,1:3], dat$treatment, dat$ps, strata = 5)
Check Balance: Single Covariate
.pull-left[
PSAgraphics::box.psa(dat$x1, dat$treatment, dat$strata5)
] .pull-right[
PSAgraphics::cat.psa(dat$x3, dat$treatment, dat$strata5)
]
PS Weights for Understanding Treatment Effects
Given that the distribution of treatment and control observations across the propensity score range are not the same, there are a number of alternative estimates of treatment effect. We will explore three additional esimates in addition to the classic average treatment effect.
dat <- dat |> mutate( ate_weight = psa::calculate_ps_weights(treatment, ps, estimand = 'ATE'), att_weight = psa::calculate_ps_weights(treatment, ps, estimand = 'ATT'), atc_weight = psa::calculate_ps_weights(treatment, ps, estimand = 'ATC'), atm_weight = psa::calculate_ps_weights(treatment, ps, estimand = 'ATM') ) dat |> head(n = 4)
Average Treatment Effect (ATE)
$$ATE = E(Y_1 - Y_0 | X) = E(Y_1|X) - E(Y_0|X)$$
ggplot() + geom_histogram(data = dat[dat$treatment == 1,], aes(x = ps, y = after_stat(count)), bins = 50, alpha = 0.5) + geom_histogram(data = dat[dat$treatment == 1,], aes(x = ps, weight = ate_weight, y = after_stat(count)), bins = 50, fill = cols[2], alpha = 0.5) + geom_histogram(data = dat[dat$treatment == 0,], aes(x = ps, y = -after_stat(count)), bins = 50, alpha = 0.5) + geom_histogram(data = dat[dat$treatment == 0,], aes(x = ps, weight = ate_weight, y = -after_stat(count)), bins = 50, fill = cols[1], alpha = 0.5) + ggtitle('Average Treatment Effect (ATE)')
Average Treatment Effect Among the Treated (ATT)
$$ATT=E(Y_{1}-Y_{0}|X,C=1)=E(Y_{1}|X,C=1)-E(Y_{0}|X,C=1)$$
ggplot() + geom_histogram(data = dat[dat$treatment == 1,], aes(x = ps, y = after_stat(count)), bins = 50, alpha = 0.5) + geom_histogram(data = dat[dat$treatment == 1,], aes(x = ps, weight = att_weight, y = after_stat(count)), bins = 50, fill = cols[2], alpha = 0.5) + geom_histogram(data = dat[dat$treatment == 0,], aes(x = ps, y = -after_stat(count)), bins = 50, alpha = 0.5) + geom_histogram(data = dat[dat$treatment == 0,], aes(x = ps, weight = att_weight, y = -after_stat(count)), bins = 50, fill = cols[1], alpha = 0.5) + ggtitle('Average Treatment Effect Among the Treated (ATT)')
Average Treatment Effect Among the Control (ATC)
$$ATC = E(Y_1 - Y_0 | X = 0) = E(Y_1 | X = 0) - E(Y_0 | X = 0)$$
ggplot() + geom_histogram(data = dat[dat$treatment == 1,], aes(x = ps, y = after_stat(count)), bins = 50, alpha = 0.5) + geom_histogram(data = dat[dat$treatment == 1,], aes(x = ps, weight = atc_weight, y = after_stat(count)), bins = 50, fill = cols[2], alpha = 0.5) + geom_histogram(data = dat[dat$treatment == 0,], aes(x = ps, y = -after_stat(count)), bins = 50, alpha = 0.5) + geom_histogram(data = dat[dat$treatment == 0,], aes(x = ps, weight = atc_weight, y = -after_stat(count)), bins = 50, fill = cols[1], alpha = 0.5) + ggtitle('Average Treatment Effect Among the Control (ATC)')
Average Treatment Effect Among the Evenly Matched
$$ATM_d = E(Y_1 - Y_0 | M_d = 1)$$
ggplot() + geom_histogram(data = dat[dat$treatment == 1,], aes(x = ps, y = after_stat(count)), bins = 50, alpha = 0.5) + geom_histogram(data = dat[dat$treatment == 1,], aes(x = ps, weight = atm_weight, y = after_stat(count)), bins = 50, fill = cols[2], alpha = 0.5) + geom_histogram(data = dat[dat$treatment == 0,], aes(x = ps, y = -after_stat(count)), bins = 50, alpha = 0.5) + geom_histogram(data = dat[dat$treatment == 0,], aes(x = ps, weight = atm_weight, y = -after_stat(count)), bins = 50, fill = cols[1], alpha = 0.5) + ggtitle('Average Treatment Effect Among the Evenly Matched (ACM)')
Treatment Effects for Weighting
$$Treatment\ Effect = \frac{\sum Y_{i}Z_{i}w_{i}}{\sum Z_{i} w_{i}} - \frac{\sum Y_{i}(1 - Z_{i}) w_{i}}{\sum (1 - Z_{i}) w_{i} }$$
Where $w$ is the weight (as defined in the following sections), $Z_i$ is the treatment assignment such that $Z = 1$ is treatment and $Z = 0$ is control, and $Y_i$ is the outcome
.pull-left[ $$w_{ATE} = \frac{Z_i}{\pi_i} + \frac{1 - Z_i}{1 - \pi_i}$$ $$w_{ATT} = \frac{\pi_i Z_i}{\pi_i} + \frac{\pi_i (1 - Z_i)}{1 - \pi_i}$$ ] .pull-right[ $$w_{ATC} = \frac{(1 - \pi_i) Z_i}{\pi_i} + \frac{(1 - e_i)(1 - Z_i)}{1 - \pi_i}$$ $$w_{ATM} = \frac{min{\pi_i, 1 - \pi_i}}{Z_i \pi_i (1 - Z_i)(1 - \pi_i)}$$ ]
class: font90
Treatment Effects
.pull-left[ .font80[Average Treatment Effect]
psa::treatment_effect( treatment = dat$treatment, outcome = dat$outcome, weights = dat$ate_weight)
lm(outcome ~ treatment, data = dat, weights = dat$ate_weight)
] .pull-right[ .font80[Average Treatment Effect Among the Treated]
psa::treatment_effect( treatment = dat$treatment, outcome = dat$outcome, weights = dat$att_weight)
lm(outcome ~ treatment, data = dat, weights = dat$att_weight)
]
class: font90
Treatment Effects (cont.)
.pull-left[ .font80[Average Treatment Effect Among the Control]
psa::treatment_effect( treatment = dat$treatment, outcome = dat$outcome, weights = dat$atc_weight)
lm(outcome ~ treatment, data = dat, weights = dat$atc_weight)
] .pull-right[ .font80[Average Treatment Effect Among the Evenly Matched]
psa::treatment_effect( treatment = dat$treatment, outcome = dat$outcome, weights = dat$atm_weight)
lm(outcome ~ treatment, data = dat, weights = dat$atm_weight)
]
class: inverse, middle, center
Example: National Supported Work Demonstration
National Supported Work
The National Supported Work (NSW) Demonstration was a federally and privately funded randomized experiment done in the 1970s to estimate the effects of a job training program for disadvantaged workers.
- Participants were randomly selected to participate in the training program.
- Both groups were followed up to determine the effect of the training on wages.
- Analysis of the mean differences (unbiased given randomization), was approximately $800.
Lalonde (1986) used data from the Panel Survey of Income Dynamics (PSID) and the Current Population Survey (CPS) to investigate whether non-experimental methods would result in similar results to the randomized experiment. He found results ranging from $700 to $16,000.
National Supported Work (cont.)
Dehejia and Wahba (1999) later used propensity score matching to analyze the data. The found that,
- Comparison groups selected by Lalonde were very dissimilar to the treated group.
- By restricting the comparison group to those that were similar to the treated group, they could replicate the original NSW results.
- Using the CPS data, the range of treatment effect was between $1,559 to $1,681. The experimental results for the sample sample was approximately $1,800.
The covariates available include: age, education level, high school degree, marital status, race, ethnicity, and earning sin 1974 and 1975.
Outcome of interest is earnings in 1978.
data(lalonde, package='Matching')
class: font90
Estimating Propensity Scores
.pull-left[ Estimate propensity scores using logistic regression.
lalonde.formu <- treat ~ age + educ + black + hisp + married + nodegr + re74 + re75 glm1 <- glm(lalonde.formu, data = lalonde, family = binomial(link = 'logit'))
Get the propensity scores:
lalonde$ps <- fitted(glm1)
Define the stratification:
strata5 <- cut(lalonde$ps, quantile(lalonde$ps, seq(0, 1, 1/5)), include.lowest = TRUE, labels = letters[1:5])
] .pull-right[ .font70[
summary(glm1)
] ]
Checking Balance: Covariate Balance Plot
covars <- all.vars(lalonde.formu) covars <- lalonde[,covars[2:length(covars)]] cv.bal.psa(covars, lalonde$treat, lalonde$ps, strata = 5)
Checking Balance: Continuous Covariates
.pull-left[
box.psa(lalonde$age, lalonde$treat, strata5)
] .pull-right[
box.psa(lalonde$re74, lalonde$treat, strata5)
]
Checking Balance: Continuous Covariates (cont.)
.pull-left[
box.psa(lalonde$educ, lalonde$treat, strata5)
] .pull-right[
box.psa(lalonde$re75, lalonde$treat, strata5)
]
Checking Balance: Categorical Covariates
.pull-left[
cat.psa(lalonde$married, lalonde$treat, strata5)
] .pull-right[
cat.psa(lalonde$hisp, lalonde$treat, strata5)
]
Checking Balance: Categorical Covariates (cont.)
.pull-left[
cat.psa(lalonde$black, lalonde$treat, strata5)
] .pull-right[
cat.psa(lalonde$nodegr, lalonde$treat, strata5)
]
Loess Regression
psadf <- data.frame(ps = lalonde$ps, Y = lalonde$re78, Tr = lalonde$treat) loess_plot(ps = psadf[psadf$Y < 30000,]$ps, outcome = psadf[psadf$Y < 30000,]$Y, treatment = as.logical(psadf[psadf$Y < 30000,]$Tr))
Stratification
.pull-left[
psa::stratification_plot(ps = psadf$ps, treatment = psadf$Tr, outcome = psadf$Y, n_strata = 5)
] .pull-right[
psa::stratification_plot(ps = psadf$ps, treatment = psadf$Tr, outcome = psadf$Y, n_strata = 10)
]
Stratification (cont.)
.pull-left[
strata5 <- cut(lalonde$ps, quantile(lalonde$ps, seq(0, 1, 1/5)), include.lowest = TRUE, labels = letters[1:5]) circ.psa(lalonde$re78, lalonde$treat, strata5)
]
.pull-right[
strata10 <- cut(lalonde$ps, quantile(lalonde$ps, seq(0, 1, 1/10)), include.lowest = TRUE, labels = letters[1:10]) circ.psa(lalonde$re78, lalonde$treat, strata10)
]
class: font60
Stratification (cont.)
.pull-left[
circ.psa(lalonde$re78, lalonde$treat, strata5)
]
.pull-right[
circ.psa(lalonde$re78, lalonde$treat, strata10)
]
Matching
rr <- Match(Y = lalonde$re78, Tr = lalonde$treat, X = lalonde$ps, M = 1, estimand = 'ATT', ties = FALSE) summary(rr)
Visualizing Matching Results
matches <- data.frame(Treat = lalonde[rr$index.treated,'re78'], Control = lalonde[rr$index.control,'re78']) granovagg.ds(matches[,c('Control','Treat')], xlab = 'Treat', ylab = 'Control')
Balance for Matching
psa::MatchBalance(df = lalonde, formu = lalonde.formu, formu.Y = update.formula(lalonde.formu, re78 ~ .), M = 1, estimand = 'ATT', ties = FALSE) |> plot()
Balance for Matching (cont.)
psa::MatchBalance(df = lalonde, formu = lalonde.formu, formu.Y = update.formula(lalonde.formu, re78 ~ .), exact.covs = c('nodegr'), #<< M = 1, estimand = 'ATT', ties = FALSE) |> plot()
class: inverse, middle, center
Sensitivity Analysis
Sensitivity Analysis
- An observational study is free of hidden bias if the propensity scores for each subject depend only on the observed covariates.
- That is, the p-value is valid if there are no unobserved confounders.
- However, there are very likely covariates that would better model treatment. These introduce hidden bias.
- Hidden bias exists if two subjects have the same covariates, but different propensity scores.
$X_a = X_b$ but ${ \pi }{ a }\neq { \pi }{ b }$ for some a and b.
Sensitivity Analysis
Each person in the treatment is matched to exactly one person in the control. The odds of being in the treatment for persons a and b are:
$O_a = \frac{ \pi_a }{ 1 - \pi_a }$ and $O_b = \frac{ \pi_b }{ 1 - \pi_b }$
The ratio of these odds, $\Gamma$, measures the bias after matching.
$$\Gamma =\frac { { O }{ a } }{ { O }{ b } } =\frac { { { \pi }{ a } / ( }{ 1-{ \pi }{ a }) } }{ { { \pi }{ b } / (1-{ \pi }{ b }) } }$$
This is the ratio of the odds the treated unit being in the treatment group to the matched control unit being in the treatment group.
Sensitivity Analysis
Sensitivity analysis tests whether the results hold for various ranges of $\Gamma$. That is, we test how large the differences in $\pi$ (i.e. propensity scores) would have to be to change our basic inference. Let $p_a$ and $p_b$ be the probability of each unit of the matched pair being treated, conditional on exactly one being treated. For example:
- If $\Gamma = 1$, the treatment and control unit within each pair has the same value of treatment assignment ( $p_a = 0.5$ and $p_b = 0.5$).
- If $\frac{1}{2} \le \Gamma \le 2$, no unit can be more than twice as likely as its match to get treated ( $0.33 \le p_a$, $p_b \le 0.66$).
- If $\frac{1}{3} \le \Gamma \le 3$, no unit can be more than three times as likely as its match to get treated ( $0.25 \le p_a$, $p_b \le 0.75$)
To get the bounds:
$$ \frac{1}{\Gamma +1 } \le p_a, p_b \le \frac{\Gamma}{\Gamma +1} $$
Wilcoxon Signed Rank Test
-
Drop pairs where the matches have the same outcome.
-
Calculate the difference in outcomes within each pair.
-
Rank the pairs from smallest absolute difference to largest absolute difference (i.e. the smallest = 1).
-
Take the sum of the ranks where the treated unit had the higher outcome.
$$W=\left| \sum { 1 }^{ { N }{ r } }{ sgn({ x }{ T,i }-{ x }{ C,i })\cdot { R }{ i } } \right|$$ Where $N$ is the number of ranked pairs; $R_i$ is the rank for pair r; $x{T,i}$ and $x_{C,i}$ are the outcomes for the $i^{th}$ treated and control pair, respectively.
Sensitivity Analysis
The process for sensitivity analysis:
-
Select a series of values for $\Gamma$. For social science research, values between 1 and 2 is an appropriate start.
-
For each $\Gamma$, estimate the p-values to see how the p-values increase for larger values of $\Gamma$.
-
For binary outcomes, use McNemar's test, for all others use Wilcoxon sign rank test and the Hodges-Lehmann point estimate. See Keele (2010) for more information.
Sensitivity Analysis
Children of parents who had worked in a factory where lead was used in making batteries were matched by age, exposure to traffic, and neighborhood with children whose parents did not work in lead-related industries. Whole blood was assessed for lead content yielding measurements in mg/dl
require(rbounds) psens(lalonde$re78[rr$index.treated], lalonde$re78[rr$index.control], Gamma = 2, GammaInc = 0.1)
class: inverse, middle, center
Bootstrapping Propensity Score Analysis
Boostrapping
-
Bootstrapping was first introduced by Efron (1979) in Bootstrap Methods: Another Look at the Jackknife.
-
Estimates confidence of statistics by resampling with replacement.
-
The bootstrap sample provides an estimate of the sampling distribution.
For PSA, sensitivity analysis is only well defined for matched samples.
Rosenbaum (2012) suggested that one way to test for sensitivity of model selection is to test the null hypothesis twice.
Bootstrapping Propensity Score Analysis
The PSAboot implements bootstrapping for propensity score analysis.
-
A stratified bootstrap sample is drawn to ensure the ratio of treatment-to-control observations is the same (i.e. sampling with replacement is done for the treatment and control observations is done separately). Note that the
control.ratioandtreat.ratioparameters allow for under sampling in the case of imbalanced data. -
For each bootstrap sample balance statistics and treatment effects are estimated using each method (five by default).
-
Overall treatment effect with confidence interval is estimated from the bootstrap samples.
PSA Bootstrapping Example
library(PSAboot) psaboot <- PSAboot(Tr = lalonde$treat, Y = lalonde$re78, X = lalonde, formu = lalonde.formu) summary(psaboot)
Checking Balance
psaboot_bal <- balance(psaboot) plot(psaboot_bal)
Plotting all Boostrap Samples
plot(psaboot)
Boxplot of Boostrapping PSA
boxplot(psaboot)
Comparing Across PSA Methods
matrixplot(psaboot)
class: inverse, middle, center
Matching of Non-Binary Treatments
Matching of Non-Binary Treatments
require(TriMatch) data(tutoring) formu <- ~ Gender + Ethnicity + Military + ESL + EdMother + EdFather + Age + Employment + Income + Transfer + GPA tutoring.tpsa <- trips(tutoring, tutoring$treat, formu) tutoring.matched.n <- trimatch(tutoring.tpsa, method=OneToN, M1=5, M2=3)
- The
TriMatchpackage provides functions for finding matched triplets. - Estimates propensity scores for three separate logistic regression models (one for each pair of groups, that is, treat1-to-control, treat2-to-control, and treat1-to-treat2).
- Finds matched triplets that minimize the total distance (i.e. sum of the standardized distance between propensity scores within the three models). within a caliper.
- Provides multiple methods for determining which matched triplets are retained:
- Optimal which attempts to retain all treatment units.
- Full which retains all matched triplets within the specified caliper (.25 by default as suggested by Rosenbaum).
- Analog of the one-to-many for matched triplets. Specify how many times each treat1 and treat2 unit can be matched.
- Unique which allows each unit to be matched once, and only once.
- Functions for conducting repeated measures ANOVA and Freidman Ranksum Tests are provided.
Example: Tutoring
Students can opt to utilize tutoring services to supplement math courses. Of those who used tutoring services, approximately 58% of students used the tutoring service once, whereas the remaining 42% used it more than once. Outcome of interest is course grade.
- Military Active military status.
- Income Income level.
- Employment Employment level.
- NativeEnglish Is English their native language
- EdLevelMother Education level of their mother.
- EdLevelFather Education level of their father.
- Ethnicity American Indian or Alaska Native, Asian, Black or African American, Hispanic, Native Hawaiian or Other Pacific Islander, Two or more races, Unknown, White
- Gender Male, Female
- Age Age at course start.
- GPA Student GPA at the beginning of the course.
New Student Outreach: Covariates
Newly enrolled students received outreach contacts until they registered for a course or six months have passed, whichever came first. Outreach was conducted by two academic advisors and a comparison group was drawn from students who enrolled prior to the start of the outreach program. Outcome of interest is number of credits attempted within the first seven months of enrollment.
PSA for Non-Binary Treatments
The TriMatch algorithm works as follows:
- Estimate three separate propensity score models for each pair of groups (i.e. Control-to-Treat1, Control-to-Treat2, Treat1-to-Treat2).
- Determine the matching order. The default is to start with the largest of two treatments, then the other treatment, followed by the control.
- For each unit in group 1, find all units from group 2 within a certain threshold (i.e. difference between PSs is within a specified caliper).
- For each unit in group 2, find all units from group 3 within a certain threshold.
- Calculate the distance (difference) between each unit 3 found and the original unit 1. Eliminate candidates that exceed the caliper.
- Calculate a total distance (sum of the three distances) and retain the smallest unique M group 1 units (by default M=2)
Matching Triplets
plot(tutoring.matched.n, rows=c(50), draw.segments=TRUE)
Checking Balance
multibalance.plot(tutoring.tpsa, grid=TRUE)
Results
boxdiff.plot(tutoring.matched.n, tutoring$Grade)
class: inverse, middle, center
Multilevel Propensity Score Analysis
Multilevel PSA
The use of PSA for clustered, or multilevel data, has been limited (Thoemmes \& Felix, 2011). Bryer and Pruzek (2012, 2013) have introduced an approach to analyzing multilevel or clustered data using stratification methods and implemented in the multilevelPSA R package.
- Exact and partially exact matching methods implicitly adjust for clustering. That is, the covariates chosen to exactly match are, in essence, clustering variables.
- Exact matching only applies to phase I of PSA. How are the clusters related to outcome of interest.
The multilevelPSA uses stratification methods (e.g. quintiles, classification trees) by:
- Estimate separate propensity scores for each cluster.
- Identify strata within each cluster (e.g. leaves of classification trees, quintiles).
- Estimate ATE (or ATT) within each cluster.
- Aggregate estimated ATE to provide an overall ATE estimate.
- Several functions to summarize and visualize results and check balance.
The Programme of International Student Assessment
-
International assessment conducted by the Organization for Economic Co-operation and Development (OECD).
-
Assesses students towards the end of secondary school (approximately 15-year-old children) in math, reading, and science.
-
Collects a robust set of background information from students, parents, teachers, and schools.
-
Assess both private and public school students in many countries.
-
We will use PISA to estimate the effects of private school attendance on PISA outcomes.
Phase I of Multilevel PSA
The multilevelPSA provides two functions, mlpsa.ctree and mlpsa.logistic, that will estimate propensity scores using classification trees and logistic regression, respectively. Since logistic regression requires a complete dataset (i.e. no missing values), we will use classification trees in this example.
data(pisana) data(pisa.colnames) data(pisa.psa.cols) student = pisana mlctree = mlpsa.ctree(student[,c('CNT','PUBPRIV',pisa.psa.cols)], formula=PUBPRIV ~ ., level2='CNT') student.party = getStrata(mlctree, student, level2='CNT') student.party$mathscore = apply( student.party[,paste0('PV', 1:5, 'MATH')], 1, sum) / 5
To assess what covariates were used in each tree model, as well as the relative importance, we can create a heat map of covariate usage by level.
Covariate Heat Map
tree.plot(mlctree, level2Col=student$CNT, colLabels=pisa.colnames[,c('Variable','ShortDesc')])
Phase II of Multilevel PSA
The mlpsa function will compare the outcome of interest.
results.psa.math = mlpsa(response=student.party$mathscore, treatment=student.party$PUBPRIV, strata=student.party$strata, level2=student.party$CNT, minN=5) results.psa.math$overall.wtd results.psa.math$overall.ci results.psa.math$level2.summary[,c('level2','Private','Private.n', 'Public','Public.n','diffwtd','ci.min','ci.max')]
Multilevel PSA Assessment Plot
The multilevel PSA assessment plot is an extension of the circ.psa plot in PSAgraphics introduced by Helmreich and Pruzek (2009).
plot(results.psa.math)
Multilevel PSA Difference Plot
mlpsa.difference.plot(results.psa.math, sd=mean(student.party$mathscore, na.rm=TRUE))
class: inverse, middle, center
Shiny Application
Shiny Application
psa::psa_shiny()
knitr::include_graphics('images/shiny_screenshot.png', dpi=350)
class: inverse, right, middle, hide-logo
Thank You!
r icons::fontawesome("paper-plane") jason.bryer@cuny.edu
r icons::fontawesome("github") @jbryer
r icons::fontawesome('mastodon') @jbryer@vis.social
r icons::fontawesome("link") psa.bryer.org
r icons::fontawesome("link") github.com/jbryer/psa
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.