Fit a Synthetic Control Method Model — scm

Unified formula interface for Synthetic Control and related causal inference methods. The formula syntax is:

Usage

scm_fit(
  formula,
  data,
  method = c("scm", "sdid", "gsc", "mc", "tasc", "si"),
  ...,
  predictors = NULL,
  covariates = NULL,
  v_selection = c("insample", "oos"),
  donor_mspe_threshold = Inf,
  lambda_pen = NULL,
  v_optim = c("auto", "coord_descent", "bfgs", "multistart"),
  qp_solver = c("auto", "active_set", "wolfe"),
  v_window = NULL,
  nu = NULL,
  fixedeff = NULL,
  pool_cohorts = TRUE,
  weights = NULL
)

Arguments

formula

A Formula object, e.g. y ~ D | unit + time.

data

A data.frame in long format (one row per unit-time).

method

One of "scm", "sdid", "gsc", "mc", "tasc", "si".

...

Estimator-specific arguments forwarded to the chosen method: lambda plus cv_folds, cv_nlambda and cv_seed ("mc"), r ("gsc", "tasc"), force ("gsc"), em_iter, em_tol, fix_A and intercept ("tasc"), zeta2 and zeta_t ("sdid"), k ("si"), control_group ("scm", "sdid", "gsc", "si"), and scale_predictors, cov_method and xi ("scm"). A name the chosen estimator does not take is an error rather than a silently ignored argument. Every argument after ... must be named in full: this is what lets lambda reach method = "mc" instead of partially matching lambda_pen.

control_group chooses each cohort's donors in a staggered fit of "scm", "sdid", "gsc" or "si". "never_treated" (default) uses the units never treated in the panel, the pure controls of Clarke et al. (2024, Algorithm 1); at least two are required. "clean" is deprecated: it adds the units adopting after the cohort, but each cohort's effect runs to the end of the panel, so those units' post-adoption outcomes enter its counterfactual and bias the effect towards zero. "mc" and "tasc" impute from every untreated cell and need no never-treated units.

cov_method chooses how covariates enter a staggered "scm" fit. "partial_out" (default) residualises the outcome on them first, the device of Clarke et al. (2024, S.2.2). "balance" is the formulation of Ben-Michael, Feller & Rothstein (2022, S.5.2, Equation 11): each covariate's pre-treatment average enters the partially pooled objective as a second imbalance measure, weighted against the lagged outcomes by xi. It requires the partially pooled solver, so it cannot be combined with nu = "separate", donor_mspe_threshold, lambda_pen or v_selection = "oos". xi defaults to the paper's rule, the sample variance of the donor units' outcomes before the last adoption date.

Three tuning parameters that govern model complexity are selected from the data when left at their NULL default, each by the rule its source paper specifies: r for "gsc" by the leave-one-period-out cross-validation of Xu (2017, Algorithm 1), lambda for "mc" by the K-fold cross-validation of Athey et al. (2021), and k for "si" by the optimal hard threshold of Gavish & Donoho (2014). Supplying a value skips the search. The r of "tasc" is a different quantity – the dimension of the state vector its Kalman smoother carries, not a factor count estimated from a donor block – and keeps its fixed default of 2.

force adds the additive fixed effects of the extended model in Xu (2017, Remark 3) to "gsc": "unit" a unit effect, "time" a time effect, "two-way" both, alongside a grand mean. "none" (default) fits the factors and covariates alone, the model of Xu's Equation (1), whose factors then have to absorb any unit or time level. The control units are within-transformed before the factors are extracted, and each treated unit's own effect is estimated with its loadings from its pre-treatment periods. With additive effects the factor count r may be zero, and its cross-validation starts there.

"tasc" runs the EM algorithm of Rho et al. (2026, Algorithm 7) until the fitted values settle, to a relative change of em_tol (default 1e-6) within at most em_iter iterations (default 2000), and warns if they do not. intercept = TRUE (default) gives each unit its own intercept in the observation equation. Without one the states have to carry each unit's level, and on panels far from zero the EM then drifts along directions the likelihood barely distinguishes, so the effect depends on when it stops; intercept = FALSE fits the paper's model as written.

predictors

A list() of pred() specifications that define the predictor matrix for SCM (see Abadie et al. 2010, S.2.3). Each pred() entry aggregates one or more variables over a time window. Pass NULL (default) to use all pre-treatment outcome periods as predictors; a one-line message states this default when it applies. A specification consisting solely of the outcome variable at each single pre-treatment period (one pred() per period, jointly covering the full pre-treatment window) defines the same predictor matrix and is therefore fitted through the same outcomes-only path as NULL, returning an identical fit; supply v_optim = "multistart" explicitly to force the predictor-path optimiser instead. That window is whatever data carries before the treatment date, not the span of the pred() list: adding earlier periods to data (to widen a plot, say) leaves the same list short of full coverage and moves the fit onto the predictor path, so a message reports the mismatch when it happens. Restricting the periods the fit is judged on is a separate choice from predictors – subset data, or keep data whole and pass v_window. pred() windows are expected to lie inside the pre-treatment periods present in data: times the panel does not carry, and times at or after the treatment date, each raise a warning because the row then aggregates something other than what the window reads as. Applies to method = "scm" only. Predictor rows are scaled by their standard deviation across all units before optimisation, matching the Synth reference implementation (ADH 2011, JSS); pass scale_predictors = FALSE to disable. When several units adopt at the same date they are pooled into one cohort, and each predictor row of the treated vector is their average under the same weights the outcome rows are pooled with, so both sides of the fit describe one target series.

covariates

An optional named list of additional time-varying covariates to partial out before estimation. Each element is a character string naming a column in data. Supported for method = "sdid", "gsc", and staggered method = "scm". Partialling out is the device of Clarke et al. (2024, S.2.2), which SDID needs because it has no way to balance a covariate. SCM does: Abadie, Diamond & Hainmueller (2010, S.2.3) put covariates into the V-weighted predictor distance, and Ben-Michael, Feller & Rothstein (2022, S.5.2) add a separate balance term for them. A sharp SCM fit therefore rejects covariates and asks for predictors = list(pred(...)) instead; the staggered path borrows the partial-out only because predictors is unavailable there.

v_selection

V matrix selection method for method = "scm". "insample" (default) follows Abadie et al. (2010): V is chosen by minimising in-sample pre-treatment MSPE. "oos" follows Abadie (2021) S.3.2 / ADH (2015): the pre-treatment window is split into a training half and a validation half. In the default outcomes-as-predictors case, candidate W(V) are fitted on training-half outcomes only, V* minimises the validation-half MSPE, and the final W* is refit with V* on the outcomes of the last floor(T_pre/2) pre-treatment periods (so v_weights has floor(T_pre/2) entries). With user-supplied predictors, the predictor matrix is fixed and only the MSPE evaluation window is restricted to the validation half; lag your pred() windows to the training period for a fully out-of-sample exercise.

donor_mspe_threshold

Donor pool filtering threshold (Abadie 2021 S.4). For method = "scm" only. Each donor's individual pre-treatment MSPE (using that donor alone as the counterfactual) is divided by the minimum such MSPE across all donors. Donors whose ratio exceeds this threshold are excluded from estimation. Inf (default) disables filtering.

lambda_pen

Penalised SCM parameter (Abadie & L'Hour 2021, JASA). For method = "scm" only. NULL (default) runs standard unpenalised SCM. "auto" selects the penalty via out-of-sample pre-treatment MSPE on the same validation window as v_selection = "oos". A non-negative number uses that value directly.

v_optim

Outer V-optimisation method for method = "scm". The outer problem (choose V so that the implied donor weights W(V) minimise pre-treatment outcome MSPE) is non-convex, and a single local search can settle in a poor basin when predictors are supplied. "auto" (default) therefore selects "multistart" for predictor-based fits and "coord_descent" for outcomes-only fits (where the single start is empirically reliable). "multistart" runs a deterministic multi-start search: a fixed start set (uniform V, one-hot V per predictor, and 100 fixed-seed random draws) is screened at one inner QP each, the leaders are polished by coordinate descent, and the winner is refined by a Nelder-Mead pass. Its solution is never worse (in pre-treatment loss) than "coord_descent", at roughly the cost of a handful of single-start fits, and is fully reproducible (no RNG state is consumed). "coord_descent" is the classic single-start coordinate descent with an 11-point grid, and is also the outcomes-only and staggered engine and the multi-start never-worse reference. "bfgs" (a single-start L-BFGS-B) is deprecated and will be removed in a future release: it has no advantage over "multistart" (which dominates it on predictor fits) or "coord_descent". mspe_ratio_pval() mirrors a multi-start fit in its placebo refits so the permutation test stays symmetric.

qp_solver

Inner-QP solver for method = "scm". "wolfe" is the Wolfe (1976) min-norm-point method, which exploits the fact that the inner QP is a projection onto the convex hull of the donor columns in the k-dimensional predictor space: by Caratheodory's theorem an optimum supported on at most k + 1 donors always exists, and Wolfe returns such a solution. "active_set" is the warm-started active-set method. Both return an exact optimum; they differ in which optimum when the QP is degenerate (fewer predictors than donors in the support), where the optimal set is a face rather than a point. "wolfe" is preferred there because its solution is sparse and reproducible, while "active_set" breaks those ties on round-off and can return different weights for the same fit under a different build. "auto" (default) is "wolfe" where that degeneracy can arise – sharp fits with predictors, and outcomes-only fits with fewer than half as many pre-treatment periods as donors – and "active_set" for the other outcomes-only fits, whose QP has a unique optimum in practice, so both solvers return the same weights and the active set gets there several times faster. It is also "active_set" on the paths that run their own solver: v_selection = "oos", lambda_pen, v_optim = "bfgs", and staggered adoption. Naming "wolfe" explicitly on one of those paths is an error rather than a silent downgrade. The resolved value is recorded in fit$qp_solver.

v_window

Optional vector of pre-treatment time values (matching the time index in data) over which the outer V optimisation evaluates the pre-treatment fit, for method = "scm" (sharp fits only). NULL (default) evaluates on all pre-treatment periods. The window restricts only the outer evaluation loss: predictor matrices (or the full pre-treatment outcome rows in the outcomes-only case) still enter the inner QP unchanged, and the reported loss and mspe_ratio_pval() MSPE components always cover the full pre-treatment window. Cannot be combined with v_selection = "oos", which manages its own train/validation split. Use it to judge the fit on part of the pre-treatment period while keeping the rest of data available for plotting and post-estimation.

nu

How strongly the adoption cohorts are pooled in staggered SCM fits (Ben-Michael, Feller & Rothstein 2022, JRSS-B). Cohort weight vectors are chosen jointly to minimise nu * (normalised pooled pre-treatment imbalance)^2 + (1 - nu) * (normalised per-cohort imbalance)^2, trading the fit of each cohort against the fit of their average, which is what the aggregate ATT rests on. "auto" (the default) is the paper's heuristic, the ratio of the pooled to the average per-cohort imbalance of the separate solution; it is also what augsynth, the authors' implementation, uses when nu is unset. A number in [0, 1] sets it directly, with 0 separate and 1 fully pooled – the paper advises against both extremes. "separate" is 0 reached through the older per-cohort path, which keeps the V optimisation and donor filtering that the pooled solver has no room for; it is also what the default falls back to when donor_mspe_threshold, lambda_pen, or v_selection = "oos" rules the pooled solver out. Balance diagnostics are stored in fit$pooling.

Staggered panels only. Under a single adoption date there are no cohorts to pool across, and several units sharing one date are already fully pooled into that cohort: S.4.1 of the paper shows the across-unit variation that would justify doing otherwise is zero when the adoption date is shared.

fixedeff

Unpenalised additive level terms. For method = "mc" it estimates unit and time fixed effects outside the nuclear norm, as in Athey et al. (2021, S.8.1), and defaults to TRUE; FALSE recovers the S.4 baseline where the level is absorbed into the penalised matrix. For method = "scm" it defaults to TRUE on staggered panels, matching augsynth and the S.6 finding that the intercept-shifted estimator "dominates either of the alternatives"; it does not apply to sharp fits. If TRUE, staggered SCM demeans every unit by its own pre-treatment mean within each cohort before fitting (intercept shift; Ben-Michael, Feller & Rothstein 2022, Section 5.1), which turns the estimator into a weighted difference-in-differences and typically improves fit when outcome levels differ across units. The reported Y_synth is shifted back to the raw outcome scale. Works with both the default and the partially pooled path. For method = "scm" on staggered panels only.

pool_cohorts

For staggered method = "scm": whether treated units sharing an adoption date are averaged into one block before the weights are fitted. TRUE (default) is the Appendix A.2 hierarchy of Ben-Michael, Feller & Rothstein (2022), whose S.4.1 shows there is no across-unit variation to exploit within a cohort. FALSE is their S.4 formulation, one weight vector per treated unit, which is what the rank-based permutation test in scm_rank_pval() needs. augsynth makes the unit-level form its default; the paper's own appendix does not. It applies to both staggered paths – the partially pooled solver and nu = "separate" – so it still takes effect where donor_mspe_threshold, lambda_pen or v_selection = "oos" have sent the fit down the latter. Sharp fits have a single cohort and reject it.

weights

Optional named numeric vector over the treated units giving the mass each carries in the aggregate ATT, for method = "scm". NULL (default) weights them equally, which for a shared adoption date is the N_treated * T_post weighting used throughout. Abadie (2021, S.8) notes that with regions or countries a population- or GDP-weighted average effect is often the quantity of interest; supply those figures here. The mass is used consistently: a block's target series is the correspondingly weighted average of its treated units, that same mass is what the block counts for in the partially pooled objective (in place of the cohort size of Ben-Michael, Feller & Rothstein 2022, Appendix A.2) and in its T_post-multiplied aggregation weight, and scm_inference() resamples each unit at its own mass. Only the ratios matter – the vector is rescaled internally – so equal weights reproduce the unweighted fit exactly.

Value

An object of classes c("coresynth_<method>", "coresynth"). Fits with staggered adoption additionally inherit from "coresynth_staggered", and multi-arm SI fits from "coresynth_multiarm"; S3 methods such as tidy() and augment() dispatch on these subclasses. All methods return at minimum:

method: estimator name
estimate: average treatment effect (ATT)
times: time index vector
T_pre: number of pre-treatment periods
Y_treat: treated unit outcome series
gap: treatment effect series (Y_treat - counterfactual)

Staggered method = "scm" fits add cohort_estimates, cohort_fits, and the settings the fit resolved: nu (the number the partially pooled solver used, or "separate" for the per-cohort path), fixedeff, pool_cohorts and qp_solver.

Details

outcome ~ treatment | unit_id + time_id

Examples

# Synthetic balanced panel: 10 units over 20 periods, unit 1 treated
# after period 15.
set.seed(1)
panel <- expand.grid(unit = 1:10, year = 1:20)
panel$treated <- as.integer(panel$unit == 1 & panel$year > 15)
panel$gdp <- panel$unit + 0.5 * panel$year +
  rnorm(nrow(panel)) + 3 * panel$treated

fit <- scm_fit(gdp ~ treated | unit + year, data = panel, method = "sdid")
summary(fit)
#> === coresynth summary ===
#> Method : SDID 
#> Periods : T_pre = 15 | T_post = 5 
#> ATT estimate: 3.49975 
#> Unit weights (non-zero donors):
#>      2      3      4      5      6      8      9     10 
#> 0.0886 0.1557 0.0275 0.1072 0.0966 0.1351 0.1746 0.2147 

# \donttest{
# Visualise the estimated gap (requires ggplot2)
plot(fit, type = "gap")

# }