Basic Synthetic Control
Alex Marsella
2022-10-11
This method was originally developed to compare a region in Spain to less than 20 control regions. In my job market paper, I am really comparing 8 treated tracts to hundreds of untreated tracts. There are newer methods for this kind of data where we have many units and multiple treated (such as microsynth from Robbins and Davenport (2020)), but for the sake of class, I will be applying SCM in the “classic” way. That is to say, I will do a bit of data manipulation to treat my treated group as one unified unit. It should be noted that synthetic control is weaker when there are limited pre-treatment periods. If I want to use ACS 5-year variables as covariates, which I will in this exercise, I will be limited to 2010-2013 as my years to “train” the model. This would be scrutinized heavily if my paper rested solely on this method to claim causation.
library(curl)
download.file("https://raw.githubusercontent.com/alexmarsella/alexmarsella.github.io/master/assets/aggrework.rds", "aggrework.rds", method="curl")
agg <- readRDS("aggrework.rds")There are many ways to go about using my data in the “old-school” synthetic control method that is made for one treated unit. One way would be to aggregate my treated tracts into one averaged (not additive, or else it will be extremely unlike every other tract, being insanely large) unit. So first, I will do some fancy data processing to ultimately turn my 8 treated tracts into one tract called “Tract 0”.
library(dplyr)
library(sf) #this is loaded because I have geometry in `agg` , without it, you might get errors
agg %>%
filter(fullzone==1) %>%
group_by(year) %>%
summarize(vcap=mean(vcap), white=mean(white),black=mean(black),
hisp=mean(hisp),nofathercap=mean(nofathercap),male15to29cap=mean(male15to29cap),
lessthanhighcap=mean(lessthanhighcap),highcap=mean(highcap),somecollegecap=mean(somecollegecap),
bachelorcap=mean(bachelorcap),pcinc=mean(pcinc),nowork12cap=mean(nowork12cap)) -> agg2
agg2$tract <- 0
subset(agg, agg$fullzone==0) -> agg3
bind_rows(agg3,agg2) -> agg4
agg4 %>%
group_by(year, tract) %>%
summarize(vcap=mean(vcap), white=mean(white),black=mean(black),
hisp=mean(hisp),nofathercap=mean(nofathercap),male15to29cap=mean(male15to29cap),
lessthanhighcap=mean(lessthanhighcap),highcap=mean(highcap),somecollegecap=mean(somecollegecap),
bachelorcap=mean(bachelorcap),pcinc=mean(pcinc),nowork12cap=mean(nowork12cap)) -> vsumdatayear
balanced <- vsumdatayear %>%
group_by(tract) %>%
filter(n_distinct(year) == 10)
balanced <- as.data.frame(balanced)
balanced <- na.omit(balanced)Finally, we will use a more modern package called tidysynth as opposed to the classic Synth package to run this synthetic control. In your homework, you are using Synth. It should be noted that some of this output is going to be very messy because it is going to try to list 375 tracts in one graph on multiple occasions. Don’t worry too much about this as it’s just this way for the sake of the tutorial. You will note, however, that the Promise Zone passes the placebo test, having the third highest Post/Pre MSPE ratio, giving it a p-value of 0.008.
library(tidysynth)
balanced$tract <- ifelse(balanced$tract==0, "Promise Zone", balanced$tract)
balanced_out <-
balanced %>%
# initial the synthetic control object
synthetic_control(outcome = vcap, # outcome
unit = tract, # unit index in the panel data
time = year, # time index in the panel data
i_unit = "Promise Zone", # unit where the intervention occurred
i_time = 2014, # time period when the intervention occurred
generate_placebos=T # generate placebo synthetic controls (for inference)
) %>%
# Generate the aggregate predictors used to fit the weights
generate_predictor(time_window = 2010:2013,
white= mean(white, na.rm = T),
black= mean(black, na.rm = T),
hisp = mean(hisp, na.rm = T),
nofathercap = mean(nofathercap, na.rm = T),
male15to29cap = mean(male15to29cap, na.rm = T),
lessthanhighcap = mean(lessthanhighcap, na.rm = T),
somecollegecap = mean(somecollegecap, na.rm = T),
bachelorcap = mean(bachelorcap, na.rm = T),
highcap = mean(highcap , na.rm = T),
pcinc = mean(pcinc , na.rm = T)) %>%
# average violent crime in 2013, right before treatment begins.
generate_predictor(time_window = 2013,
vcap = mean(vcap, na.rm = T)) %>%
# Generate the fitted weights for the synthetic control
generate_weights(optimization_window = 2010:2013, # time to use in the optimization task
margin_ipop = .02,sigf_ipop = 3,bound_ipop = 3 # optimizer options
) %>%
# Generate the synthetic control
generate_control()
balanced_out %>% plot_trends()
balanced_out %>% plot_weights()
balanced_out %>% grab_balance_table()## # A tibble: 11 × 4
## variable `Promise Zone` `synthetic_Promise Zone` donor_sample
## <chr> <dbl> <dbl> <dbl>
## 1 bachelorcap 0.101 0.109 0.140
## 2 black 0.733 0.731 0.435
## 3 highcap 0.310 0.319 0.338
## 4 hisp 0.0321 0.0363 0.109
## 5 lessthanhighcap 0.230 0.231 0.196
## 6 male15to29cap 0.192 0.192 0.125
## 7 nofathercap 0.575 0.576 0.442
## 8 pcinc 13565. 14328. 23691.
## 9 somecollegecap 0.249 0.249 0.218
## 10 white 0.181 0.183 0.426
## 11 vcap 36.4 36.4 27.7balanced_out %>% plot_placebos()
balanced_out %>% plot_mspe_ratio()
balanced_out %>% grab_signficance() ## # A tibble: 368 × 8
## unit_name type pre_mspe post_mspe mspe_ratio rank fishers_ex…¹ z_score
## <chr> <chr> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
## 1 139 Donor 8.76 481. 54.9 1 0.00272 8.26
## 2 147 Donor 7.51 279. 37.2 2 0.00543 5.42
## 3 Promise Zone Treated 4.44 150. 33.9 3 0.00815 4.88
## 4 353.02 Donor 0.198 6.44 32.6 4 0.0109 4.68
## 5 242 Donor 6.75 199. 29.4 5 0.0136 4.17
## 6 170 Donor 10.2 274. 26.8 6 0.0163 3.74
## 7 290 Donor 0.617 15.0 24.4 7 0.0190 3.36
## 8 145 Donor 33.0 785. 23.8 8 0.0217 3.26
## 9 323 Donor 4.93 116. 23.5 9 0.0245 3.21
## 10 78 Donor 11.8 270. 22.9 10 0.0272 3.13
## # … with 358 more rows, and abbreviated variable name ¹fishers_exact_pvaluebalanced_out %>% grab_unit_weights() ## # A tibble: 367 × 2
## unit weight
## <chr> <dbl>
## 1 1 0.000157
## 2 10.01 0.0000980
## 3 10.02 0.000107
## 4 100 0.00210
## 5 101 0.00108
## 6 102 0.0392
## 7 103 0.00130
## 8 104 0.00877
## 9 11.01 0.000273
## 10 11.02 0.000232
## # … with 357 more rowsbalanced_out %>% grab_predictor_weights## # A tibble: 11 × 2
## variable weight
## <chr> <dbl>
## 1 bachelorcap 0.0124
## 2 black 0.172
## 3 highcap 0.0228
## 4 hisp 0.134
## 5 lessthanhighcap 0.0915
## 6 male15to29cap 0.115
## 7 nofathercap 0.0429
## 8 pcinc 0.0446
## 9 somecollegecap 0.0934
## 10 white 0.132
## 11 vcap 0.140It should also be noted that the pre-treatment fit is not great, and warrants more experimentation to find a better fit. I would not rest my causal inference on this synthetic control specification. In my paper, I actually opt for the more robust “Synthetic Difference-in-Differences”, but that is a bit beyond the scope of this course.
As an exercise on your own, perhaps you can think of a criterion on which we could truncate our control group and/or use an alternative set of covariates for the matching process. Or maybe, we the Urban Institute guide we examined did, experiment with different pre-treatment values of the outcome variable, such as the average, or the 2010 value, or the 2010 and 2013 value, etc.
In fact, let’s try it one more time, but this time we’ll use the pre-treatment average of violent crime instead of 2013 only.
balanced_out2 <-
balanced %>%
# initial the synthetic control object
synthetic_control(outcome = vcap, # outcome
unit = tract, # unit index in the panel data
time = year, # time index in the panel data
i_unit = "Promise Zone", # unit where the intervention occurred
i_time = 2014, # time period when the intervention occurred
generate_placebos=T # generate placebo synthetic controls (for inference)
) %>%
# Generate the aggregate predictors used to fit the weights
generate_predictor(time_window = 2010:2013,
white= mean(white, na.rm = T),
black= mean(black, na.rm = T),
hisp = mean(hisp, na.rm = T),
nofathercap = mean(nofathercap, na.rm = T),
male15to29cap = mean(male15to29cap, na.rm = T),
lessthanhighcap = mean(lessthanhighcap, na.rm = T),
somecollegecap = mean(somecollegecap, na.rm = T),
bachelorcap = mean(bachelorcap, na.rm = T),
highcap = mean(highcap , na.rm = T),
pcinc = mean(pcinc , na.rm = T)) %>%
# average violent crime in 2013, right before treatment begins.
generate_predictor(time_window = 2010:2013,
vsum = mean(vcap, na.rm = T)) %>%
# Generate the fitted weights for the synthetic control
generate_weights(optimization_window = 2010:2013, # time to use in the optimization task
margin_ipop = .02,sigf_ipop = 3,bound_ipop = 3 # optimizer options
) %>%
# Generate the synthetic control
generate_control()
balanced_out2 %>% plot_trends()
balanced_out2 %>% plot_weights()
balanced_out2 %>% grab_balance_table()## # A tibble: 11 × 4
## variable `Promise Zone` `synthetic_Promise Zone` donor_sample
## <chr> <dbl> <dbl> <dbl>
## 1 bachelorcap 0.101 0.103 0.140
## 2 black 0.733 0.729 0.435
## 3 highcap 0.310 0.315 0.338
## 4 hisp 0.0321 0.0385 0.109
## 5 lessthanhighcap 0.230 0.232 0.196
## 6 male15to29cap 0.192 0.191 0.125
## 7 nofathercap 0.575 0.576 0.442
## 8 pcinc 13565. 15014. 23691.
## 9 somecollegecap 0.249 0.251 0.218
## 10 white 0.181 0.182 0.426
## 11 vsum 39.0 39.0 30.3balanced_out2 %>% plot_placebos()
balanced_out2 %>% plot_mspe_ratio()
balanced_out2 %>% grab_signficance() ## # A tibble: 368 × 8
## unit_name type pre_mspe post_mspe mspe_ratio rank fishers_ex…¹ z_score
## <chr> <chr> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
## 1 242 Donor 3.26 207. 63.4 1 0.00272 8.92
## 2 125 Donor 1.43 76.2 53.1 2 0.00543 7.38
## 3 78 Donor 4.05 213. 52.6 3 0.00815 7.29
## 4 145 Donor 18.9 598. 31.7 4 0.0109 4.16
## 5 147 Donor 16.4 401. 24.5 5 0.0136 3.09
## 6 267 Donor 2.91 68.5 23.5 6 0.0163 2.94
## 7 278 Donor 4.19 93.0 22.2 7 0.0190 2.74
## 8 Promise Zone Treated 7.24 146. 20.2 8 0.0217 2.45
## 9 9800 Donor 163. 3275. 20.1 9 0.0245 2.43
## 10 122.03 Donor 36.7 733. 20.0 10 0.0272 2.41
## # … with 358 more rows, and abbreviated variable name ¹fishers_exact_pvaluebalanced_out2 %>% grab_unit_weights() ## # A tibble: 367 × 2
## unit weight
## <chr> <dbl>
## 1 1 0.000179
## 2 10.01 0.000155
## 3 10.02 0.000166
## 4 100 0.00160
## 5 101 0.00114
## 6 102 0.0273
## 7 103 0.00134
## 8 104 0.00329
## 9 11.01 0.000360
## 10 11.02 0.000286
## # … with 357 more rowsbalanced_out2 %>% grab_predictor_weights## # A tibble: 11 × 2
## variable weight
## <chr> <dbl>
## 1 bachelorcap 0.124
## 2 black 0.173
## 3 highcap 0.0616
## 4 hisp 0.0826
## 5 lessthanhighcap 0.132
## 6 male15to29cap 0.0247
## 7 nofathercap 0.0479
## 8 pcinc 0.0331
## 9 somecollegecap 0.0499
## 10 white 0.195
## 11 vsum 0.0754The pre-treatment fit is actually worse for this one! Generally speaking, I would not rest my causal inference on either of these synthetic control specifications, the fits just aren't very good. In my experience, convincing people of causation with synthetic control does not always rest solely on the p-values, but of the p-values conditional on a solid pre-treatment match between the eventually treated unit and its synthetic partner.