Getting started with bacipair

Why this package exists

bacipair is the audit-grade methodology stack behind a Carbon Market Register (CER) Australian Carbon Credit Unit (ACCU) Human-Induced Regeneration (HIR) monitoring, reporting and verification (MRV) audit. Before any real register data is touched, the synthetic-control extension must first recover the oracle average treatment effect on the treated (ATT) on simulated staggered-treatment ground truth, at or below a pre-registered bias floor of 0.15 percentage points (pp).

This vignette runs that ground-truth-recovery check on a small simulated panel. The audit-grade version of the same check (12 replicates at n_pairs = 200) achieved a mean bias of 0.0347 pp against the floor of 0.15 pp for the estimator this package ships (Section 3 below); the run shown here is a single small-n replicate and its own numbers, printed inline, will differ.

1. Estimand and estimators

The target is the panel-average ATT, the mean effect of treatment on the treated units in the periods after they are treated:

\[ \mathrm{ATT} \;=\; \mathbb{E}\bigl[Y_{it}(1) - Y_{it}(0) \mid D_i = 1,\; t \ge T_i\bigr], \]

where \(Y_{it}(1)\) and \(Y_{it}(0)\) are the treated and untreated potential outcomes for unit \(i\) at time \(t\), \(D_i\) indicates treatment, and \(T_i\) is unit \(i\)’s treatment year. Because \(Y_{it}(0)\) is never observed for a treated unit after \(T_i\), every estimator below is a different way of constructing a counterfactual for it.

Paired difference-in-differences (baci_paired_did()). On the within-pair-and-arm-demeaned design, the fitted regression is the standard two-way fixed-effects specification

\[ Y_{ipt} \;=\; \alpha_{ip} \;+\; \gamma_t \;+\; \tau \, D_{ip}\!\times\!\mathrm{Post}_{pt} \;+\; \varepsilon_{ipt}, \]

with a pair-and-arm fixed effect \(\alpha_{ip}\), a year fixed effect \(\gamma_t\), and \(\tau\), the treated-arm-by-post-treatment interaction coefficient, as the ATT. Both fixed-effect blocks are absorbed by within-(pair, arm) demeaning before the interaction term is estimated by ordinary least squares, so att_hat is \(\hat{\tau}\) from that demeaned regression. Standard errors are the CR1 cluster-robust sandwich (Cameron and Miller, 2015):

\[ \widehat{V} \;=\; c \cdot (X'X)^{-1}\left(\sum_{g=1}^{G} X_g' u_g u_g' X_g\right)(X'X)^{-1}, \qquad c \;=\; \frac{G}{G-1}\cdot\frac{N-1}{N-K}, \]

clustered by proponent (default) or by pair, with \(G\) clusters, \(N\) observations and \(K\) estimated coefficients.

Synthetic control (synth_control_baci()). For each treated unit, a donor weight vector \(w\) over the pool of other units’ comparison arms is chosen to minimise squared distance to the treated unit’s pre-treatment trajectory and static predictors, subject to a simplex constraint:

\[ \hat{w} \;=\; \operatorname*{arg\,min}_{w \,\in\, \Delta^{k-1}} \; \lVert X w - y \rVert_2^2, \qquad \Delta^{k-1} \;=\; \Bigl\{ w \in \mathbb{R}^k : w_j \ge 0,\; \textstyle\sum_j w_j = 1 \Bigr\}, \]

where \(X\) stacks the donors’ pre-treatment outcome trajectory and standardised static predictors (baseline cover, elevation, rainfall climatology) and \(y\) is the treated unit’s corresponding vector. The package solves this by projected gradient descent onto the simplex, using the exact Euclidean projection of Wang and Carreira-Perpinan (2013). The per-period ATT is then the observed outcome minus the fitted counterfactual, \(\widehat{Y_{it}(0)} = X_t \hat{w}\).

2. Simulate a staggered-treatment BACI panel

The simulator carries an explicit oracle: the panel-average ATT is exactly mean(ite) over treated post-treatment observations. Every distribution and seed is recorded.

library(bacipair)

sim <- simulate_staggered_panel(
  seed    = 20260520L,
  n_pairs = 80L
)

cat("Oracle panel-average ATT:",
    round(sim$oracle_panel_att, 3), "pp\n")
#> Oracle panel-average ATT: 0.797 pp
cat("Treatment-year distribution:\n")
#> Treatment-year distribution:
print(table(sim$pair_table$treatment_year))
#> 
#> 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022 
#>   12   11    3    5    5    8    8    6    9   13
cat("Proponent cluster sizes (top 10):\n")
#> Proponent cluster sizes (top 10):
print(head(sort(table(sim$pair_table$proponent_id),
                decreasing = TRUE), 10L))
#> 
#>  1  2  3  6  4  7 37 11 14 15 
#> 24 13  6  4  3  3  3  2  2  2

3. Naive paired-DiD

The matched 1:1 comparator with TWFE is the natural starting point and the foil for the synthetic-control estimator. Because elevation and rainfall-climatology covary with treatment timing in this data-generating process, the paired-DiD coefficient carries a measurable residual bias.

did <- baci_paired_did(sim$panel)

cat(sprintf("Paired-DiD ATT_hat = %.3f pp  (cluster-robust SE = %.3f)\n",
            did$att_hat, did$se_cluster))
#> Paired-DiD ATT_hat = 0.682 pp  (cluster-robust SE = 0.163)
cat(sprintf("95%% CI: [%.3f, %.3f]\n", did$ci_lo, did$ci_hi))
#> 95% CI: [0.361, 1.002]
cat(sprintf("Oracle: %.3f pp  ->  bias = %+0.3f pp\n",
            sim$oracle_panel_att,
            did$att_hat - sim$oracle_panel_att))
#> Oracle: 0.797 pp  ->  bias = -0.115 pp

Averaged over 12 seeds at n_pairs = 120, this estimator’s bias is -0.256 pp against a 0.70 pp oracle – 37 per cent of the oracle effect, and roughly twice its own cluster-robust standard error (0.138 pp). That bias is why Section 4’s synthetic control exists as the primary estimator, with paired-DiD kept only as the naive foil.

4. Time-varying synthetic-control comparator

The donor pool is every other pair’s comparison arm; per-treated-unit weights are simplex-constrained and chosen to minimise the joint objective of Section 1 over pre-treatment outcomes and standardised static predictors (baseline cover, elevation, rainfall climatology).

sc <- synth_control_baci(sim$panel)

cat(sprintf("SC panel ATT_hat = %.3f pp\n", sc$panel_att_hat))
#> SC panel ATT_hat = 0.961 pp
cat(sprintf("SC pre-fit RMSE  = %.3f pp\n", sc$pre_fit_rmse))
#> SC pre-fit RMSE  = 0.870 pp
cat(sprintf("SC bias vs oracle = %+0.3f pp\n",
            sc$panel_att_hat - sim$oracle_panel_att))
#> SC bias vs oracle = +0.165 pp

At package defaults the synthetic-control estimator nets out most of the elevation and climatology coupling that biases paired-DiD in Section 3: on the audit-grade run (12 replicates, n_pairs = 200) its mean bias is 0.0347 pp, against the pre-registered floor of 0.15 pp. The estimator shown in this vignette is the uncorrected variant that ships in the package – it has no covariate-regression bias correction and no pre-trend de-trending step, and it is that uncorrected variant’s numbers quoted above, not a corrected extension. On this vignette’s small-n single-replicate run, expect a wider envelope than the audit-grade figure.

resid_rank <- order(vapply(sc$fits, `[[`, numeric(1L), "pre_residual"))
best_k  <- resid_rank[1L]
worst_k <- resid_rank[length(resid_rank)]

op <- par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
for (k in c(best_k, worst_k)) {
  fk <- sc$fits[[k]]
  pd <- fk$panel[order(get("year"))]
  plot(pd$year, pd$woody_cover, type = "b", pch = 16,
       xlab = "Year", ylab = "Woody cover (pp)",
       main = sprintf("Pair %s (pre-fit RMSE = %.2f)",
                       fk$pair_id, fk$pre_residual),
       cex.main = 0.85,
       ylim = range(c(pd$woody_cover, pd$counterfactual)))
  lines(pd$year, pd$counterfactual, lty = 2, col = "firebrick")
  abline(v = fk$treatment_year, lty = 3, col = "grey40")
  legend("topleft", legend = c("Treated", "Synthetic control"),
         lty = c(1, 2), pch = c(16, NA), col = c("black", "firebrick"),
         bty = "n", cex = 0.75)
}

Pre-treatment fit and post-treatment counterfactual gap for the best- and worst-fitting treated units, ranked by pre-treatment RMSE. The synthetic control (dashed) tracks the treated unit (solid) closely before the vertical treatment-year line in both panels; after it, the gap between the two lines is that unit’s estimated treatment effect.

par(op)

The best-fitting unit’s synthetic control tracks its pre-treatment trajectory closely and opens a clear, stable post-treatment gap; the worst-fitting unit’s pre-fit RMSE is an order of magnitude larger, so its post-treatment gap is a correspondingly noisier estimate of that unit’s treatment effect – the pre-fit RMSE is a per-unit diagnostic for how much to trust the post-treatment number next to it.

pu <- sc$per_unit_att
plot_range <- range(c(pu$oracle_unit, pu$att_hat_unit))
plot(pu$oracle_unit, pu$att_hat_unit, pch = 16, col = "steelblue",
     xlab = "Oracle unit ATT (pp)", ylab = "Estimated unit ATT (pp)",
     xlim = plot_range, ylim = plot_range)
abline(0, 1, lty = 2, col = "grey40")

Per-unit synthetic-control ATT estimate against the simulator’s oracle unit-level treatment effect. Points on the dashed 1:1 line are exactly recovered; the vertical distance of a point from the line is that unit’s estimation error.

Most units sit close to the 1:1 line; the panel-average bias reported above is the mean vertical offset of this cloud, and it is small next to the unit-to-unit scatter – the estimator is closer to unbiased in aggregate than it is precise for any one treated unit, which is the usual synthetic-control trade-off between panel-level and unit-level accuracy.

donor_ids <- sort(unique(sim$panel[get("arm") == "comparison",
                                    get("pair_id")]))
treated_ids <- vapply(sc$fits, `[[`, numeric(1L), "pair_id")
w_mat <- matrix(0, nrow = length(sc$fits), ncol = length(donor_ids),
                dimnames = list(as.character(treated_ids),
                                as.character(donor_ids)))
for (k in seq_along(sc$fits)) {
  fk <- sc$fits[[k]]
  used_donors <- donor_ids[donor_ids != fk$pair_id]
  w_mat[k, as.character(used_donors)] <- fk$weights
}
image(t(w_mat), axes = FALSE,
      xlab = "Donor pair (comparison arm)", ylab = "Treated unit",
      col = hcl.colors(20, "YlOrRd", rev = TRUE))

Donor-weight matrix: each row is a treated unit and each column a comparison-arm donor pair, shaded by the simplex weight that unit’s synthetic control places on that donor. Darker cells mark heavier weight.

Weight is concentrated on a small subset of donors per treated unit rather than spread evenly across the whole pool – the simplex constraint lets the optimiser find a sparse, targeted comparator per unit instead of averaging over every donor, which is what makes each unit’s counterfactual identifiable from the pre-treatment fit in the first figure.

5. Proponent-level cluster bootstrap

One proponent runs many adjacent projects (Macintosh et al., 2024), so the right inferential unit is the proponent cluster. The bootstrap resamples whole proponents with replacement and refits the chosen estimator.

boot <- proponent_cluster_bootstrap(
  sim$panel,
  estimator = "paired_did",
  n_boot    = 60L,
  seed      = 20260521L
)

cat(sprintf("Bootstrap (paired-DiD): SE = %.3f pp,  95%% CI [%.3f, %.3f]\n",
            boot$se, boot$ci_lo, boot$ci_hi))
#> Bootstrap (paired-DiD): SE = 0.156 pp,  95% CI [0.493, 1.059]

6. Goodman-Bacon decomposition (diagnostic)

The TWFE coefficient is a weighted average of all 2x2 DiD comparisons implied by the treatment-timing structure. The Goodman-Bacon (2021) decomposition exposes those weights; large weight on Treated vs Untreated and small, non-negative weight on Earlier vs Later Treated is the “clean TWFE” signature.

if (have_bacondecomp) {
  gb <- goodman_bacon_decomp(sim$panel)
  print(gb$weight_summary)
  cat(sprintf("Implied TWFE estimate from 2x2 weights: %.3f pp\n",
              gb$twfe_estimate))
} else {
  message("bacondecomp not installed - skipping decomposition.")
}
#>                       type    weight weighted_estimate
#> 1 Earlier vs Later Treated 0.4539027         1.1331389
#> 2 Later vs Earlier Treated 0.5460973        -0.7097637
#> Implied TWFE estimate from 2x2 weights: 0.423 pp

7. Oracle-recovery log

write_oracle_recovery_log() writes the oracle-vs-estimate comparison computed above to a JSON file in the schema bacipair_targets_plan() uses as its terminal target, so a targets pipeline and this vignette report identical numbers from identical inputs.

if (have_jsonlite) {
  log_path <- file.path(tempdir(), "oracle-recovery-vignette.json")
  write_oracle_recovery_log(
    panel             = sim,
    synth_control_fit = sc,
    paired_did_fit    = did,
    cluster_boot      = boot,
    path              = log_path
  )

  log <- jsonlite::fromJSON(log_path)
  str(log[c("oracle_panel_att", "sc_panel_att_hat", "sc_minus_oracle",
            "did_att_hat", "did_minus_oracle",
            "bootstrap_se", "bootstrap_ci_lo", "bootstrap_ci_hi")])
} else {
  message("jsonlite not installed - skipping log demonstration.")
}
#> List of 8
#>  $ oracle_panel_att: num 0.797
#>  $ sc_panel_att_hat: num 0.961
#>  $ sc_minus_oracle : num 0.165
#>  $ did_att_hat     : num 0.682
#>  $ did_minus_oracle: num -0.115
#>  $ bootstrap_se    : num 0.156
#>  $ bootstrap_ci_lo : num 0.493
#>  $ bootstrap_ci_hi : num 1.06

8. Drop-in `targets` skeleton

bacipair_targets_plan() returns a list of tar_target_raw() calls ready to drop into a _targets.R:

if (have_targets) {
  plan <- bacipair_targets_plan(
    simulate_args = list(seed = 20260520L, n_pairs = 80L),
    boot_args     = list(n_boot = 60L, seed = 20260521L)
  )
  print(names(plan))
} else {
  message("targets not installed - skipping plan demonstration.")
}
#> [1] "panel"               "synth_control_fit"   "paired_did_fit"     
#> [4] "cluster_boot"        "goodman_bacon"       "oracle_recovery_log"

Reproducibility

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] bacipair_0.1.0 rmarkdown_2.32
#> 
#> loaded via a namespace (and not attached):
#>  [1] vctrs_0.7.3         cli_3.6.6           knitr_1.52         
#>  [4] rlang_1.3.0         xfun_0.61           otel_0.2.0         
#>  [7] processx_3.9.0      targets_1.12.0      data.table_1.18.6.1
#> [10] jsonlite_2.0.0      glue_1.8.1          prettyunits_1.2.0  
#> [13] bacondecomp_0.1.1   buildtools_1.0.0    backports_1.5.1    
#> [16] htmltools_0.5.9     maketools_1.3.2     sys_3.4.3          
#> [19] ps_1.9.3            sass_0.4.10         tibble_3.3.1       
#> [22] evaluate_1.0.5      jquerylib_0.1.4     base64url_1.4      
#> [25] fastmap_1.2.0       yaml_2.3.12         lifecycle_1.0.5    
#> [28] compiler_4.6.1      igraph_2.3.4        codetools_0.2-20   
#> [31] pkgconfig_2.0.3     digest_0.6.39       R6_2.6.1           
#> [34] tidyselect_1.2.1    pillar_1.11.1       magrittr_2.0.5     
#> [37] callr_3.8.0         bslib_0.12.0        tools_4.6.1        
#> [40] secretbase_1.3.1    cachem_1.1.0