knitr::opts_chunk$set( fig.width = 7, fig.height = 4.5, fig.align = "center", # cache = TRUE, autodep = TRUE )
Based on the basque
dataset and the corresponding example of package Synth
, this vignette illustrates the usage of package MSCMT
.
Thereby, parts of @Abadie2003 are reproduced in this example.
The MSCMT
package expects its input data to be a list
of matrices, where each matrix corresponds to one variable of interest.
The column names of the matrices correspond to the names of all units (treated unit and control units) and have to be identical across all elements of the list.
The row names correspond to the instants of time where data is available, they may well be different between variables.
Currently, all variables must either contain annual, quarterly, monthly, weekly or daily data.
The row names must therefore match one of the following patterns (where d
corresponds to a single digit):
dddd
for annual dates, e.g. 2016
,ddddQd
for quarterly dates, e.g. 2016Q1
for the first quarter of 2016,dddd-dd
for monthly dates, e.g. 2016-03
for March 2016.ddddWdd
for weekly dates, e.g. 2016W03
for the third week of 2016.dddd-dd-dd
for daily dates, e.g. 2016-03-17
.A list of matrices is a very useful, but quite unusual data format.
Therefore, a helper function called listFromLong
, which generates a list of matrices from
a more common data representation, the so-called 'long' format, has been included in package MSCMT
.
With listFromLong
, migrating from the Synth
package to MSCMT
is simple, because package Synth
(more precisely, its dataprep
function) expects its input data to be in
'long' format.
The basque
dataset in package Synth
in 'long' format has the following structure:
library(Synth) data(basque) colnames(basque) basque[703,]
With listFromLong
from package MSCMT
, the conversion to 'list' format is a simple task.
The names of the parameters correspond to the names of the parameters for function dataprep
of package Synth
:
library(MSCMT) Basque <- listFromLong(basque, unit.variable="regionno", time.variable="year", unit.names.variable="regionname") names(Basque) head(Basque$gdpcap)
The 'list' representation allows for simple, reproducible data preparations. To reproduce the analysis of @Abadie2003, the following preparations are sufficient:
# define the sum of all cases school.sum <- with(Basque,colSums(school.illit + school.prim + school.med + school.high + school.post.high)) # combine school.high and school.post.high in a single class Basque$school.higher <- Basque$school.high + Basque$school.post.high # calculate ratios and multiply by number of observations to obtain percentages from totals for (item in c("school.illit", "school.prim", "school.med", "school.higher")) Basque[[item]] <- 6 * 100 * t(t(Basque[[item]]) / school.sum)
In package Synth
, data preparation and model formulation are combined in function dataprep
, whereas the estimation is separated from the model formulation.
The approach of package MSCMT
is different: in order to facilitate the estimation of (potentially many) different models based on the same dataset, the model formulation has been
moved to function mscmt
, which also does the model estimation.
To allow for the methodological extensions (Multivariate Synthetic Control Methods using Time-Series), the model syntax had to be made more flexible. A model formulation now consists of
treatment.identifier
),controls.identifier
),2 x L
-matrix (parameter times.dep
), where L
dependent variables,2 x K
-matrix (parameter times.pred
), whereK
predictor variables,K
aggregation functions (parameter agg.fns
) for the predictor variables.
If missing, all predictor variables are considered to be time series, which corresponds to the function name "id"
.
Whenever the result of an aggregation function has length exceeding 1, the resulting data is considered as a time series, too.Note that times.dep
and times.pred
may contain duplicate column names for further increased flexibility of the model specification.
The following model specification reproduces the model in @Abadie2003:
treatment.identifier <- "Basque Country (Pais Vasco)" controls.identifier <- setdiff(colnames(Basque[[1]]), c(treatment.identifier, "Spain (Espana)")) times.dep <- cbind("gdpcap" = c(1960,1969)) times.pred <- cbind("school.illit" = c(1964,1969), "school.prim" = c(1964,1969), "school.med" = c(1964,1969), "school.higher" = c(1964,1969), "invest" = c(1964,1969), "gdpcap" = c(1960,1969), "sec.agriculture" = c(1961,1969), "sec.energy" = c(1961,1969), "sec.industry" = c(1961,1969), "sec.construction" = c(1961,1969), "sec.services.venta" = c(1961,1969), "sec.services.nonventa" = c(1961,1969), "popdens" = c(1969,1969)) agg.fns <- rep("mean", ncol(times.pred))
After preparing the data and formulating the model, the model estimation is done with function mscmt
.
Apart from the function parameters concerning data and model specification, there are further parameters which, e.g.,
A simple estimation (without a placebo study) produces the following console output (if parameter verbose
is TRUE
, which is the default):
res <- mscmt(Basque, treatment.identifier, controls.identifier, times.dep, times.pred, agg.fns, seed=1)
Package MSCMT
ships with an S3
method for print
, which gives a nice human-readable summary of the estimation results:
res
While there is a basic S3
method for plot
, it is strongly recommended to use package ggplot2
and the corresponding S3
method for ggplot
contained in package MSCMT
.
With the results of a simple estimation, two types of plots are available, the first being a comparison of original and synthesized data for the treated unit.
The variable to be plotted can be selected with parameter what
, by default the (first) dependent variable is being chosen:
library(ggplot2) ggplot(res, type="comparison")
The second type of plot is a plot of the gaps, ie. the differences between original and synthesized data:
ggplot(res, type="gaps")
It is possible to plot several variables by providing a vector for argument what
. Pre-defined sets of variables named "dependents"
, "predictors"
, and "all"
can be selected with parameter what.set
.
ggplot(res, what=c("gdpcap","invest","school.higher","sec.energy"), type="comparison")
Placebo studies are performed by simply setting the function argument placebo
to TRUE
. By default, the original treated unit is not added to the donor pool of the (original) control units, but this can be changed with parameter placebo.with.treated
.
A remarkable speed-up (depending on the number of control units) can be achieved for placebo studies by making use of a cluster, which can be set up with the help of package parallel
.
The simplest form of a cluster is a local cluster, which makes the power of multi-core cpus available for the (lengthy) computations involved with placebo studies.
Setting up a local cluster is very easy, see the example below.
The argument cl
of function mscmt
can be used to specify the cluster to be used for placebo studies.
The only drawback of using a cluster for placebo studies is losing the verbose output from the individual (one for each unit) SCM estimations, which
includes a lack of progress information for the whole placebo study.
Nevertheless, the speed-up should compensate for this drawback in all applications where a placebo study is meaningful.
Although clusters should be shut down automatically when the corresponding (master) R
process is finished, function stopCluster
can (and should) be used to shut down the cluster manually.
In the following example, a (local) cluster with two nodes is used for the estimation.
library(parallel) cl <- makeCluster(2) resplacebo <- mscmt(Basque, treatment.identifier, controls.identifier, times.dep, times.pred, agg.fns, cl=cl, placebo=TRUE, seed=1) stopCluster(cl)
Object resplacebo
now contains single SCM estimations for each unit as well as aggregated information concerning original data, synthesized data, and gaps for all units.
The individual SCM estimations can be accessed separately (as list elements with names corresponding to the units' names).
With the following plot, one can inspect whether there is some effect for Catalonia:
ggplot(resplacebo[["Cataluna"]], type="comparison")
Several functions in package MSCMT
are able to make use of the results of the placebo study as a whole.
One example are so-called placebo plots, by setting the plot type to "placebo.gaps"
(the default for results of a placebo study):
ggplot(resplacebo)
For statistical inference based on the results of a placebo study, the literature has developed so-called 'placebo tests', which have similarities to permutation tests. Two of these are
Most often, not all control units can be synthesized with an acceptable fit in a placebo study, resulting in large pre-treatment gaps. Of course, large post-treatment gaps are expected for these units, but since these gaps are rather caused from lack of fit than from an existing treatment effect, excluding such units is strongly advisable while investigating the effect for the (original) treated unit.
Excluding control units with large pre-treatment errors is usually done by limiting the ratio of a control unit's pre-treatment (r)mspe to the treated unit's (r)mspe.
Control units with large pre-treatment errors can easily be excluded from placebo plots:
ggplot(resplacebo, exclude.ratio=5, ratio.type="mspe")
The p-values of per-period placebo tests can be calculated via function pvalue
or plotted via ggplot
with plot.type="p.values"
:
pvalue(resplacebo, exclude.ratio=5, ratio.type="mspe", alternative="less") ggplot(resplacebo, exclude.ratio=5, ratio.type="mspe", type="p.value", alternative="less")
Calculating the aggregated treatment effect and testing its significance can be done with function did
of package MSCMT
:
did(resplacebo, range.post=c(1970,1990), exclude.ratio=5, alternative="less")
In this vignette, the basic workflow of preparing the data, defining and estimating the model, and evaluating the results of a (simple) SCM application with package MSCMT
has been illustrated.
Many features and options of package MSCMT
remained untouched, because this would have led far beyond the scope of this simple example.
There are more specialized vignettes which illustrate
synth
of package Synth
,In future releases of this package, additional vignettes will (probably) be added to illustrate
ggplot
method.Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.