| Title: | Fast and Unified Synthetic Control Methods |
|---|---|
| Description: | A unified 'Formula' interface to the Synthetic Control Method (SCM) and related panel-data causal inference estimators: Synthetic Difference-in-Differences (SDID), Generalized Synthetic Control (GSC), Matrix Completion (MC), Time-Aware Synthetic Control (TASC), and Synthetic Interventions (SI), together with an experimental-design variant. Computational bottlenecks (quadratic programming, singular value decomposition, and Kalman filtering) are implemented in 'C++' via 'RcppArmadillo'. Methods are described in Abadie, Diamond and Hainmueller (2010) <doi:10.1198/jasa.2009.ap08746>, Arkhangelsky, Athey, Hirshberg, Imbens and Wager (2021) <doi:10.1257/aer.20190159>, Xu (2017) <doi:10.1017/pan.2016.2>, Athey, Bayati, Doudchenko, Imbens and Khosravi (2021) <doi:10.1080/01621459.2021.1891924>, and Agarwal, Shah and Shen (2025) <doi:10.1287/opre.2025.1590>. |
| Authors: | Yosuke Abe [aut, cre] |
| Maintainer: | Yosuke Abe <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 1.1.0.9000 |
| Built: | 2026-10-31 01:14:01 UTC |
| Source: | https://github.com/yo5uke/coresynth |
Applies a ridge-regression-based bias correction to a fitted SCM object, following Ben-Michael, Feller & Rothstein (2021, JASA). The corrected estimator is:
augment_scm(fit, lambda_ridge = NULL)augment_scm(fit, lambda_ridge = NULL)
fit |
A |
lambda_ridge |
Ridge penalty (non-negative). |
tau_aug = tau_SCM + (m_tr_post - sum_j W_j * m_j_post)
where m_i_post = Y_pre_i' beta_hat is the ridge outcome model prediction for unit i's mean post-treatment outcome, and beta_hat is estimated by ridge regression across control units.
A list with:
att_aug: Augmented ATT estimate
delta: Bias correction term (m_tr_post - sum_j W_j m_j_post)
att_scm: Original SCM ATT for comparison
lambda_ridge: Ridge penalty used
beta_hat: Ridge regression coefficients (length T_pre)
Implements the permutation-based conformal inference procedure of
Chernozhukov, Wuthrich & Zhu (2021, JASA). The test inverts a sharp null
by imputing the treated post-treatment
counterfactual as , re-estimating the counterfactual
proxy on all periods (imposing the null), and computing a
moving-block permutation p-value from the estimated residuals. A confidence
interval is obtained by test inversion over a grid of candidate .
conformal_inference( fit, tau0 = 0, q = 1, alternative = c("two.sided", "greater", "less"), ci = TRUE, level = 0.95, grid = NULL, n_grid = 200L, grid_mult = 4, ... )conformal_inference( fit, tau0 = 0, q = 1, alternative = c("two.sided", "greater", "less"), ci = TRUE, level = 0.95, grid = NULL, n_grid = 200L, grid_mult = 4, ... )
fit |
A |
tau0 |
Null value of the ATT for the reported p-value (default 0). |
q |
Exponent of the |
alternative |
|
ci |
Logical; construct a confidence interval by test inversion
(default |
level |
Confidence level for the interval (default 0.95). |
grid |
Optional numeric vector of candidate |
n_grid |
Number of grid points when |
grid_mult |
Half-width multiplier when |
... |
Unused. |
Supported for sharp (single-cohort) fits with method in
c("scm", "sdid", "gsc", "mc", "si"). Staggered, multi-arm, and tasc
fits are not supported (use sdid_inference(), gsc_inference(), or
si_inference() instead).
A list of class c("conformal_inference", "coresynth_inference")
with estimate, se (NA; conformal has no SE), p_value (at tau0),
ci_lower, ci_upper, method ("conformal"), n_controls,
alternative, staggered (FALSE), plus tau0, q, grid, and
p_grid (p-values along the grid). Compatible with tidy() / glance().
Chernozhukov, V., Wuthrich, K., & Zhu, Y. (2021). An Exact and Robust Conformal Inference Method for Counterfactual and Synthetic Controls. Journal of the American Statistical Association, 116(536), 1849-1864.
Generates a comprehensive, standardized JSON record covering all six
estimators. Suitable for reproducibility workflows (Xu & Yang 2026) and
downstream tooling. Pass the result of mspe_ratio_pval() or gsc_boot()
via the inference argument to include inference results.
export_json(x, file = "coresynth_results.json", inference = NULL, digits = 6L)export_json(x, file = "coresynth_results.json", inference = NULL, digits = 6L)
x |
A |
file |
Output file path. Default |
inference |
Optional list from |
digits |
Number of significant digits applied to numeric values (default 6L). |
Invisibly, the R list that was (or would be) serialized.
One-row summary of a coresynth_inference (or sdid_inference) object.
## S3 method for class 'coresynth_inference' glance(x, ...)## S3 method for class 'coresynth_inference' glance(x, ...)
x |
An inference object. |
... |
Unused. |
A one-row data.frame with columns method, n_controls,
staggered, estimate, std.error, p.value, conf.low,
conf.high, alternative, n_boot_valid.
Uncertainty for the ATT of a sharp GSC fit, conditional on the covariates, factors and loadings, by resampling the residuals of the estimated IFE model (Xu 2017, S.3.2 and Algorithm 2). It is the procedure the paper proposes when few units are treated, where resampling the treated units themselves cannot reach their noise.
gsc_boot(fit, B = 499L, alpha = 0.05, seed = NULL)gsc_boot(fit, B = 499L, alpha = 0.05, seed = NULL)
fit |
A |
B |
Number of bootstrap replications, and of prediction errors drawn for the treated units (default 499L). |
alpha |
Significance level for the confidence interval (default 0.05). |
seed |
RNG seed for reproducibility (default NULL). |
Each replicate rebuilds the panel without a treatment effect from the fitted values and re-estimates it. Residuals are resampled as whole time series, which keeps their serial correlation. The control units take the in-sample residuals of the IFE fit. The treated units take prediction errors instead, which are larger: before the bootstrap, a control unit is set aside as if it were treated, the model is refitted on a resample of the remaining controls, and the left-out unit's observed series minus its predicted counterfactual is one prediction error. Covariates enter through the fitted values and are re-estimated in each replicate.
A list with:
p_value: Two-sided p-value: mean(|ATT*| >= |ATT_obs|)
ci_lower: Lower bound of the (1-alpha)*100% CI for the ATT,
ATT_obs - q(1 - alpha/2) where q are quantiles of boot_dist
ci_upper: Upper bound, ATT_obs - q(alpha/2)
se: Bootstrap standard error
boot_dist: Numeric vector of length B (bootstrap ATT* values)
att_obs: Observed ATT from the original fit
Implements Xu (2017) IFE model with optional covariate adjustment and optional additive fixed effects (Xu 2017, Remark 3). When X_co has p > 0 slices, alternates the two steps of Bai (2009): factors: truncated SVD of Y_tilde = Y_co - X_co * beta beta: inverse(sum_i X_i' M_F X_i) sum_i X_i' M_F Y_i, with the loadings concentrated out through the annihilator M_F of the factors When X_co has 0 slices (default), falls back to the plain 3-step estimator. Additive effects are swept out of the control panel by the within transformation first; on a balanced panel this solves the joint least squares problem exactly, and the effects are recovered from the means of the residual.
gsc_ife_cpp( Y_co, Y_tr_pre, r, X_co, X_tr_pre, max_iter = 1000L, tol = 1e-08, force = "none" )gsc_ife_cpp( Y_co, Y_tr_pre, r, X_co, X_tr_pre, max_iter = 1000L, tol = 1e-08, force = "none" )
Y_co |
Control units outcome matrix (T x N_co) |
Y_tr_pre |
Treated units pre-treatment outcomes (T_pre x N_tr) |
r |
Number of latent factors (must be <= min(T, N_co)); 0 is allowed with additive effects |
X_co |
Time-varying covariate cube (T x N_co x p). Pass an empty cube (0 slices) for the covariate-free estimator. |
X_tr_pre |
Time-varying covariate cube for treated units in the pre-treatment window (T_pre x N_tr x p). Required for correct Step 2 loading estimation per Xu (2017): lambda_hat is estimated from Y_tr_pre - X_tr_pre * beta (covariate- demeaned). Pass an empty cube (0 slices) to skip demeaning (backward-compatible, but biased when beta != 0). |
max_iter |
Maximum iterations of the covariate loop (default 1000) |
tol |
Convergence tolerance on relative beta change (default 1e-8) |
force |
Additive fixed effects: |
A list with components:
F: estimated time factors (T x r).
L_co: control-unit factor loadings (N_co x r).
L_tr: treated-unit factor loadings (N_tr x r).
Y_tr_hat: estimated treated-unit counterfactual outcomes (T x N_tr),
additive effects included, covariate contribution excluded.
Y_co_hat: fitted control outcomes on the same basis (T x N_co).
singular_values: singular values from the final truncated SVD.
beta: estimated covariate coefficients (p x 1), empty when no
covariates are supplied.
mu, alpha_co, alpha_tr, xi: grand mean, unit effects of the
control and treated units, and time effects (zero where force
leaves them out).
n_iter: iterations of the covariate loop (0 without covariates).
converged: whether the covariate loop met tol within max_iter
(always TRUE without covariates).
Estimates SE and confidence intervals for the ATT by re-estimating the fit on resampled units: a bootstrap of the treated and the control units, drawn separately, or a jackknife that drops one unit at a time. Works for both sharp and staggered GSC fits.
gsc_inference( fit, method = c("bootstrap", "jackknife", "jackknife_global"), n_boot = 499L, level = 0.95, alternative = c("two.sided", "greater", "less"), seed = NULL )gsc_inference( fit, method = c("bootstrap", "jackknife", "jackknife_global"), n_boot = 499L, level = 0.95, alternative = c("two.sided", "greater", "less"), seed = NULL )
fit |
A |
method |
|
n_boot |
Number of bootstrap replications (default 499L; ignored for jackknife). |
level |
Confidence level (default 0.95). |
alternative |
|
seed |
RNG seed for reproducibility (default NULL). |
Xu (2017, S.3.2) notes that resampling units is valid when many units are
treated and proposes the parametric bootstrap of gsc_boot() for when few
are: the treated units' own noise varies across replicates only through
which treated units are drawn. With a single treated unit there is no such
variation, so this function stops and points to gsc_boot() and
conformal_inference().
On a staggered fit the bootstrap draws the treated units from all cohorts
together, so cohort sizes vary across replicates, and draws the controls
once per replicate for every cohort; cohorts are rebuilt from the draw and
combined with their N_treated * T_post weights. "jackknife" drops each
cohort's own donors and treated units and combines the per-cohort
variances by the delta method, which ignores the covariance that shared
donors induce between cohorts; "jackknife_global" drops every unit from
the whole panel in turn and keeps it. Both jackknives drop treated units,
so every cohort needs at least two (Clarke et al. 2024, Algorithm 3).
Note: gsc_boot() performs a parametric bootstrap under H0 (hypothesis
testing). gsc_inference() provides non-parametric SE and CIs suitable
for inference about the ATT magnitude.
A list of class coresynth_inference.
Implements the Kalman filter (forward pass) and Rauch-Tung-Striebel smoother (backward pass) for the state-space model in Rho et al. (2026):
kalman_smoother_cpp(Y, W, A, C, Q, R, z0, P0)kalman_smoother_cpp(Y, W, A, C, Q, R, z0, P0)
Y |
Observed data matrix (N x T). Use NA for unobserved entries. |
W |
Observation / loading matrix (N x r) |
A |
State transition matrix (r x r). Pass diag(r) for random-walk dynamics. |
C |
State drift vector (r x 1) |
Q |
State noise covariance (r x r) |
R |
Observation noise covariance (N x N, diagonal in practice) |
z0 |
Initial state mean (r x 1) |
P0 |
Initial state covariance (r x r) |
State: z(t+1) = A z(t) + C + eta(t), eta(t) ~ N(0, Q) Observation: y_t = W z_t + eps_t, eps_t ~ N(0, R)
Observation rows with NA (treated post-intervention) are automatically dropped at each time step so only control-unit rows update the filter.
The P update uses the numerically stable Joseph form: P(t|t) = (I - K W_obs) P(t|t-1) (I - K W_obs)^T + K R_obs K^T
A list with z_smooth, P_smooth, P_cross, z_pred, z_upd. P_cross is an r x r x (T-1) cube. Slice t (C++ 0-indexed, t=0,...,T-2) stores P(t+1, t | T) (0-indexed), i.e. P(t+2, t+1 | T) in 1-indexed Shumway-Stoffer notation. Formula: P(t+1|T) * J_t^T (eq. 6.68-6.69).
Iteratively re-estimates the synthetic control excluding one contributing donor at a time, holding the predictor weights V fixed at their baseline values (Abadie, Diamond & Hainmueller 2015, footnote 20). The spread of the leave-one-out ATT estimates shows how much the result hinges on any single donor.
loo_donors(fit, weight_threshold = 1e-06)loo_donors(fit, weight_threshold = 1e-06)
fit |
A sharp |
weight_threshold |
Only donors whose baseline weight exceeds this
value are dropped (removing a zero-weight donor cannot change the fit).
Default |
For penalised fits (lambda_pen used), the same penalty is re-applied in
each leave-one-out QP.
A list with:
att_original: baseline ATT
results: data.frame with one row per excluded donor
(donor, weight, att_loo)
att_range: range of the leave-one-out ATTs
placebo_in_time(), mspe_ratio_pval()
Computes the Abadie et al. (2010) / Abadie (2021) permutation p-value. For each control unit, a leave-one-out synthetic control is fitted.
mspe_ratio_pval( fit, mspe_threshold = 0, max_iter = 100L, tol = 1e-04, use_covariates = NULL, alternative = c("two.sided", "greater", "less"), mspe_prune = Inf, statistic = c("ratio", "post_mspe") )mspe_ratio_pval( fit, mspe_threshold = 0, max_iter = 100L, tol = 1e-04, use_covariates = NULL, alternative = c("two.sided", "greater", "less"), mspe_prune = Inf, statistic = c("ratio", "post_mspe") )
fit |
A |
mspe_threshold |
Minimum pre-treatment MSPE for including a control unit in the two-sided test. Ignored for one-sided tests. Default: 0 (no filtering). |
max_iter |
Passed to |
tol |
Passed to |
use_covariates |
Controls which predictor specification the placebo
refits use. Default |
alternative |
Direction of the alternative hypothesis:
|
mspe_prune |
Discard placebo runs whose pre-treatment MSPE exceeds
this multiple of the treated unit's, following Abadie, Diamond &
Hainmueller (2010, Figures 5-7), who use 20, 5 and 2. A run the synthetic
control could not fit before treatment carries no information about how
rare a large post-treatment gap is, so leaving it in the permutation
distribution only adds noise. |
statistic |
Which quantity to rank, for two-sided tests.
|
When alternative = "two.sided" (default), the test statistic is the
post/pre MSPE ratio, following Abadie et al. (2010). When
alternative = "greater" or "less", the test statistic is the signed
average post-treatment gap (ATT), giving a one-sided permutation test as
recommended by Abadie (2021) S.3.5 for improved power when the direction
of the treatment effect is known.
The placebo refits mirror the treated fit's outer optimiser and evaluation
window: a fit estimated with the multi-start outer search (v_optim = "multistart", or the "auto" default with predictors) or with a
v_window runs every placebo unit through the same configuration, keeping
the permutation statistic exchangeable across units.
An object of class scm_placebo (a list) with:
p_value: Permutation p-value between 0 and 1
mspe_ratio_treated: MSPE_post / MSPE_pre for the treated unit (two.sided only)
mspe_ratios_all: Named numeric vector (treated first, then controls); two.sided only
placebo_effects: Named N_co-vector of placebo ATT estimates
treated_effect: ATT estimate for the treated unit
n_placebo_used: Number of control units used
n_placebo_pruned: Number discarded by mspe_prune
mspe_prune, statistic: the settings the p-value was computed under
gaps: T x N_co matrix of placebo gap paths (unit minus its synthetic
control over all periods), for the Abadie et al. (2010) Figure 4-7 plot
treated_gap: T-vector of the treated unit's gap path
mspe_pre_treated, mspe_pre_placebo: Pre-treatment MSPEs used for
the relative pruning rule in plot.scm_placebo()
times, T_pre: Time axis metadata for plotting
plot.scm_placebo() for the placebo gap and MSPE ratio plots.
Re-estimates the synthetic control after artificially backdating the
treatment to a pre-treatment period, following Abadie, Diamond &
Hainmueller (2015) and Abadie & Vives-i-Bastida (2022, principle 7:
"out-of-sample validation is key"). Only pre-treatment data enter the
exercise, so the placebo gap after t0_placebo is uncontaminated by the
actual intervention. A credible design shows no sizable divergence at the
backdated treatment time.
placebo_in_time(fit, t0_placebo = NULL)placebo_in_time(fit, t0_placebo = NULL)
fit |
A sharp |
t0_placebo |
Backdated treatment period as a 1-based position in
|
The refit uses the outcomes of periods 1..t0_placebo as predictors
(the predictors = NULL default), regardless of how the original fit was
specified, because user-supplied pred() windows cannot be lagged
automatically (ADH 2015 lag their predictors by hand).
A list with:
t0_placebo: the backdated treatment period used
times: time values of the pre-treatment window
unit_weights: placebo donor weights
Y_treat, Y_synth, gap: series over the pre-treatment window
placebo_att: mean placebo gap over (t0_placebo, T_pre]
fit_rmspe: RMSPE over the placebo fitting window 1..t0_placebo
eval_rmspe: RMSPE over the placebo post window (t0_placebo, T_pre]
mspe_ratio_pval() for in-space placebos, loo_donors() for
donor-robustness checks.
Returns the tidy data.frame that plot() draws for a given type, so the
underlying series, weights, or placebo paths can be inspected, joined into a
table, or re-plotted directly. plot() stays the quick path; plot_data()
is the handle for anyone who wants to relabel simplified series names, feed
the numbers into their own figure, or postprocess them further.
plot_data(x, ...) ## Default S3 method: plot_data(x, ...) ## S3 method for class 'coresynth' plot_data( x, type = c("trend", "gap", "weights", "pred_weights"), align = FALSE, top_n = Inf, show_donors = 0, ... ) ## S3 method for class 'scm_placebo' plot_data(x, type = c("gaps", "ratios"), mspe_prune = NULL, ...) ## S3 method for class 'sdid_event_study' plot_data(x, ci_type = NULL, ...)plot_data(x, ...) ## Default S3 method: plot_data(x, ...) ## S3 method for class 'coresynth' plot_data( x, type = c("trend", "gap", "weights", "pred_weights"), align = FALSE, top_n = Inf, show_donors = 0, ... ) ## S3 method for class 'scm_placebo' plot_data(x, type = c("gaps", "ratios"), mspe_prune = NULL, ...) ## S3 method for class 'sdid_event_study' plot_data(x, ci_type = NULL, ...)
x |
A |
... |
Passed to methods (unused by the current methods). |
type |
For a |
align |
For |
top_n |
For |
show_donors |
For |
mspe_prune |
For a |
ci_type |
For a |
The frame mirrors what the matching plot(x, type = ...) call shows, with
two deliberate departures that make it a better data source:
Plain column names (time, value, series, weight, ...) are used
instead of the dotted convention of augment(), since this is data to
manipulate rather than model-augmented observations.
The cosmetic "drop donors with weight below 1e-4" filter that
plot(type = "weights") applies is not used here: every donor is
returned (use top_n to subset), so the frame is the complete set of
weights.
Only the arguments that change which rows or values appear are accepted
(align, top_n, show_donors, mspe_prune); purely cosmetic arguments
of plot() (colors, labels, linetypes, vline, fill, ...) have no
data counterpart and are not part of this interface.
A tidy data.frame. Columns by type:
"trend"time, value, series ("Treated" /
"Synthetic Control"); with show_donors > 0, also "Donors" rows
and a unit column.
"gap"time, gap (treated minus synthetic control).
"weights"unit, weight; SDID fits add a panel column
("omega" unit weights, "lambda" time weights), with unit holding
the pre-period label for "lambda" rows.
"pred_weights"predictor, weight (sharp SCM only).
"gaps"time, gap, unit (NA for the treated series),
series ("Treated" / "Placebo (donor pool)").
"ratios"unit, ratio, post_mspe, series – both
candidate statistics, so the frame answers either
mspe_ratio_pval() statistic.
plot.coresynth(), plot.scm_placebo()
set.seed(1) panel <- expand.grid(unit = 1:10, year = 1:20) panel$treated <- as.integer(panel$unit == 5 & 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 = "scm") head(plot_data(fit, type = "trend")) plot_data(fit, type = "gap") plot_data(fit, type = "weights") # Relabel the simplified series names, then plot it yourself df <- plot_data(fit, type = "trend") df$series <- sub("Synthetic Control", "Synthetic Unit 5", df$series) ggplot2::ggplot(df, ggplot2::aes(time, value, color = series)) + ggplot2::geom_line()set.seed(1) panel <- expand.grid(unit = 1:10, year = 1:20) panel$treated <- as.integer(panel$unit == 5 & 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 = "scm") head(plot_data(fit, type = "trend")) plot_data(fit, type = "gap") plot_data(fit, type = "weights") # Relabel the simplified series names, then plot it yourself df <- plot_data(fit, type = "trend") df$series <- sub("Synthetic Control", "Synthetic Unit 5", df$series) ggplot2::ggplot(df, ggplot2::aes(time, value, color = series)) + ggplot2::geom_line()
Plot a coresynth model
## S3 method for class 'coresynth' plot( x, type = c("trend", "gap", "weights", "pred_weights"), colors = NULL, labels = NULL, linetypes = NULL, vline = list(), vline_offset = .VLINE_OFFSET_DEFAULT, hline = list(), fill = NULL, top_n = Inf, align = FALSE, show_donors = 0, ... )## S3 method for class 'coresynth' plot( x, type = c("trend", "gap", "weights", "pred_weights"), colors = NULL, labels = NULL, linetypes = NULL, vline = list(), vline_offset = .VLINE_OFFSET_DEFAULT, hline = list(), fill = NULL, top_n = Inf, align = FALSE, show_donors = 0, ... )
x |
A |
type |
One of |
colors |
For |
labels |
For |
linetypes |
For |
vline |
Aesthetic overrides for the vertical treatment-time line, as a
list passed to |
vline_offset |
For |
hline |
Aesthetic overrides for the horizontal zero line in |
fill |
For |
top_n |
For |
align |
For |
show_donors |
For |
... |
Ignored. |
A ggplot2 plot object.
set.seed(1) panel <- expand.grid(unit = 1:10, year = 1:20) panel$treated <- as.integer(panel$unit == 5 & 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 = "scm") plot(fit, type = "trend") plot(fit, type = "gap") plot(fit, type = "weights") plot(fit, type = "weights", top_n = 5) # Predictor (V) weights: which predictors the fit leans on plot(fit, type = "pred_weights") # Overlay the five largest donors behind the treated/synthetic series plot(fit, type = "trend", show_donors = 5) # SDID: align the synthetic series on the lambda-weighted pre-period level fit_sdid <- scm_fit(gdp ~ treated | unit + year, data = panel, method = "sdid") plot(fit_sdid, type = "trend", align = TRUE) plot(fit_sdid, type = "gap", align = TRUE) # Customize series colors, legend text, line types, and reference lines plot(fit, type = "trend", colors = c(treated = "black"), labels = c(treated = "Unit 5"), vline = list(color = "red", linetype = "dashed")) # Draw both series solid, or restyle the gap line plot(fit, type = "trend", linetypes = c(synthetic = "solid")) plot(fit, type = "gap", linetypes = "dashed") # Move the treatment line onto the first post-treatment period, # pin it to an absolute time, or drop it entirely plot(fit, type = "trend", vline_offset = 0) plot(fit, type = "trend", vline = list(xintercept = 15.5)) plot(fit, type = "trend", vline = FALSE)set.seed(1) panel <- expand.grid(unit = 1:10, year = 1:20) panel$treated <- as.integer(panel$unit == 5 & 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 = "scm") plot(fit, type = "trend") plot(fit, type = "gap") plot(fit, type = "weights") plot(fit, type = "weights", top_n = 5) # Predictor (V) weights: which predictors the fit leans on plot(fit, type = "pred_weights") # Overlay the five largest donors behind the treated/synthetic series plot(fit, type = "trend", show_donors = 5) # SDID: align the synthetic series on the lambda-weighted pre-period level fit_sdid <- scm_fit(gdp ~ treated | unit + year, data = panel, method = "sdid") plot(fit_sdid, type = "trend", align = TRUE) plot(fit_sdid, type = "gap", align = TRUE) # Customize series colors, legend text, line types, and reference lines plot(fit, type = "trend", colors = c(treated = "black"), labels = c(treated = "Unit 5"), vline = list(color = "red", linetype = "dashed")) # Draw both series solid, or restyle the gap line plot(fit, type = "trend", linetypes = c(synthetic = "solid")) plot(fit, type = "gap", linetypes = "dashed") # Move the treatment line onto the first post-treatment period, # pin it to an absolute time, or drop it entirely plot(fit, type = "trend", vline_offset = 0) plot(fit, type = "trend", vline = list(xintercept = 15.5)) plot(fit, type = "trend", vline = FALSE)
Plot an scm_design object
## S3 method for class 'scm_design' plot(x, type = c("outcome", "gap"), ...)## S3 method for class 'scm_design' plot(x, type = c("outcome", "gap"), ...)
x |
An |
type |
|
... |
Currently ignored. |
A ggplot object: for type = "outcome", the synthetic treated and
synthetic control outcome series; for type = "gap", the estimated
treatment effect over the experimental periods with split-conformal
confidence intervals. The object is returned for printing or further
customisation.
Draws the trade-off traced by scm_balance_frontier(): the imbalance of the
treated average, which the aggregate ATT rests on, against the root mean
square of the per-block imbalances, which the block-level effects rest on.
Separate SCM sits at one end of the curve and fully pooled SCM at the other,
with the heuristic the fit used marked along it.
## S3 method for class 'scm_frontier' plot(x, relative = TRUE, label_nu = TRUE, colors = NULL, ...)## S3 method for class 'scm_frontier' plot(x, relative = TRUE, label_nu = TRUE, colors = NULL, ...)
x |
An |
relative |
Plot the imbalances normalised by the separate solution (default), which puts both axes on comparable scales, rather than in the outcome's own units. |
label_nu |
Annotate the drawn points with their |
colors |
Optional named overrides for the |
... |
Ignored. |
The curve is normally strongly convex, which is the point of looking at it:
a small move away from either end usually buys a large reduction in the
other imbalance for very little, so the interesting values of nu are the
ones near the bend.
A ggplot2 plot object.
set.seed(1) panel <- expand.grid(unit = 1:15, year = 1:20) adopt <- c(`1` = 12, `2` = 15, `3` = 17) panel$treated <- as.integer(!is.na(adopt[as.character(panel$unit)]) & panel$year >= adopt[as.character(panel$unit)]) panel$treated[is.na(panel$treated)] <- 0L panel$gdp <- panel$unit + 0.4 * panel$year + rnorm(nrow(panel)) + 2 * panel$treated fit <- scm_fit(gdp ~ treated | unit + year, data = panel, method = "scm") plot(scm_balance_frontier(fit, nu_grid = c(0, 0.25, 0.5, 0.75, 1)))set.seed(1) panel <- expand.grid(unit = 1:15, year = 1:20) adopt <- c(`1` = 12, `2` = 15, `3` = 17) panel$treated <- as.integer(!is.na(adopt[as.character(panel$unit)]) & panel$year >= adopt[as.character(panel$unit)]) panel$treated[is.na(panel$treated)] <- 0L panel$gdp <- panel$unit + 0.4 * panel$year + rnorm(nrow(panel)) + 2 * panel$treated fit <- scm_fit(gdp ~ treated | unit + year, data = panel, method = "scm") plot(scm_balance_frontier(fit, nu_grid = c(0, 0.25, 0.5, 0.75, 1)))
Visualizes the placebo study returned by mspe_ratio_pval(), following
Abadie, Diamond & Hainmueller (2010, Section 3.4).
## S3 method for class 'scm_placebo' plot( x, type = c("gaps", "ratios"), mspe_prune = NULL, colors = NULL, labels = NULL, linetypes = NULL, vline = list(), vline_offset = .VLINE_OFFSET_DEFAULT, hline = list(), ... )## S3 method for class 'scm_placebo' plot( x, type = c("gaps", "ratios"), mspe_prune = NULL, colors = NULL, labels = NULL, linetypes = NULL, vline = list(), vline_offset = .VLINE_OFFSET_DEFAULT, hline = list(), ... )
x |
A |
type |
One of |
mspe_prune |
Exclude placebo units whose
pre-treatment MSPE exceeds |
colors |
A named vector overriding series colors, e.g.
|
labels |
A named vector overriding the legend text of individual
series, e.g. |
linetypes |
Only for |
vline |
Only for |
vline_offset |
Only for |
hline |
Only for |
... |
Ignored. |
type = "gaps" overlays the treated unit's gap path (treated minus
synthetic control) on the placebo gap paths obtained by reassigning the
intervention to each donor unit (ADH 2010, Figure 4). Placebo units whose
synthetic control fits poorly before treatment carry no information about
the rarity of a large post-treatment gap, so ADH exclude units whose
pre-treatment MSPE exceeds a multiple of the treated unit's: 20, 5, and 2
in their Figures 5-7 (mspe_prune).
type = "ratios" shows the statistic behind the two-sided permutation
p-value, one point per unit (ADH 2010, Figure 8): the post/pre-treatment
MSPE ratio, or the post-treatment MSPE where the test was run with
statistic = "post_mspe". The ratio needs no pruning cutoff by
construction, which is why mspe_ratio_pval() leaves mspe_prune at
Inf; where the test did prune, the same runs are dropped here so that the
points and the p-value in the subtitle describe one set of placebo runs.
A ggplot2 plot object.
set.seed(1) panel <- expand.grid(unit = 1:10, year = 1:20) panel$treated <- as.integer(panel$unit == 5 & 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 = "scm") placebo <- mspe_ratio_pval(fit) # Treated gap overlaid on the donor-pool placebo gaps (ADH 2010, Fig. 4) plot(placebo, type = "gaps") # Prune poorly fitting placebos and relabel the legend plot(placebo, type = "gaps", mspe_prune = 5, labels = c(treated = "Unit 5")) # Set the line type of the placebo paths off against the treated one plot(placebo, type = "gaps", linetypes = c(placebo = "dotted")) # Move the treatment line onto the first post-treatment period plot(placebo, type = "gaps", vline_offset = 0) # Post/pre-treatment MSPE ratios (ADH 2010, Fig. 8) plot(placebo, type = "ratios")set.seed(1) panel <- expand.grid(unit = 1:10, year = 1:20) panel$treated <- as.integer(panel$unit == 5 & 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 = "scm") placebo <- mspe_ratio_pval(fit) # Treated gap overlaid on the donor-pool placebo gaps (ADH 2010, Fig. 4) plot(placebo, type = "gaps") # Prune poorly fitting placebos and relabel the legend plot(placebo, type = "gaps", mspe_prune = 5, labels = c(treated = "Unit 5")) # Set the line type of the placebo paths off against the treated one plot(placebo, type = "gaps", linetypes = c(placebo = "dotted")) # Move the treatment line onto the first post-treatment period plot(placebo, type = "gaps", vline_offset = 0) # Post/pre-treatment MSPE ratios (ADH 2010, Fig. 8) plot(placebo, type = "ratios")
Draws the period-by-period treatment effect estimates and bootstrap
confidence intervals returned by sdid_event_study(), together with a
horizontal zero reference line and a vertical treatment-time reference line.
When the numeric time index has gaps (e.g. missing survey years such as
2011–2012 between 2010 and 2013), axis tick marks default to the observed
periods (x$times) so unobserved years are never labeled on the axis, and
vline_offset interpolates across the adjacent observed periods (-0.5
places the vertical line midway between the last pre-treatment period and
the first post-treatment period; the default -1 places it on the last
pre-treatment period).
## S3 method for class 'sdid_event_study' plot( x, ci_type = NULL, colors = NULL, connect = FALSE, vline = list(), vline_offset = .VLINE_OFFSET_DEFAULT, hline = list(), ... )## S3 method for class 'sdid_event_study' plot( x, ci_type = NULL, colors = NULL, connect = FALSE, vline = list(), vline_offset = .VLINE_OFFSET_DEFAULT, hline = list(), ... )
x |
An |
ci_type |
Which confidence interval to draw: |
colors |
Optional single color string for the point estimates and
error bars (default |
connect |
Logical; connect adjacent period estimates with a line
(default |
vline |
Aesthetic overrides for the vertical treatment-time line, as a
list passed to |
vline_offset |
Where to draw the vertical treatment line, in periods
relative to the first post-treatment period (default |
hline |
Aesthetic overrides for the horizontal zero line, as a list
passed to |
... |
Ignored. |
A ggplot2 plot object.
sdid_event_study(), plot_data()
Creates a single predictor specification for use in scm_fit() with
method = "scm". Pass a list() of pred() calls as the predictors
argument to define the full covariate matrix.
pred(vars, times, op = "mean")pred(vars, times, op = "mean")
vars |
Character vector of variable names. All variables share the
same |
times |
Numeric/integer vector of time values to aggregate over.
Values are matched against the time index of the panel passed to
|
op |
Aggregation operator applied to each variable over |
A pred_spec object (a named list with class "pred_spec").
scm_fit() for the predictors argument that consumes a list()
of pred_spec objects.
# Three variables averaged over the same window pred(c("lnincome", "retprice", "age15to24"), 1980:1988) # Single variable at a specific year pred("cigsale", 1975) # Single variable averaged over a range pred("beer", 1984:1988) # Abadie, Diamond & Hainmueller (2010) California Prop 99 style: combine # several covariates aggregated over different windows plus three outcome # lags at specific years. The resulting list is passed to # scm_fit(..., predictors = predictors). predictors <- list( pred(c("lnincome", "retprice", "age15to24"), 1980:1988), pred("beer", 1984:1988), pred("cigsale", 1988), pred("cigsale", 1980), pred("cigsale", 1975) ) predictors# Three variables averaged over the same window pred(c("lnincome", "retprice", "age15to24"), 1980:1988) # Single variable at a specific year pred("cigsale", 1975) # Single variable averaged over a range pred("beer", 1984:1988) # Abadie, Diamond & Hainmueller (2010) California Prop 99 style: combine # several covariates aggregated over different windows plus three outcome # lags at specific years. The resulting list is passed to # scm_fit(..., predictors = predictors). predictors <- list( pred(c("lnincome", "retprice", "age15to24"), 1980:1988), pred("beer", 1984:1988), pred("cigsale", 1988), pred("cigsale", 1980), pred("cigsale", 1975) ) predictors
Traces how the two pre-treatment imbalances of partially pooled SCM trade
off as the pooling parameter nu varies: q_pool, the imbalance of the
average of the treated units, which the aggregate ATT rests on, and
q_sep, the root mean square of the per-block imbalances, which the
block-level effects rest on. Both are the cohort-level measures of
Ben-Michael, Feller & Rothstein (2022, Appendix A.2), in which a block
enters through the sum over its treated units, so a block of
units counts times towards q_pool and times towards
q_sep. Supplying weights to scm_fit() replaces the unit count with
the mass it gives them, and the frontier is traced on that same objective.
scm_balance_frontier(fit, nu_grid = seq(0, 1, length.out = 11L))scm_balance_frontier(fit, nu_grid = seq(0, 1, length.out = 11L))
fit |
A staggered |
nu_grid |
Values of |
Ben-Michael, Feller & Rothstein (2022) call this plot "an important tool for
using partially pooled SCM in practice" (S.4.2): the curve is strongly
convex, so a little pooling usually buys a large reduction in one imbalance
for a small rise in the other, and where that bargain stops being worth it
is a judgement the data can inform but not settle. Their nu heuristic,
which scm_fit() uses by default, is one point on this curve.
A fit made with cov_method = "balance" is traced on the objective that
includes its covariate term, and the frontier reports that term's two
measures in q_cov_sep and q_cov_pool. The plot still draws the outcome
imbalances, which are then only part of what the fit trades off.
A data.frame of class scm_frontier, one row per nu, with
nu, q_sep, q_pool, their values normalised by the separate solution
(q_sep_rel, q_pool_rel), estimate, and heuristic marking the row
closest to the fit's own nu. A fit made with cov_method = "balance"
trades off a covariate imbalance as well, so its frontier carries
q_cov_sep and q_cov_pool too, and the two normalisation constants are
then the combined measure of the paper's footnote 7 rather than the
outcome imbalance alone – which is why q_sep_rel and q_pool_rel are
not 1 at nu = 0 there. Plot it with plot().
Ben-Michael, E., Feller, A., & Rothstein, J. (2022). Synthetic controls with staggered adoption. JRSS-B, 84(2), 351-381.
set.seed(1) panel <- expand.grid(unit = 1:15, year = 1:20) adopt <- c(`1` = 12, `2` = 15, `3` = 17) panel$treated <- as.integer(!is.na(adopt[as.character(panel$unit)]) & panel$year >= adopt[as.character(panel$unit)]) panel$treated[is.na(panel$treated)] <- 0L panel$gdp <- panel$unit + 0.4 * panel$year + rnorm(nrow(panel)) + 2 * panel$treated fit <- scm_fit(gdp ~ treated | unit + year, data = panel, method = "scm") fr <- scm_balance_frontier(fit, nu_grid = c(0, 0.5, 1)) frset.seed(1) panel <- expand.grid(unit = 1:15, year = 1:20) adopt <- c(`1` = 12, `2` = 15, `3` = 17) panel$treated <- as.integer(!is.na(adopt[as.character(panel$unit)]) & panel$year >= adopt[as.character(panel$unit)]) panel$treated[is.na(panel$treated)] <- 0L panel$gdp <- panel$unit + 0.4 * panel$year + rnorm(nrow(panel)) + 2 * panel$treated fit <- scm_fit(gdp ~ treated | unit + year, data = panel, method = "scm") fr <- scm_balance_frontier(fit, nu_grid = c(0, 0.5, 1)) fr
Selects which units to assign to the treatment arm (and which to the control
arm) in a planned experiment, following Abadie and Zhao (2026). Both sets
of units are chosen by minimising the distance between their weighted-average
predictor vectors and the population-average predictor vector ,
so the resulting estimates are less susceptible to post-randomisation bias
than pure random assignment.
scm_design( data, outcome, unit, time, T0, T_fit = NULL, m_min = 1L, m_max = 1L, f = NULL, predictors = NULL, design = c("base", "weakly_targeted", "unit_level"), beta = 1, xi = 1, alpha = 0.05, normalize = TRUE, max_subsets = 100000L )scm_design( data, outcome, unit, time, T0, T_fit = NULL, m_min = 1L, m_max = 1L, f = NULL, predictors = NULL, design = c("base", "weakly_targeted", "unit_level"), beta = 1, xi = 1, alpha = 0.05, normalize = TRUE, max_subsets = 100000L )
data |
Long-format data frame (one row per unit–time). |
outcome |
Name of the outcome column. |
unit |
Name of the unit identifier column. |
time |
Name of the time identifier column. |
T0 |
Last pre-experimental period (a value present in the time
column). Periods after |
T_fit |
Number of fitting periods, counted from the start of
the pre-experimental phase. Defaults to |
m_min |
Minimum number of units assigned to treatment (default 1). |
m_max |
Maximum number of units assigned to treatment (default 1). |
f |
Named numeric vector of population weights |
predictors |
A |
design |
Design formulation: |
beta |
Trade-off parameter |
xi |
Trade-off parameter |
alpha |
Significance level for confidence intervals (default 0.05). |
normalize |
If |
max_subsets |
Maximum number of treatment-set candidates to evaluate before switching to random sampling (default 100 000). |
Three design formulations are available:
"base" (eq. 7): both the synthetic treated and the synthetic control
independently target the population average .
"weakly_targeted" (eq. 9): the treated and control weights jointly
minimise the distance of the synthetic treated to plus
beta times the distance between the synthetic control and the synthetic
treated.
"unit_level" (eq. 10): each treated unit gets its own synthetic
control; the treated weights trade off matching against xi
times each treated unit's own synthetic-control fit, and the aggregate
control weight is the w-weighted combination of the per-unit controls
(eq. 11).
Inference uses "blank periods" — pre-experimental periods whose outcomes were
not used to estimate the weights. Set T_fit strictly smaller than the
number of pre-experimental periods to enable the permutation test and split-
conformal confidence intervals from Section 3 of Abadie and Zhao (2026).
An object of class "scm_design" with components:
treated_units: unit identifiers selected for treatment
control_units: unit identifiers in the control pool
w: J-length weight vector for the synthetic treated unit (sums to 1)
v: J-length weight vector for the synthetic control unit (sums to 1)
tau_hat: estimated treatment effects for each experimental period
p_value: permutation p-value (NA when blank periods are unavailable)
ci_lower, ci_upper: per-period split-conformal confidence interval
Y_synth_tr, Y_synth_co: synthetic treated/control series (all periods)
estimate: ATT (mean of tau_hat)
Abadie, A. and Zhao, J. (2026). "Synthetic Controls for Experimental Design." MIT Working Paper.
Unified formula interface for Synthetic Control and related causal inference methods. The formula syntax is:
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 )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 )
formula |
A |
data |
A |
method |
One of |
... |
Estimator-specific arguments forwarded to the chosen method:
Three tuning parameters that govern model complexity are selected from
the data when left at their
|
predictors |
A |
covariates |
An optional named |
v_selection |
V matrix selection method for |
donor_mspe_threshold |
Donor pool filtering threshold (Abadie 2021 S.4).
For |
lambda_pen |
Penalised SCM parameter (Abadie & L'Hour 2021, JASA).
For |
v_optim |
Outer V-optimisation method for |
qp_solver |
Inner-QP solver for |
v_window |
Optional vector of pre-treatment time values (matching the
time index in |
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
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 |
pool_cohorts |
For staggered |
weights |
Optional named numeric vector over the treated units giving
the mass each carries in the aggregate ATT, for |
outcome ~ treatment | unit_id + time_id
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.
# 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) # Visualise the estimated gap (requires ggplot2) plot(fit, type = "gap")# 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) # Visualise the estimated gap (requires ggplot2) plot(fit, type = "gap")
Standard errors, confidence intervals and p-values for the aggregate ATT of a staggered SCM fit, by the two procedures of Ben-Michael, Feller & Rothstein (2022).
scm_inference( fit, method = c("jackknife", "wild_bootstrap"), n_boot = 1000L, level = 0.95, alternative = c("two.sided", "greater", "less"), seed = NULL )scm_inference( fit, method = c("jackknife", "wild_bootstrap"), n_boot = 1000L, level = 0.95, alternative = c("two.sided", "greater", "less"), seed = NULL )
fit |
A staggered SCM fit from |
method |
|
n_boot |
Number of bootstrap draws ( |
level |
Confidence level. Default 0.95. |
alternative |
Direction of the alternative hypothesis for the
p-value: |
seed |
Optional RNG seed ( |
method = "jackknife" (default) is the leave-one-unit-out jackknife of
their Appendix A.1: every unit in the panel – treated units and donors
alike – is dropped in turn, the weights (and intercepts, for fixedeff
fits) are refitted with the pooling parameter nu held at the fit's value,
and the aggregate ATT is recomputed. The variance is
, and the
interval and p-value use the normal approximation. Because the donors are
resampled too, the interval reflects the noise in the synthetic controls,
which every treated unit of a cohort shares.
method = "wild_bootstrap" is the multiplier bootstrap of their Section
5.3, which adapts Otsu & Rai (2017). The ATT is written as a sum over every
unit of the panel, and each unit's term is perturbed by its own
golden-ratio two-point multiplier while the weights stay fixed: a treated
unit's term is its effect's deviation from the ATT plus its residual, a
donor's is its residual times the weight the synthetic controls put on it.
Residuals are against a two-way (unit and period) model of the untreated
outcome. Following Otsu & Rai's Equation (9), only the treated units'
heterogeneity term is centred at the ATT; centring every unit, as the
paper's Equation (13) reads, would make the interval grow with the size of
the effect. The draws cost no refitting, so this is the choice for panels
too large for the jackknife. It is conservative where the outcome has
structure the two-way model does not absorb (unit-specific trends): in
simulation a nominal 95% interval covered 96-99.5% across donor pools of
8 to 40 and 3 to 24 treated units, and 93% with 2 treated units.
The jackknife refits the model once per unit, which grows quickly with the panel: about 0.2 s at 50 units, 0.7 s at 100 and 3 s at 200, but over a minute at 400. Donors that carry no weight in any block are not refitted (dropping one leaves the solution unchanged), and each refit starts from the full fit's weights.
Both work with both staggered SCM paths (the partially pooled default and
nu = "separate") and honour the aggregation weights N_treated x T_post,
or the weights given to scm_fit().
The jackknife treats the outcome as fixed: a fit with covariates
partialled out keeps the partialled-out outcome, and one with
cov_method = "balance" keeps its scaled covariate matrix, rather than
re-estimating either without the dropped unit. A dropped unit that was the
only member of its cohort removes that cohort from the estimate, as in the
paper's .
A coresynth_inference object with the standard fields
(estimate, se, p_value, ci_lower, ci_upper, method,
staggered, n_controls, alternative), compatible with
tidy.coresynth_inference() and glance.coresynth_inference().
n_treated records the number of treated units. The wild bootstrap adds
the draws in boot_ests, centred on the estimate; the jackknife adds the leave-one-out estimates
in jack_ests (named by the dropped unit, NA where a refit was not
possible).
Ben-Michael, E., Feller, A., & Rothstein, J. (2022). Synthetic controls with staggered adoption. JRSS-B, 84(2), 351-381.
Otsu, T., & Rai, Y. (2017). Bootstrap inference of matching estimators for average treatment effects. JASA, 112(520), 1720-1732.
set.seed(1) dat <- expand.grid(time = 1:20, id = paste0("u", 1:12)) dat$y <- rnorm(nrow(dat)) + as.numeric(factor(dat$id)) dat$d <- as.integer( (dat$id %in% c("u1", "u2") & dat$time > 10) | (dat$id %in% c("u3", "u4") & dat$time > 14) ) fit <- scm_fit(y ~ d | id + time, data = dat, method = "scm") scm_inference(fit)set.seed(1) dat <- expand.grid(time = 1:20, id = paste0("u", 1:12)) dat$y <- rnorm(nrow(dat)) + as.numeric(factor(dat$id)) dat$d <- as.integer( (dat$id %in% c("u1", "u2") & dat$time > 10) | (dat$id %in% c("u3", "u4") & dat$time > 14) ) fit <- scm_fit(y ~ d | id + time, data = dat, method = "scm") scm_inference(fit)
Solves the inner-loop QP for SCM: given a fixed diagonal metric matrix V, finds donor weights W on the simplex minimising the V-weighted covariate loss. The returned weights are a KKT-verified exact optimum whenever the active-set solver converges, with accelerated projected gradient as a fallback.
scm_inner_weights_cpp(X0, X1, V_diag, wolfe = FALSE)scm_inner_weights_cpp(X0, X1, V_diag, wolfe = FALSE)
X0 |
Covariate matrix for control units (k x N_co) |
X1 |
Covariate vector for the treated unit (k x 1) |
V_diag |
Diagonal of the metric matrix V (k x 1, non-negative, need not sum to 1) |
wolfe |
If |
Donor weight vector W (N_co x 1) on the unit simplex
For each control unit, treats it as pseudo-treated and fits SCM weights from the remaining N_co-1 donors. Returns MSPE components for constructing MSPE-ratio permutation p-values in R.
scm_placebo_cpp( Y_pre, Y_post, max_iter = 100L, tol = 1e-04, z_rows = NULL, wolfe = FALSE )scm_placebo_cpp( Y_pre, Y_post, max_iter = 100L, tol = 1e-04, z_rows = NULL, wolfe = FALSE )
Y_pre |
Control pre-treatment outcomes (T_pre x N_co) |
Y_post |
Control post-treatment outcomes (T_post x N_co) |
max_iter |
Outer coordinate-descent iterations (default 100) |
tol |
Convergence tolerance for V updates (default 1e-4) |
z_rows |
Optional 1-based pre-period row indices of the outer
evaluation window (the |
wolfe |
If |
A list with:
mspe_pre: N_co-vector of pre-treatment MSPE per placebo unit
mspe_post: N_co-vector of post-treatment MSPE per placebo unit
effects: N_co-vector of mean post-period gap per placebo unit
gaps: (T_pre + T_post) x N_co matrix of placebo gap paths
Covariate-spec counterpart of scm_placebo_cpp(): for each control unit,
treats it as pseudo-treated with its own predictor column X0[, i] and
fits the nested V/W optimisation against the remaining donors' predictors
X0[, -i], evaluating the prediction loss on pre-treatment outcomes.
Each leave-one-out problem is identical to a scm_weights_cpp() call on
the same submatrices; iterations are independent and run in parallel
under OpenMP.
scm_placebo_x_cpp( X0, Y_pre, Y_post, max_iter = 100L, tol = 1e-04, z_rows = NULL, multistart = FALSE, wolfe = FALSE )scm_placebo_x_cpp( X0, Y_pre, Y_post, max_iter = 100L, tol = 1e-04, z_rows = NULL, multistart = FALSE, wolfe = FALSE )
X0 |
Predictor matrix for control units (k x N_co), on the same
scale as the treated fit (SD-scaled when |
Y_pre |
Control pre-treatment outcomes (T_pre x N_co) |
Y_post |
Control post-treatment outcomes (T_post x N_co) |
max_iter |
Outer coordinate-descent iterations (default 100) |
tol |
Convergence tolerance for V updates (default 1e-4) |
z_rows |
Optional 1-based pre-period row indices of the outer
evaluation window (the |
multistart |
If |
wolfe |
If |
A list with:
mspe_pre: N_co-vector of pre-treatment MSPE per placebo unit
mspe_post: N_co-vector of post-treatment MSPE per placebo unit
effects: N_co-vector of mean post-period gap per placebo unit
gaps: (T_pre + T_post) x N_co matrix of placebo gap paths
A placebo unit whose solver fails yields NaN entries.
Permutation inference on the aggregate effect when more than one unit is treated, following Abadie & L'Hour (2021, S.3.2). Each treated unit contributes its own test statistic; the identities of the treated units are then reassigned at random among all units, and the sum of the ranks of the real treated units' statistics is compared with its permutation distribution.
scm_rank_pval( fit, B = 500L, statistic = c("mspe_ratio", "effect", "abs_effect"), alternative = c("two.sided", "greater", "less"), seed = 1L )scm_rank_pval( fit, B = 500L, statistic = c("mspe_ratio", "effect", "abs_effect"), alternative = c("two.sided", "greater", "less"), seed = 1L )
fit |
A |
B |
Number of treatment reassignments. |
statistic |
Per-unit statistic to rank. |
alternative |
|
seed |
Seed for the reassignment draw, so the p-value is reproducible without consuming the caller's RNG state. |
Ranking is what makes the test usable here. mspe_ratio_pval() compares one
treated unit against placebo runs on single donors, which breaks down as
soon as the treated series is an average: the average is less noisy than any
single donor, so its pre/post MSPE ratio is inflated and the test rejects
far too often. Ranks sidestep that, because every statistic entering the
comparison, real or permuted, is computed on a single unit.
The fit must expose per-unit effects, so a staggered fit has to be made with
pool_cohorts = FALSE. A sharp fit with several treated units pools them by
design, and its per-unit statistics are computed here from the stored
per-unit series.
An object of class c("scm_rank_test", "coresynth_inference") with
p_value, statistic (the observed rank sum), rank_sums (its
permutation distribution), unit_statistics, n_treated, B and
alternative.
Abadie, A., & L'Hour, J. (2021). A penalized synthetic control estimator for disaggregated data. JASA, 116(536), 1817-1834.
mspe_ratio_pval() for the single-treated-unit test.
set.seed(1) panel <- expand.grid(unit = 1:12, year = 1:20) panel$treated <- as.integer(panel$unit <= 3 & panel$year > 12) panel$gdp <- panel$unit + 0.4 * panel$year + rnorm(nrow(panel)) + 2 * panel$treated fit <- scm_fit(gdp ~ treated | unit + year, data = panel, method = "scm") scm_rank_pval(fit, B = 50)set.seed(1) panel <- expand.grid(unit = 1:12, year = 1:20) panel$treated <- as.integer(panel$unit <= 3 & panel$year > 12) panel$gdp <- panel$unit + 0.4 * panel$year + rnorm(nrow(panel)) + 2 * panel$treated fit <- scm_fit(gdp ~ treated | unit + year, data = panel, method = "scm") scm_rank_pval(fit, B = 50)
Jointly optimises donor weights W (on the simplex) and the diagonal metric matrix V via coordinate descent on the pre-treatment prediction MSPE, following Abadie, Diamond & Hainmueller (2010).
scm_weights_cpp( X0, X1, Z0, Z1, max_iter = 100L, tol = 1e-04, t_train = -1L, z_rows = NULL, multistart = FALSE, cheap_face = FALSE, wolfe = FALSE )scm_weights_cpp( X0, X1, Z0, Z1, max_iter = 100L, tol = 1e-04, t_train = -1L, z_rows = NULL, multistart = FALSE, cheap_face = FALSE, wolfe = FALSE )
X0 |
Covariate matrix for control units (k x N_co, typically pre-treatment outcomes) |
X1 |
Covariate vector for the treated unit (k x 1) |
Z0 |
Outcome matrix for control units in the pre-treatment window (T_pre x N_co) |
Z1 |
Outcome vector for the treated unit in the pre-treatment window (T_pre x 1) |
max_iter |
Maximum coordinate-descent iterations (default 100) |
tol |
Convergence tolerance on MSPE improvement (default 1e-4) |
t_train |
Validation-window split for V selection. -1 (default): V selected on the full Z window (in-sample). Positive: rows t_train..(T_pre-1) of Z form the validation window used to select V (W is fitted on the full X throughout); after selecting V*, W is refit and the reported loss uses the full Z window. |
z_rows |
Optional 1-based row indices of Z defining the evaluation
window for the outer V optimisation (the |
multistart |
If |
cheap_face |
If |
wolfe |
If |
When t_train > 0, V is selected by minimising MSPE on a validation
window (rows t_train..T_pre-1 of Z) while W is fitted on the full
predictor matrix X. This is appropriate when X is a fixed predictor
matrix that contains no validation-period outcome information (the
user-supplied predictors case). For the outcomes-only case the proper
Abadie (2021) S.3.2 train/validation split is implemented in R
(.scm_oos_outcomes()): candidate W(V) are fitted on training-half
outcomes only, by passing the training rows as X and the validation
rows as Z to this function with t_train = -1.
A list with:
W: Donor weight vector (N_co x 1) on the unit simplex
V: Optimal metric diagonal (k x 1, normalised to sum to 1)
loss: Final pre-treatment prediction loss (full pre-treatment window)
Given unit weights omega and time weights lambda, computes the SDID estimator as a weighted two-way difference:
sdid_estimate_cpp(Y_pre_co, Y_post_co, Y_pre_tr, Y_post_tr, omega, lambda)sdid_estimate_cpp(Y_pre_co, Y_post_co, Y_pre_tr, Y_post_tr, omega, lambda)
Y_pre_co |
Control pre-treatment outcomes (T_pre x N_co) |
Y_post_co |
Control post-treatment outcomes (T_post x N_co) |
Y_pre_tr |
Treated pre-treatment outcomes (T_pre x 1) |
Y_post_tr |
Treated post-treatment outcomes (T_post x 1) |
omega |
Unit weights (N_co x 1) |
lambda |
Time weights (T_pre x 1) |
tau_sdid = (Y_tr_post_mean - Y_tr_pre_wt) - (Y_co_post_wt - Y_co_pre_wt)
A single numeric value: the SDID treatment-effect estimate
tau_sdid.
Computes period-by-period treatment effects (event-study series) for a sharp
Synthetic Difference-in-Differences fit, aligned on the
-weighted pre-treatment gap so that the post-treatment mean
of the period effects equals the SDID ATT (fit$estimate) exactly, together
with bootstrap standard errors and confidence intervals (both normal
approximation and percentile) for every period and for the overall ATT from
the same bootstrap draws.
sdid_event_study( fit, cluster = NULL, data = NULL, n_boot = 200L, level = 0.95, ci_type = c("normal", "percentile"), seed = NULL )sdid_event_study( fit, cluster = NULL, data = NULL, n_boot = 200L, level = 0.95, ci_type = c("normal", "percentile"), seed = NULL )
fit |
A sharp |
cluster |
Cluster specification for block/cluster bootstrap resampling (e.g. prefecture or state):
Unnamed positional vectors are rejected to prevent silent misalignment between row/column order and unit IDs. |
data |
Optional |
n_boot |
Integer number of bootstrap replications (default |
level |
Confidence level for the intervals (default |
ci_type |
Which confidence interval type ( |
seed |
Optional integer seed for reproducible bootstrap resampling. |
Resampling scheme: As in sdid_inference() with method = "bootstrap"
(Arkhangelsky et al. 2021, Algorithm 2), sdid_event_study() resamples
both treated and control units with replacement – here at the cluster
level when cluster is supplied, or at the unit level when
cluster = NULL – and reports period-by-period effects as well as the
ATT. Duplicated units in a
bootstrap draw enter as distinct unit columns (equivalent to assigning fresh
unit IDs per draw), and , (unless user-fixed in
scm_fit()), unit weights , and time weights are
all re-estimated on each draw. Draws in which either the resampled treated
group or the resampled control group has zero units (or fewer than two
control units) are discarded; the counts of valid and discarded draws are
reported in n_boot_valid and n_boot_dropped.
Pre-treatment period coefficients are returned and plotted for visual
inspection of pre-trends only; no joint Wald test is performed because the
pre-treatment gaps are constrained by the and
weights and their covariance matrix can be singular.
An object of class "sdid_event_study" (a named list) containing:
terms: A data.frame with one row per period (time, event_time,
period ("pre" / "post"), estimate, se, ci_norm_lower,
ci_norm_upper, ci_pct_lower, ci_pct_upper, conf.low, conf.high).
mean(terms$estimate[terms$period == "post"]) equals estimate (ATT).
estimate: Overall SDID ATT (fit$estimate).
se: Bootstrap standard error of the ATT (sd(boot_att)).
p_value: Two-sided normal-approximation p-value for the ATT.
ci_norm_lower, ci_norm_upper: Normal-approximation CI for the ATT.
ci_pct_lower, ci_pct_upper: Percentile CI for the ATT.
ci_lower, ci_upper: ATT CI bounds matching ci_type.
ci_type: Resolved default CI type ("normal" or "percentile").
level: Confidence level.
n_boot: Requested number of bootstrap draws.
n_boot_valid: Number of valid bootstrap draws completed.
n_boot_dropped: Number of draws discarded (empty treated/control group
or solver failure).
clustered: Logical; TRUE when cluster was supplied.
n_clusters: Number of unique clusters (or total units when unclustered).
n_treated: Number of treated units in fit.
n_controls: Number of control units in fit.
times: Time index vector (fit$times).
T_pre: Number of pre-treatment periods (fit$T_pre).
boot_gaps: Numeric matrix of dimension n_boot_valid x T holding the
aligned period-by-period gap series for each valid bootstrap draw
(colnames are as.character(times)), or NULL when n_boot == 0.
boot_att: Numeric vector of length n_boot_valid holding the ATT for
each valid draw (identical to rowMeans(boot_gaps[, (T_pre + 1):T])), or
NULL when n_boot == 0.
scm_fit(), sdid_inference(), plot.sdid_event_study()
set.seed(1) panel <- expand.grid(unit = 1:12, year = 2001:2010) panel$pref <- paste0("pref_", (panel$unit - 1) %% 4 + 1) panel$d <- as.integer(panel$unit <= 3 & panel$year >= 2007) panel$y <- panel$unit + 0.5 * (panel$year - 2000) + rnorm(nrow(panel), sd = 0.3) - 1.5 * panel$d fit <- scm_fit(y ~ d | unit + year, data = panel, method = "sdid") es <- sdid_event_study(fit, cluster = "pref", data = panel, n_boot = 50, seed = 42) es broom::tidy(es) plot(es)set.seed(1) panel <- expand.grid(unit = 1:12, year = 2001:2010) panel$pref <- paste0("pref_", (panel$unit - 1) %% 4 + 1) panel$d <- as.integer(panel$unit <= 3 & panel$year >= 2007) panel$y <- panel$unit + 0.5 * (panel$year - 2000) + rnorm(nrow(panel), sd = 0.3) - 1.5 * panel$d fit <- scm_fit(y ~ d | unit + year, data = panel, method = "sdid") es <- sdid_event_study(fit, cluster = "pref", data = panel, n_boot = 50, seed = 42) es broom::tidy(es) plot(es)
Computes standard errors and p-values for a SDID estimate using one of the procedures of Arkhangelsky et al. (2021) and Clarke et al. (2024): the unit bootstrap (Algorithm 2), the jackknife (Algorithm 3), or the placebo permutation (Algorithm 4), plus, for staggered fits, a global jackknife across every unit of all cohorts.
sdid_inference( fit, method = c("placebo", "bootstrap", "jackknife", "jackknife_global"), n_boot = 200L, level = 0.95, alternative = c("two.sided", "greater", "less"), seed = NULL, ... )sdid_inference( fit, method = c("placebo", "bootstrap", "jackknife", "jackknife_global"), n_boot = 200L, level = 0.95, alternative = c("two.sided", "greater", "less"), seed = NULL, ... )
fit |
A |
method |
Inference method: |
n_boot |
Number of bootstrap replications (only for |
level |
Confidence level for the interval (all methods). |
alternative |
Direction of the alternative hypothesis: |
seed |
Integer seed for reproducibility (only for |
... |
Not used. Passing |
method = "bootstrap" draws N units with replacement from the whole
panel, treated and control alike, redrawing a sample with no treated unit
or fewer than two controls, and runs the whole estimator again on it: the
covariate partial-out when the fit has covariates, ,
(unless fixed in scm_fit()), and
. A unit drawn twice enters as two units. On a staggered
fit the cohorts are rebuilt from the treated units drawn, so their sizes
vary across draws.
method = "jackknife" drops one unit at a time, treated or control,
and recomputes the weighted difference-in-differences with the fitted
(renormalised over the remaining donors) and
held fixed; on a staggered fit it does so cohort by cohort
and combines the cohort variances by the delta method, while
"jackknife_global" drops each unit from the whole panel at once and keeps
the covariance between cohorts.
Both resample the treated units, so they need more than one: with a single
treated unit its own noise never varies across replicates (the jackknife is
not even defined), and the function stops. The jackknife on a staggered fit
further needs two treated units in every cohort. Use "placebo" for a
single treated unit; it learns the noise level from the controls and so
assumes their noise matches the treated units'.
For method = "placebo", the p-value is the permutation p-value, while
the standard error is the dispersion of the placebo distribution
(Clarke et al. 2024, Algorithm 4) and the confidence interval is the
normal approximation around the estimate with that SE. The placebo SE
assumes the treated unit's noise is comparable to the control units';
interpret it with caution when the donor pool is small.
sdid_event_study() runs the same bootstrap at the level of clusters of
units and reports period-by-period effects alongside the ATT.
A list with:
estimate: The SDID point estimate.
se: Standard error (placebo: placebo-distribution SD, Algorithm 4).
p_value: Permutation or normal-approximation p-value.
ci_lower, ci_upper: Confidence interval bounds.
method: The inference method used.
n_controls: Number of control units.
alternative: The alternative hypothesis direction.
placebo_effects: Named vector of LOO placebo effects (placebo only).
boot_ests: Bootstrap estimate distribution (bootstrap only).
sdid_event_study() for period-by-period event-study effects and
cluster bootstrap resampling.
For each control unit, treats it as the "pseudo-treated" unit and estimates the leave-one-out SDID effect. The distribution of these placebo effects provides a permutation-based null distribution for inference.
sdid_placebo_cpp(Y_pre, Y_post, time_weights, zeta2)sdid_placebo_cpp(Y_pre, Y_post, time_weights, zeta2)
Y_pre |
Control units pre-treatment outcomes (T_pre x N_co) |
Y_post |
Control units post-treatment outcomes (T_post x N_co) |
time_weights |
Lambda weights for pre-treatment periods (T_pre x 1) |
zeta2 |
Ridge penalty (same as used in the main estimate) |
A numeric vector of length N_co. Each element is the
leave-one-out placebo SDID effect obtained by treating that control unit
as the pseudo-treated unit; the vector serves as a permutation-based null
distribution for inference.
Solves the time-weight QP (with implicit intercept lambda_0 concentrated out):
sdid_time_weights_cpp(Y_pre_co, Y_post_target, zeta_t)sdid_time_weights_cpp(Y_pre_co, Y_post_target, zeta_t)
Y_pre_co |
Pre-treatment outcomes for control units, row-demeaned (T_pre x N_co) |
Y_post_target |
Post-treatment mean per control unit, demeaned (N_co x 1) |
zeta_t |
Ridge penalty for time weights (paper: 1e-6 * sigma_hat) |
min over lambda in Delta_pre: ||Y_post_target - Y_pre_co^T lambda||^2 + zeta_t^2 * N_co * ||lambda||^2
The caller is responsible for pre-demeaning Y_pre_co (row-wise) and Y_post_target (subtract the cross-unit mean) to concentrate out lambda_0, as described in Arkhangelsky et al. (2021) Algorithm 1, Eq. (2.3).
A numeric vector of length T_pre holding the SDID time weights
lambda (non-negative and summing to one).
Solves the regularized QP: min over omega in Delta: sum_t (sum_i omega_i Y_it - Y_tr_t)^2 + zeta^2 * T_pre * ||omega||^2
sdid_unit_weights_cpp(Y_pre, Y_tr_pre, zeta2)sdid_unit_weights_cpp(Y_pre, Y_tr_pre, zeta2)
Y_pre |
Pre-treatment outcome matrix for control units (T_pre x N_co) |
Y_tr_pre |
Pre-treatment outcome vector for treated unit (T_pre x 1), averaged if multiple |
zeta2 |
Ridge penalty parameter (zeta^2). The code internally multiplies by T_pre per the paper. |
This corresponds to equation (5) in Arkhangelsky et al. (2021).
A numeric vector of length N_co holding the SDID unit weights
omega (non-negative and summing to one).
Estimates SE and confidence intervals for the ATT by re-estimating the
fit on resampled units: a bootstrap of the treated and the control units,
drawn separately, or a jackknife that drops one unit at a time. Works for
sharp, staggered and multi-arm SI fits. Each refit keeps the rank k the
fit used.
si_inference( fit, method = c("bootstrap", "jackknife", "jackknife_global"), n_boot = 499L, level = 0.95, alternative = c("two.sided", "greater", "less"), seed = NULL )si_inference( fit, method = c("bootstrap", "jackknife", "jackknife_global"), n_boot = 499L, level = 0.95, alternative = c("two.sided", "greater", "less"), seed = NULL )
fit |
A |
method |
|
n_boot |
Number of bootstrap replications (default 499L; ignored for jackknife). |
level |
Confidence level (default 0.95). |
alternative |
|
seed |
RNG seed for reproducibility (default NULL). |
Holding the treated units fixed would leave out their own noise, which is
most of the variance when few units are treated, so they are resampled
too; with a single treated unit (in any arm) there is nothing to resample
and the function stops. For a sharp fit with one treated unit,
conformal_inference() gives a valid interval.
The controls are drawn once per replicate and shared by every arm and
every cohort, as the SVD basis is a property of the shared donor pool.
Treated units are drawn within their arm. On a staggered fit they are
pooled across cohorts, so cohort sizes vary across replicates, and the
cohorts are rebuilt from the draw and combined with their
N_treated * T_post weights. "jackknife" drops each cohort's own donors
and treated units and combines the per-cohort variances by the delta
method, which ignores the covariance that shared donors induce between
cohorts; "jackknife_global" (staggered fits only) drops every unit from
the whole panel in turn and keeps it. Both jackknives drop treated units,
so every cohort needs at least two (Clarke et al. 2024, Algorithm 3).
A list of class coresynth_inference.
Implements the SI-PCR estimator of Agarwal et al. (2025). Uses the top-k SVD of pre-treatment control outcomes to find donor weights that predict each treated unit's pre-treatment trajectory, then applies those weights to post-treatment control outcomes.
si_pcr_cpp(Y_pre_co, Y_post_co, Y_pre_tr, k)si_pcr_cpp(Y_pre_co, Y_post_co, Y_pre_tr, k)
Y_pre_co |
Pre-treatment control outcomes (T_pre x N_co) |
Y_post_co |
Post-treatment control outcomes (T_post x N_co) |
Y_pre_tr |
Pre-treatment treated outcomes (T_pre x N_tr) |
k |
Number of SVD components to retain |
A list with:
W: Donor weight matrix (N_co x N_tr)
Y_hat: Counterfactual post-treatment outcomes (T_post x N_tr)
Solves: min_L (1/2) ||O o (Y - L)||F^2 + lambda * ||L||* via iterative SVD soft-thresholding (Mazumder, Hastie, Tibshirani 2010). Note: lambda is NOT normalized by |O|.
soft_impute_cpp(Y, O, lambda, max_iter = 1000L, tol = 1e-05)soft_impute_cpp(Y, O, lambda, max_iter = 1000L, tol = 1e-05)
Y |
Observed outcome matrix (T x N). Unobserved entries should be 0. |
O |
Binary mask matrix (T x N): 1 = observed, 0 = missing (treated post). |
lambda |
Nuclear norm penalty (soft-threshold on singular values). |
max_iter |
Maximum iterations. |
tol |
Convergence tolerance (relative Frobenius norm change). |
This is soft_impute_fe_cpp() with the fixed effects switched off.
A numeric matrix of the same dimension as Y (T x N): the
completed low-rank matrix L that minimises the soft-thresholded
nuclear-norm objective.
Solves
min over L, a, b of (1/2) * squaredFrobenius(O o (Y - L - a1' - 1b')) + lambda * nuclearNorm(L)
by alternating a back-fitting update of the fixed effects with one
soft-thresholded SVD step for L (Mazumder, Hastie & Tibshirani 2010;
Athey et al. 2021, S.4 Eq. 4.3).
soft_impute_fe_cpp( Y, O, lambda, estimate_fe = TRUE, max_iter = 1000L, tol = 1e-05, fe_sweeps = 5L )soft_impute_fe_cpp( Y, O, lambda, estimate_fe = TRUE, max_iter = 1000L, tol = 1e-05, fe_sweeps = 5L )
Y |
Observed outcome matrix (T x N). Unobserved entries should be 0. |
O |
Binary mask matrix (T x N): 1 = observed, 0 = missing (treated post). |
lambda |
Nuclear norm penalty (soft-threshold on singular values). |
estimate_fe |
Estimate unpenalised row and column effects. |
max_iter |
Maximum iterations. |
tol |
Convergence tolerance (relative Frobenius norm change). |
fe_sweeps |
Back-fitting sweeps per iteration. The mask makes the two effects non-orthogonal, so one closed-form pass does not solve for both. |
The row and column effects are the unpenalised time and unit fixed effects
of Athey et al. (2021, S.4 Eq. 4.3; see also S.8.1): "we regularize the
estimates of L*, but do not wish to regularize the estimates of the fixed
effects". Leaving them inside L shrinks the panel's level towards zero
along with its factor structure, which biases the imputed counterfactual
whenever the level is far from zero. estimate_fe = false reproduces the
pure low-rank baseline Y = L + e (Mazumder et al. 2010, Eq. 2).
Both effects absorb an additive constant, so a and b are identified
only up to a shift between them; their sum, which is what enters the
completed matrix, is unaffected.
A list with L (the penalised low-rank part), a (row effects),
b (column effects), and M, the completed matrix L + a 1' + 1 b'.
Tensor Unfolding (Matricization) for Synthetic Interventions
tensor_unfold_cpp(T_cube, mode)tensor_unfold_cpp(T_cube, mode)
T_cube |
A 3D array (cube) of dimensions (n1, n2, n3) |
mode |
The mode to unfold along (1, 2, or 3) |
A numeric matrix: the mode-mode unfolding (matricization) of
T_cube, with dimensions n1 x (n2 * n3), n2 x (n1 * n3), or
n3 x (n1 * n2) for mode 1, 2, or 3 respectively.
Coerces a coresynth_inference or sdid_inference object to a one-row
tidy data.frame with broom-style column names so it can be combined
with regression output for paper tables.
## S3 method for class 'coresynth_inference' tidy(x, conf.int = TRUE, ...)## S3 method for class 'coresynth_inference' tidy(x, conf.int = TRUE, ...)
x |
A |
conf.int |
Logical. Include |
... |
Unused. |
A one-row data.frame with columns term, estimate, std.error,
statistic, p.value, conf.low, conf.high, method, alternative,
n_controls, staggered.
tidy() returns the period-by-period event-study table (one row per time
period, with the -aligned point estimate, bootstrap
standard error, and both normal-approximation and percentile confidence
intervals). glance() returns a one-row summary for the overall ATT
alongside the valid and discarded bootstrap replication counts.
## S3 method for class 'sdid_event_study' tidy(x, ci_type = NULL, conf.int = TRUE, ...) ## S3 method for class 'sdid_event_study' glance(x, ci_type = NULL, ...)## S3 method for class 'sdid_event_study' tidy(x, ci_type = NULL, conf.int = TRUE, ...) ## S3 method for class 'sdid_event_study' glance(x, ci_type = NULL, ...)
x |
An |
ci_type |
Optional override ( |
conf.int |
Logical; include |
... |
Unused. |
For tidy(), a data.frame with one row per period and columns
time, event_time, period, estimate, std.error, conf.low,
conf.high, ci_norm_lower, ci_norm_upper, ci_pct_lower,
ci_pct_upper. For glance(), a one-row data.frame with columns
estimate, std.error, p.value, conf.low, conf.high,
ci_norm_lower, ci_norm_upper, ci_pct_lower, ci_pct_upper,
ci_type, level, clustered, n_clusters, n_treated, n_controls,
T_pre, T_post, n_boot, n_boot_valid, n_boot_dropped.
Accessor generics that return the outcome series stored in a fitted
coresynth object under a uniform interface, regardless of the estimation
method:
treated_outcomes(x, ...) ## S3 method for class 'coresynth' treated_outcomes(x, na.rm = FALSE, ...) synthetic_outcomes(x, ...) ## S3 method for class 'coresynth' synthetic_outcomes(x, na.rm = FALSE, ...) ## S3 method for class 'coresynth_staggered' synthetic_outcomes(x, na.rm = FALSE, ...) ## S3 method for class 'coresynth_tasc' synthetic_outcomes(x, na.rm = FALSE, ...) donor_outcomes(x, ...) ## S3 method for class 'coresynth_scm' donor_outcomes(x, ...) ## S3 method for class 'coresynth_sdid' donor_outcomes(x, ...) ## S3 method for class 'coresynth_si' donor_outcomes(x, ...) ## S3 method for class 'coresynth_gsc' donor_outcomes(x, ...) ## S3 method for class 'coresynth_mc' donor_outcomes(x, ...) ## S3 method for class 'coresynth_tasc' donor_outcomes(x, ...) ## S3 method for class 'coresynth' donor_outcomes(x, ...)treated_outcomes(x, ...) ## S3 method for class 'coresynth' treated_outcomes(x, na.rm = FALSE, ...) synthetic_outcomes(x, ...) ## S3 method for class 'coresynth' synthetic_outcomes(x, na.rm = FALSE, ...) ## S3 method for class 'coresynth_staggered' synthetic_outcomes(x, na.rm = FALSE, ...) ## S3 method for class 'coresynth_tasc' synthetic_outcomes(x, na.rm = FALSE, ...) donor_outcomes(x, ...) ## S3 method for class 'coresynth_scm' donor_outcomes(x, ...) ## S3 method for class 'coresynth_sdid' donor_outcomes(x, ...) ## S3 method for class 'coresynth_si' donor_outcomes(x, ...) ## S3 method for class 'coresynth_gsc' donor_outcomes(x, ...) ## S3 method for class 'coresynth_mc' donor_outcomes(x, ...) ## S3 method for class 'coresynth_tasc' donor_outcomes(x, ...) ## S3 method for class 'coresynth' donor_outcomes(x, ...)
x |
A |
... |
Passed to methods. |
na.rm |
Logical; passed to the per-period averaging over multiple
treated units (default |
treated_outcomes(): the treated unit's observed outcome series
(length ). When several units are treated it is their per-period
average – weighted by weights where a sharp fit was given them, and
unweighted across cohorts for a staggered fit, whose treated units do not
share a treatment date.
synthetic_outcomes(): the estimated counterfactual series
(length ), i.e. the synthetic control or model-fitted outcome.
donor_outcomes(): the matrix of observed donor
(control unit) outcomes over all periods.
Each accessor returns NULL when the requested series is not stored in
the fit. Staggered-adoption fits keep their counterfactual per cohort (in
fit$cohort_fits), because collapsing it to one series would average
cohorts whose treatment starts at different times, so
synthetic_outcomes() returns NULL for them; use augment() for the
per-cohort series.
For treated_outcomes() and synthetic_outcomes(), a numeric
vector of length , or NULL. For donor_outcomes(), a
numeric matrix (donors in columns, named when unit
names are available), or NULL.
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 = "scm") y1 <- treated_outcomes(fit) # observed treated series y1_0 <- synthetic_outcomes(fit) # synthetic counterfactual Yco <- donor_outcomes(fit) # donor outcome matrix all.equal(y1 - y1_0, unname(fit$gap))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 = "scm") y1 <- treated_outcomes(fit) # observed treated series y1_0 <- synthetic_outcomes(fit) # synthetic counterfactual Yco <- donor_outcomes(fit) # donor outcome matrix all.equal(y1 - y1_0, unname(fit$gap))