## Below we provide two examples of how to use synth
## First example using artifcial data:
## First, load the dataframe synth.data
## this dataframe has
## 8 units:
## 1 treated, 6 control (2,7,13,17,29,32,38), treated is no 7
## 3 predictors (X1, X2, X3)
## 21 time periods (1980 - 2000)
## a unit.names.variable column ("names"
## and an outcome variable column (Y)
## all columns have column names
data(synth.data)
## run dataprep using character strings
## imagine treatment took place in 1991.
## as predictors, we want the average values of
## X1, X2, X3 over 1984 to 1989.
## and we want to track our synthetic control unit
## to track the treated unit's 1980-1991 outcomes
## as well as possible.
## we want to set it up to plot outcomes from 1984-1996
dataprep.out <-
dataprep(
foo = synth.data,
predictors = c("X1", "X2", "X3"),
predictors.op = "mean",
dependent = "Y",
unit.variable = "unit.num",
time.variable = "year",
special.predictors = list(
list("Y", 1991, "mean"),
list("Y", 1985, "mean"),
list("Y", 1980, "mean")
),
treatment.identifier = 7,
controls.identifier = c(29, 2, 13, 17, 32, 38),
time.predictors.prior = c(1984:1989),
time.optimize.ssr = c(1984:1990),
unit.names.variable = "name",
time.plot = 1984:1996
)
## then we run the synth command to identify optimal
## weights-- the optimal way to aggregate information about all the
## control units in such a way as to create the best possiblesynthetic
## control unit for the treated.
synth.out <- synth(dataprep.out,method = "BFGS")
## then we want to summarize the results, both numerically and graphically
## there are two ways to summarzie the results
## we can either access the output from synth.out directly
## eg.
synth.out$solution.w
# contains the unit weights or
synth.out$solution.v
## contains the predictor weights.
## the output from synth opt can be flexibly combined with
## the output from dataprep to compute other quatities of interest
## for example, the preriod by preriod discrepancies between the
## treated unit and its synthetic control unit
## can be computed by typing
gaps<-dataprep.out$Y1plot-(
dataprep.out$Y0plot
%*%
synth.out$solution.w
) ; gaps
## the Synth package also offers three convenience functions to summarize results.
## to get summary tables for all information
## (V and W weights plus balance btw. treated and synthtic control) use the
## synth.tab() command
synth.tables <- synth.tab(dataprep.res = dataprep.out,synth.res = synth.out)
print(synth.tables)
## to get summary plots for outcome trajectories
## of the treated and the synthtic control unit use the
## path.plot() and the gaps.plot() commands
## plot in levels (treated and synthetic)
path.plot(dataprep.res = dataprep.out,synth.res = synth.out)
## plot the gaps (treated - synthetic)
gaps.plot(dataprep.res = dataprep.out,synth.res = synth.out)
## Second example: The economic impact of terrorism in the
## Basque country using data from Abadie and Gardeazabal (2003)
## see the paper for details
data(basque)
# dataprep: prepare data for synth
dataprep.out <-
dataprep(
foo = basque
,predictors = c(
"school.illit","school.prim","school.med","school.high","school.post.high"
,"invest"
)
,predictors.op = c("mean")
,dependent = c("gdpcap")
,unit.variable = c("regionno")
,time.variable = c("year")
,special.predictors = list(
list("gdpcap", 1960:1969, c("mean")),
list("sec.agriculture", seq(1961,1969,2),c("mean")),
list("sec.energy", seq(1961,1969,2),c("mean")),
list("sec.industry", seq(1961,1969,2),c("mean")),
list("sec.construction", seq(1961,1969,2),c("mean")),
list("sec.services.venta", seq(1961,1969,2),c("mean")),
list("sec.services.nonventa",seq(1961,1969,2),c("mean")),
list("popdens", 1969, c("mean"))
)
,treatment.identifier = 17
,controls.identifier = c(2:16,18)
,time.predictors.prior = c(1964:1969)
,time.optimize.ssr = c(1960:1969)
,unit.names.variable = c("regionname")
,time.plot = c(1955:1997)
)
# 1. combine highest and second highest schooling category and eliminate highest category
dataprep.out$X1["school.high",] <- dataprep.out$X1["school.high",] + dataprep.out$X1["school.post.high",]
dataprep.out$X1 <- as.matrix(dataprep.out$X1[-which(rownames(dataprep.out$X1)=="school.post.high"),])
dataprep.out$X0["school.high",] <- dataprep.out$X0["school.high",] + dataprep.out$X0["school.post.high",]
dataprep.out$X0 <- dataprep.out$X0[-which(rownames(dataprep.out$X0)=="school.post.high"),]
# 2. make total and compute shares for the schooling catgeories
lowest <- which(rownames(dataprep.out$X0)=="school.illit")
highest <- which(rownames(dataprep.out$X0)=="school.high")
dataprep.out$X1[lowest:highest,] <- (100 * dataprep.out$X1[lowest:highest,]) / sum(dataprep.out$X1[lowest:highest,])
dataprep.out$X0[lowest:highest,] <- 100 * scale(
dataprep.out$X0[lowest:highest,],
center=FALSE,
scale=colSums(dataprep.out$X0[lowest:highest,])
)
# run synth
synth.out <- synth(
data.prep.obj = dataprep.out,
method = "BFGS"
)
# Get result tables
synth.tables <- synth.tab(
dataprep.res = dataprep.out,
synth.res = synth.out
)
# results tables:
print(synth.tables)
# plot results:
# path
path.plot(synth.res = synth.out,
dataprep.res = dataprep.out,
Ylab = c("real per-capita GDP (1986 USD, thousand)"),
Xlab = c("year"),
Ylim = c(0,13),
Legend = c("Basque country","synthetic Basque country"),
)
## gaps
gaps.plot(synth.res = synth.out,
dataprep.res = dataprep.out,
Ylab = c("gap in real per-capita GDP (1986 USD, thousand)"),
Xlab = c("year"),
Ylim = c(-1.5,1.5),
)Run the code above in your browser using DataLab