1  Gsynth

This chapter demonstrates the default gsynth estimator—estimator = "gsynth", (Xu 2017)— which estimates latent factors from the control group only and then imputes counterfactuals for treated units. We cover three settings: block DID, staggered DID, and unbalanced panels, with a dedicated section on inference and one on the placebo test.

For plot customization, see Chapter 3. For IFE-EM and matrix completion (the soft-deprecated paths), see Chapter 2.

ImportantConduct diagnostics

GSC imputes treated counterfactuals by projecting treated-unit loadings onto a factor space estimated from controls. The estimator can return tight-looking confidence intervals even when (a) the treated loadings sit outside the convex hull of control loadings, so the imputation is an extrapolation, or (b) the model overfits, and the parallel-trends-on-factors assumption fails. A point estimate plus a CI is not enough. Two checks should accompany every gsynth fit:

• Loading-overlap (plot(out, type = "loading.overlap")): are treated points inside the shaded control hull? Outside-hull treated units flag extrapolation; report results with that caveat or drop those units.
• Placebo test (placeboTest = TRUE, placebo.period = c(...)): does the model recover near-zero effects on held-out pre-treatment periods? A rejected placebo is evidence against the valid assumptions, and the headline ATT should not be read at face value.

Both diagnostics are demonstrated below for the simulated data and again for the staggered turnout application.

library(gsynth)
library(panelView)
data(gsynth)

1.1 Block DID

We start with a simulated dataset: 5 treated units, 45 control units, and 30 time periods. Treatment begins at Period 21 for all treated units.

head(simdata)

1.1.1 Visualizing data

Before estimation, visualize the data structure with panelView.

panelview(Y ~ D, data = simdata, index = c("id", "time"),
          pre.post = TRUE)

Outcome trajectories by treatment status:

panelview(Y ~ D, data = simdata, index = c("id", "time"),
          type = "outcome")

1.1.2 Estimation

We estimate the model using the outcome \(Y\), the treatment indicator \(D\), and two covariates \(X_1\) and \(X_2\).

The force option specifies fixed effects: "none", "unit", "time", or "two-way". The number of factors r counts factors on top of these fixed effects, so with force = "two-way", setting r to 0 gives the two-way fixed-effects model. Cross-validation (CV = TRUE) selects the number of factors within the range specified by c(0, 5). A single r with CV = TRUE searches from r to 5, so a call without r also searches 0 to 5. From v1.5.0 the CV default is rolling-window with k = 20 folds, cv.prop = 0.1, cv.buffer = 1. To reproduce the pre-v1.5.0 block-CV behavior in existing code, set cv.method = "block" and k = 5 explicitly. Cross-validation and the bootstrap use random draws; seed fixes them, so rerunning the code in this manual with the same parallel setting gives the same results (the number of cores does not change them).

NoteChanged in 1.5.0

From v1.3.1 to v1.4.0, a call with one number for r and CV = TRUE (the default) fitted r factors without a search, so a call without r fitted no factors. Now such a call searches from r to 5. Default calls can therefore give different results than in v1.4.0. To fit exactly r factors, set CV = FALSE.

out <- gsynth(Y ~ D + X1 + X2, data = simdata,
              index = c("id", "time"), force = "two-way",
              CV = TRUE, r = c(0, 5), se = TRUE,
              inference = "parametric", nboots = 1000,
              parallel = TRUE, cores = 16, seed = 02139)

The cross-validation results are stored in out$CV.out, with one row for each number of factors tried. It reports sigma2 (error variance), IC (information criterion), and MSPE (mean squared prediction error), among other measures. out$r.cv is the number of factors chosen.

out$CV.out[, c("r", "sigma2", "IC", "MSPE")]
#>      r    sigma2       IC     MSPE
#> [1,] 0 1.8486501 1.020232 2.358837
#> [2,] 1 1.5154097 1.205881 2.340072
#> [3,] 2 0.9973663 1.161299 1.914622
#> [4,] 3 0.9466397 1.472162 2.827789
#> [5,] 4 0.8941062 1.767452 3.244184
#> [6,] 5 0.8506026 2.059278 4.332982
out$r.cv
#> [1] 2

Cross-validation picks r by the rule in cv.rule. The default, "1se", picks the smallest r whose MSPE is within one standard error of the lowest MSPE.

out$est.att
#>              ATT      S.E.   CI.lower    CI.upper      p.value count
#> -19  0.392159788 0.5707540 -0.7264974  1.51081700 4.920255e-01     5
#> -18  0.276547958 0.4444736 -0.5946043  1.14770019 5.338156e-01     5
#> -17 -0.275392926 0.5128895 -1.2806379  0.72985201 5.913063e-01     5
#> -16  0.441201288 0.4279676 -0.3975998  1.28000236 3.025773e-01     5
#> -15 -0.889595124 0.4925567 -1.8549885  0.07579822 7.090636e-02     5
#> -14  0.593890957 0.4641080 -0.3157439  1.50352584 2.006719e-01     5
#> -13  0.528749012 0.4114603 -0.2776984  1.33519638 1.987732e-01     5
#> -12  0.171568737 0.5624257 -0.9307653  1.27390276 7.603270e-01     5
#> -11  0.610832288 0.4964509 -0.3621936  1.58385821 2.185480e-01     5
#> -10  0.170597468 0.4913591 -0.7924487  1.13364366 7.284448e-01     5
#> -9  -0.271891657 0.5707332 -1.3905081  0.84672480 6.337965e-01     5
#> -8   0.094842558 0.5072861 -0.8994199  1.08910502 8.516915e-01     5
#> -7  -0.651975701 0.5469977 -1.7240714  0.42012004 2.332939e-01     5
#> -6   0.573686472 0.4621732 -0.3321563  1.47952929 2.145022e-01     5
#> -5  -0.469685905 0.4945081 -1.4389040  0.49953216 3.422117e-01     5
#> -4  -0.077766449 0.5408511 -1.1378152  0.98228231 8.856700e-01     5
#> -3  -0.141784521 0.5337727 -1.1879598  0.90439075 7.905264e-01     5
#> -2  -0.157100323 0.3895251 -0.9205554  0.60635479 6.867184e-01     5
#> -1  -0.915575087 0.5220482 -1.9387707  0.10762050 7.946248e-02     5
#> 0   -0.003308833 0.3673270 -0.7232565  0.71663886 9.928129e-01     5
#> 1    1.235962010 0.6840227 -0.1046978  2.57662186 7.077755e-02     5
#> 2    1.630264312 0.5857685  0.4821792  2.77834947 5.383881e-03     5
#> 3    2.712177702 0.5726587  1.5897872  3.83456820 2.178539e-06     5
#> 4    3.466757691 0.7429665  2.0105701  4.92294531 3.069671e-06     5
#> 5    5.740132310 0.5192152  4.7224893  6.75777531 0.000000e+00     5
#> 6    5.280034526 0.5983469  4.1072962  6.45277282 0.000000e+00     5
#> 7    8.436484821 0.4493122  7.5558491  9.31712057 0.000000e+00     5
#> 8    7.839901526 0.6866684  6.4940561  9.18574696 0.000000e+00     5
#> 9    9.455114684 0.5350747  8.4063875 10.50384191 0.000000e+00     5
#> 10   9.638509457 0.4809624  8.6958404 10.58117847 0.000000e+00     5
out$est.avg
#>       ATT.avg      S.E. CI.lower CI.upper p.value
#> [1,] 5.543534 0.2625893 5.028868   6.0582       0
out$est.beta
#>        Coef       S.E.  CI.lower CI.upper p.value
#> X1 1.021890 0.02922963 0.9646012 1.079179       0
#> X2 3.052994 0.03019166 2.9938198 3.112169       0

1.1.3 Plot results

The default plot() produces a “gap” plot showing the estimated ATT by period. The true effects in simdata go from 1 to 10 (plus noise) from period 21 to 30.

plot(out)

For counterfactual plots, factor / loading plots, and customization (axes, legends, modern vs legacy visual recipe, highlight API), see Chapter 3.

1.1.4 Factors and loadings

We fit with two factors (r=2) to inspect the estimated factor structure directly:

out_r2 <- gsynth(Y ~ D + X1 + X2, data = simdata,
                 index = c("id", "time"), force = "two-way",
                 r = 2, CV = FALSE, se = FALSE,
                 parallel = FALSE)

type = "factors" plots the estimated time-varying factors, one curve per factor, and type = "loadings" shows the distribution of unit-level loadings on each factor. With force = "two-way", as here, both plots also show the fixed effects: the time effects as the curve labeled 0, and the unit effects in the panel labeled FE (include.FE = FALSE leaves them out):

plot(out_r2, type = "factors", xlab = "Time")

plot(out_r2, type = "loadings")

1.1.5 Loading-overlap diagnostic

The generalized synthetic control method imputes each treated unit’s counterfactual by projecting its loadings onto a factor space estimated from controls. When the treated-unit loadings sit outside the convex hull of control-unit loadings, the imputed counterfactual is an extrapolation rather than a weighted combination of observed control trajectories.

plot(out, type = "loading.overlap") shows treated and control loadings in the first two factor dimensions, with the convex hull of control loadings shaded. Use it as an at-fit diagnostic, paired with the placebo test below.

plot(out_r2, type = "loading.overlap")

When all treated points fall inside the shaded control hull, the GSC counterfactual is supported by interpolation across observed control trajectories. When treated points fall outside the hull, the counterfactual extrapolates beyond the observed support, and the placebo test below should be checked carefully.

1.2 Inference

gsynth() supports three inference methods via the inference argument; the default depends on the estimator:

inference =	Method	When to use	Default for
"parametric"	Parametric bootstrap with AR-aware residual draws (Xu (2017); AR-preserving fix in fect v2.3.0)	Small \(N_{tr}\); finite-sample setting	estimator = "gsynth"
"bootstrap" (alias "nonparametric")	Block bootstrap at unit level (or cluster level if cl set)	Large \(N\), \(N_{tr}\) growing with \(N\); super-population framework	estimator = "ife" and "mc"
"jackknife"	Leave-one-unit-out jackknife	Cluster-robust under leave-one-unit-out resampling; computationally cheap when \(N_{co}\) is moderate	—

1.2.1 Confidence-interval method (ci.method)

New in v1.5.0. The ci.method argument controls how the confidence intervals stored in the est.* slots (e.g. out$est.att) are formed:

• ci.method = "normal" (default): Wald intervals, \(\hat{\theta} \pm z_{1-\alpha/2}\,\widehat{\mathrm{SE}}\). Symmetric around the point estimate. Works with every inference, and it is the only interval for inference = "jackknife".
• ci.method = "basic": reflected-pivot bootstrap interval (Davison and Hinkley (1997), Sec. 5.2.1). Asymmetric; reads tail quantiles of the bootstrap distribution. Works with inference = "parametric" and "bootstrap". Its p-values agree with the basic interval (fect 2.4.7): under the bootstrap, a p-value is below alpha exactly when the interval excludes zero; with parametric inference, the p-values compare the estimate with the draws simulated under no effect. With inference = "jackknife", gsynth() stops, because the jackknife gives a standard error but no bootstrap distribution. With \(B < 1000\) replications the tail-quantile endpoints can be unstable; fect emits a warning and recommends nboots = 1000 for publication-grade CIs.

gsynth() accepts only these two values. For percentile, BC, and BCa intervals on the alternative estimands (att.cumu, aptt, log.att), call estimand(fit, type, ci.method = ...) post-fit. For a jackknife fit, estimand() also allows only ci.method = "normal", which is its default for such fits (fect 2.4.7).

1.2.2 Number of bootstrap replicates (nboots)

The default is nboots = 200, sufficient for SEs and Wald ("normal") CIs. For tail-quantile-based intervals (ci.method = "basic", "percentile", "bc", "bca"), bump to nboots = 1000+.

1.2.3 Parallel computing

parallel = TRUE (default) runs the bootstrap in parallel (and the cross-validation too, on large panels). Set cores explicitly when running on a shared machine; if cores = NULL, gsynth uses the number of available cores minus 2, at most 8. cores counts R worker processes, not CPU cores. With a multithreaded linear-algebra library (BLAS, such as OpenBLAS, MKL or Apple’s Accelerate), each worker may also start several threads, so the workers together can ask for more threads than the machine has cores. To avoid this, set the library’s thread count to 1 before starting R, for example OPENBLAS_NUM_THREADS=1 for OpenBLAS.

out_fast <- gsynth(Y ~ D + X1 + X2, data = simdata,
                   index = c("id", "time"), force = "two-way",
                   CV = TRUE, r = c(0, 5), se = TRUE,
                   inference = "parametric", nboots = 1000,
                   parallel = TRUE, cores = 16, seed = 02139)
print(out_fast)
#> Call:
#> gsynth(formula = Y ~ D + X1 + X2, data = simdata, index = c("id", 
#>     "time"), force = "two-way", r = c(0, 5), CV = TRUE, se = TRUE, 
#>     nboots = 1000, inference = "parametric", parallel = TRUE, 
#>     cores = 16, seed = 2139)
#> 
#> Average Treatment Effect on the Treated:
#>      ATT.avg   S.E. CI.lower CI.upper p.value
#> [1,]   5.544 0.2626    5.029    6.058       0
#> 
#>    ~ by Period (including Pre-treatment Periods):
#>           ATT   S.E. CI.lower CI.upper   p.value count
#> -19  0.392160 0.5708  -0.7265   1.5108 4.920e-01     5
#> -18  0.276548 0.4445  -0.5946   1.1477 5.338e-01     5
#> -17 -0.275393 0.5129  -1.2806   0.7299 5.913e-01     5
#> -16  0.441201 0.4280  -0.3976   1.2800 3.026e-01     5
#> -15 -0.889595 0.4926  -1.8550   0.0758 7.091e-02     5
#> -14  0.593891 0.4641  -0.3157   1.5035 2.007e-01     5
#> -13  0.528749 0.4115  -0.2777   1.3352 1.988e-01     5
#> -12  0.171569 0.5624  -0.9308   1.2739 7.603e-01     5
#> -11  0.610832 0.4965  -0.3622   1.5839 2.185e-01     5
#> -10  0.170597 0.4914  -0.7924   1.1336 7.284e-01     5
#> -9  -0.271892 0.5707  -1.3905   0.8467 6.338e-01     5
#> -8   0.094843 0.5073  -0.8994   1.0891 8.517e-01     5
#> -7  -0.651976 0.5470  -1.7241   0.4201 2.333e-01     5
#> -6   0.573686 0.4622  -0.3322   1.4795 2.145e-01     5
#> -5  -0.469686 0.4945  -1.4389   0.4995 3.422e-01     5
#> -4  -0.077766 0.5409  -1.1378   0.9823 8.857e-01     5
#> -3  -0.141785 0.5338  -1.1880   0.9044 7.905e-01     5
#> -2  -0.157100 0.3895  -0.9206   0.6064 6.867e-01     5
#> -1  -0.915575 0.5220  -1.9388   0.1076 7.946e-02     5
#> 0   -0.003309 0.3673  -0.7233   0.7166 9.928e-01     5
#> 1    1.235962 0.6840  -0.1047   2.5766 7.078e-02     5
#> 2    1.630264 0.5858   0.4822   2.7783 5.384e-03     5
#> 3    2.712178 0.5727   1.5898   3.8346 2.179e-06     5
#> 4    3.466758 0.7430   2.0106   4.9229 3.070e-06     5
#> 5    5.740132 0.5192   4.7225   6.7578 0.000e+00     5
#> 6    5.280035 0.5983   4.1073   6.4528 0.000e+00     5
#> 7    8.436485 0.4493   7.5558   9.3171 0.000e+00     5
#> 8    7.839902 0.6867   6.4941   9.1857 0.000e+00     5
#> 9    9.455115 0.5351   8.4064  10.5038 0.000e+00     5
#> 10   9.638509 0.4810   8.6958  10.5812 0.000e+00     5
#> 
#> Coefficients for the Covariates:
#>     Coef    S.E. CI.lower CI.upper p.value
#> X1 1.022 0.02923   0.9646    1.079       0
#> X2 3.053 0.03019   2.9938    3.112       0

1.2.4 Comparing CI methods on the same fit

out_basic <- gsynth(Y ~ D + X1 + X2, data = simdata,
                    index = c("id", "time"), force = "two-way",
                    CV = TRUE, r = c(0, 5), se = TRUE,
                    inference = "bootstrap", nboots = 1000,
                    ci.method = "basic",
                    parallel = TRUE, cores = 16, seed = 02139)
head(out_basic$est.att)
#>            ATT      S.E.   CI.lower   CI.upper    p.value count
#> -19  0.3921598 0.4322123 -0.4990681 1.18492018 0.37678030     5
#> -18  0.2765480 0.3132457 -0.3059942 0.91318466 0.33922796     5
#> -17 -0.2753929 0.4701334 -1.2007338 0.62484684 0.59933558     5
#> -16  0.4412013 0.2732394 -0.1734041 0.88720876 0.14032118     5
#> -15 -0.8895951 0.4442171 -1.6972307 0.01449076 0.05400718     5
#> -14  0.5938910 0.5898105 -0.5450310 1.74532814 0.35961463     5

The point estimate and SE are identical to a ci.method = "normal" fit on the same bootstrap replicates; the intervals and p-values differ. Asymmetry of the basic interval relative to the Wald interval signals skewness in the bootstrap distribution. With ci.method = "basic", the p-values go with the basic interval (fect 2.4.7): a p-value is below 0.05 exactly when the 95% basic interval excludes zero. At event time -15, for example, the basic interval includes zero and the p-value is above 0.05.

1.3 Placebo test

A placebo test re-fits the model after pretending some pre-treatment periods are part of the treated window, then asks whether the estimated effect on those held-out periods is close to zero. If the model is well-specified and the parallel-trends assumption holds, the placebo effect should be near zero. If it is not, the headline ATT estimate is suspect.

In v1.5.0, gsynth() accepts two new arguments for this:

• placeboTest = TRUE — turn the placebo test on.
• placebo.period = c(start, end) — which event-time periods to use as the placebo block. For example, c(-2, 0) treats event-times \(-2\) through \(0\) as placebo.

The default plot() highlights the placebo region automatically when placeboTest = TRUE was set at fit time.

out_pb <- gsynth(Y ~ D + X1 + X2, data = simdata,
                 index = c("id", "time"), force = "two-way",
                 r = 2, CV = FALSE, se = TRUE,
                 inference = "parametric", nboots = 500,
                 placeboTest = TRUE, placebo.period = c(-2, 0),
                 parallel = TRUE, cores = 16, seed = 02139)

plot(out_pb, main = "Placebo test on simdata")

The plot has two regions worth checking:

• Pre-placebo (event time \(< -2\)): the estimated effect should fluctuate around zero. Persistent deviation here means the model is fitting noise as signal, not an overlap problem.
• Placebo region (event times \(-2\) through \(0\)): the confidence interval should cover zero. If it does, the test does not reject pre-treatment zero, which is consistent with the parallel-trends assumption.

For numeric output by event time, the slot out_pb$est.placebo holds the placebo-region inference (point estimate, SE, CI), and out_pb$est.att holds the ATT estimates outside the placebo block.

out_pb$est.placebo
#>                      Coef      S.E.  CI.lower  CI.upper  p.value CI.lower(90%)
#> Placebo effect -0.4930715 0.3010401 -1.083099 0.0969563 0.101444    -0.9882385
#>                CI.upper(90%)
#> Placebo effect   0.002095436

A note on language: a placebo test that does not reject pre-treatment zero is consistent with the assumption; it does not prove the assumption holds. Conversely, a placebo test that rejects pre-treatment zero is evidence against the model’s pre-treatment fit, and the post-treatment ATT should then be reported with that caveat.

The placebo test comes back in Chapter 2 applied to the IFE-EM and matrix-completion estimators.

1.4 Staggered DID

The turnout dataset examines the effect of Election Day Registration (EDR) reforms on voter turnout in the United States, where treatment kicks in at different times.

names(turnout)
#> [1] "abb"            "year"           "turnout"        "policy_edr"    
#> [5] "policy_mail_in" "policy_motor"

1.4.1 Data structure

The balanced panel has 9 treated units with staggered treatment timing.

panelview(turnout ~ policy_edr, data = turnout,
          index = c("abb", "year"), pre.post = TRUE,
          by.timing = TRUE)

Outcome trajectories by treatment-status group:

panelview(turnout ~ policy_edr, data = turnout,
          index = c("abb", "year"), type = "outcome",
          main = "EDR Reform and Turnout",
          by.group = TRUE)

1.4.2 Estimation without factors

When no factor is assumed, the estimates are close to difference-in-differences:

out0 <- gsynth(turnout ~ policy_edr + policy_mail_in + policy_motor,
               data = turnout, index = c("abb", "year"),
               se = TRUE, inference = "parametric",
               r = 0, CV = FALSE, force = "two-way",
               nboots = 1000, seed = 02139)

plot(out0, type = "gap", xlim = c(-10, 5), ylim = c(-30, 20))

1.4.3 Estimation with factors

Cross-validation selects the number of factors:

out <- gsynth(turnout ~ policy_edr + policy_mail_in + policy_motor,
              data = turnout, index = c("abb", "year"),
              se = TRUE, inference = "parametric",
              r = c(0, 5), CV = TRUE, force = "two-way",
              nboots = 1000, seed = 02139)

1.4.4 Implied weights

out$wgt.implied (N_co × N_tr) stores the implied weights of all control units for each treated unit. Unlike classical synthetic control, the weights can be both positive and negative. It is NULL when cross-validation selects no factors (r.cv = 0).

NoteChanged in 1.5.0

With fect 2.4.7, wgt.implied uses the formula from Xu (2017). Weighting the control units’ outcomes, net of the fixed effects and covariates, gives the factor part of each treated unit’s counterfactual. The weights differ from those of earlier versions. Rows and columns are named by unit id. Here column 8 is Wisconsin ("WI").

dim(out$wgt.implied)
#> [1] 38  9
sort(out$wgt.implied[, 8])
#>            KY            FL            NM            AL            NC 
#> -0.0666150434 -0.0471286670 -0.0428120089 -0.0399577214 -0.0353501587 
#>            MD            TN            DE            VA            NV 
#> -0.0342555651 -0.0340305386 -0.0313750200 -0.0311529473 -0.0261010786 
#>            SC            MS            AR            AZ            IN 
#> -0.0214728836 -0.0183917692 -0.0173111436 -0.0150737627 -0.0147643827 
#>            GA            MO            IL            TX            OR 
#> -0.0106747218 -0.0089810911 -0.0035708332 -0.0015264494  0.0009464093 
#>            OH            WV            UT            MI            WA 
#>  0.0032483326  0.0057624783  0.0068496885  0.0115088734  0.0167074710 
#>            LA            CO            KS            OK            CA 
#>  0.0174455188  0.0184435129  0.0221358313  0.0230374294  0.0287200440 
#>            VT            NE            NJ            RI            NY 
#>  0.0345042870  0.0353460588  0.0385931797  0.0402606880  0.0429474652 
#>            SD            PA            MA 
#>  0.0469042286  0.0508528256  0.0563314639

1.4.5 Plot results

plot(out, type = "gap", xlim = c(-10, 5), ylim = c(-15, 15))

The gap plot always shows the average effect over all treated states. For one state, the counterfactual plot shows its observed outcome and its estimated counterfactual (Wisconsin):

plot(out, type = "counterfactual", id = "WI", main = "Wisconsin")

1.4.6 Factors and loadings

The estimated time-varying factors and unit-level loadings:

plot(out, type = "factors", xlab = "Year")

plot(out, type = "loadings")

Loading-overlap check on the turnout fit, in the first two factor dimensions:

plot(out, type = "loading.overlap")

Additional plot types and customization (raw, counterfactual, modern theme, highlight API) are in Chapter 3.

1.5 Unbalanced panels

gsynth can accommodate unbalanced panels. We demonstrate by randomly removing observations from the turnout dataset.

set.seed(123456)
turnout.ub <- turnout[-c(which(turnout$abb == "WY")[1:15],
                         sample(1:nrow(turnout), 50,
                                replace = FALSE)), ]

Visualize the data structure. White cells represent missing values.

panelview(turnout ~ policy_edr + policy_mail_in + policy_motor,
          data = turnout.ub, index = c("abb", "year"),
          pre.post = TRUE)

1.5.1 Estimation

out <- gsynth(turnout ~ policy_edr + policy_mail_in + policy_motor,
              data = turnout.ub, index = c("abb", "year"),
              se = TRUE, inference = "parametric",
              r = c(0, 5), CV = TRUE, force = "two-way",
              parallel = TRUE, min.T0 = 8,
              nboots = 1000, seed = 02139)

On unbalanced panels, the parametric bootstrap uses a unit-level Gaussian approximation for residual perturbations (introduced in fect v2.3.0); residuals are drawn from a Gaussian fitted to each unit’s residual structure rather than from the empirical pool, naturally accommodating missing cells. This is automatic — no API change for the user. Balanced panels see no change in behavior.

1.5.2 Missing-data plot

The "missing" plot is drawn by panelView from the data given to gsynth(). It shows each unit’s treatment status by period; white cells are missing observations. Units that gsynth() removed for having fewer than min.T0 pre-treatment periods still appear in the plot. The status used in estimation is in out$obs.missing: 1 = under treatment, 2 = not under treatment (control units, and treated units before treatment), 3 = missing, 4 = unit removed, 5 = placebo period.

plot(out, type = "missing")

Subset to treated units only:

plot(out, type = "missing", xlab = "Year", ylab = "State",
     main = "Treatment Status", id = out$id[out$tr],
     xlim = c(1920, 2012), axis.adjust = TRUE)

Gap plot with the unbalanced panel:

plot(out, type = "gap", xlim = c(-10, 5), ylim = c(-15, 15))

1.5.3 Notes on performance

1. Unbalanced panels take significantly more time than balanced panels (often 100:1 or higher) due to the EM algorithm for missing values.
2. Adding covariates also slows the algorithm because the IFE/MC model takes more iterations to converge.
3. Setting min.T0 to a positive number helps — the algorithm discards treated units with too few pre-treatment periods. A larger T_0 reduces bias and extrapolation risk. For cross-validation, min.T0 must be at least r.max + 2.

Athey, Susan, Mohsen Bayati, Nikolay Doudchenko, Guido Imbens, and Khashayar Khosravi. 2021. “Matrix Completion Methods for Causal Panel Data Models.” Journal of the American Statistical Association forthcoming.

Chen, Jiafeng, and Jonathan Roth. 2024. “Logs with Zeros? Some Problems and Solutions.” Quarterly Journal of Economics.

Davison, A. C., and D. V. Hinkley. 1997. Bootstrap Methods and Their Application. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press.

Gobillon, Laurent, and Thierry Magnac. 2016. “Regional Policy Evaluation: Interactive Fixed Effects and Synthetic Controls.” Review of Economics and Statistics 98 (3): 535–51.

Xu, Yiqing. 2017. “Generalized Synthetic Control Method: Causal Inference with Interactive Fixed Effects Models.” Political Analysis 25 (1): 57–76.