sun_abraham: thread vcov_type ∈ {classical, hc1, hc2, hc2_bm} (Phase 1b 1/8)

Adds `vcov_type` parameter to `SunAbraham`, mirroring the DiD/MPD/TWFE
chain from Phase 1a. Defaults to "hc1" (preserves prior bit-equal
behavior - SA historically hard-coded HC1). First PR of Phase 1b, which
threads `vcov_type` through the 8 standalone estimators that expose
`cluster=` but not yet `vcov_type=`.

Methodology: when `vcov_type ∈ {classical, hc2, hc2_bm}`,
`_fit_saturated_regression` auto-routes to a full-dummy saturated design
(intercept + cohort × event-time interactions + unit dummies + time
dummies). FWL preserves cohort coefficients but not the hat matrix, so
HC2 leverage and Bell-McCaffrey DOF must be computed on the full FE
projection. Mirrors TWFE Gate 1 from PR #469. Empirically matches
`lm() + sandwich::vcovHC(type="HC2")` and
`clubSandwich::vcovCR(..., type="CR2") + coef_test()$df_Satt` at
atol=1e-10 (pinned in tests/test_methodology_sun_abraham.py).

Scope limits: replicate-weight survey + hc2/hc2_bm raises
NotImplementedError (per-replicate full-dummy refit not implemented).
`vcov_type="conley"` rejected at __init__ with a deferral message
(threading conley_* params is a follow-up). Auto-cluster-at-unit is
dropped when the user opts into explicit `vcov_type="hc2"` or
`"classical"` (both one-way only); preserved for `"hc1"` and
`"hc2_bm"`.

Documented deviation from R: SA's within-transform HC1 SE differs from
`fixest::sunab()` by ~1-2% on typical panel sizes (different (n-k)
finite-sample correction). The IW aggregation is otherwise identical;
parity at atol=5e-3.

Test surface: 15 new behavioral tests in test_sun_abraham.py covering
default-vs-explicit bit-equality, all four vcov_type values
finite-and-distinct, auto-cluster drop/preserve, replicate-weight
reject, get_params/set_params, clone+repeat-fit idempotence, invalid
value rejection, cluster_var=None cascade through survey-PSU injection,
full-dummy vs within-transform HC2 divergence. 4 new R-parity tests in
test_methodology_sun_abraham.py against sandwich/clubSandwich/fixest
goldens.

New R golden scenario `sun_abraham_two_cohort` in
benchmarks/data/clubsandwich_cr2_golden.json (5 cohorts × 8 periods
panel; pins classical_se, hc2_se, cr2_bm_singleton_se+dof,
cr2_bm_unit_se+dof, sunab_hc1_event_study_e0_se).

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

Author: igerber

CHANGELOG.md

@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 
 ### Added
+- **SunAbraham `vcov_type` parameter (Phase 1b PR 1/8).** `SunAbraham(vcov_type=...)` now accepts `{"classical","hc1","hc2","hc2_bm"}` (defaults to `"hc1"`, which preserves prior behavior bit-equally - SA historically hard-coded HC1). Auto-cluster-at-unit dropped when the user opts into explicit `vcov_type="hc2"` or `vcov_type="classical"` (one-way only); preserved for `"hc1"` and `"hc2_bm"`. When `vcov_type in {"classical","hc2","hc2_bm"}`, `_fit_saturated_regression` auto-routes to a full-dummy saturated design (mirrors TWFE Gate 1 from PR #469): FWL preserves cohort coefficients but not the hat matrix, so HC2 leverage and Bell-McCaffrey Satterthwaite DOF must be computed on the full FE projection. Empirically matches R `lm()` summary classical SE, `sandwich::vcovHC(type="HC2")`, and `clubSandwich::vcovCR(..., type="CR2")` + `coef_test()$df_Satt` at atol=1e-10 (cohort SE and BM DOF pinned in `tests/test_methodology_sun_abraham.py`). `vcov_type` is now propagated to `SunAbrahamResults.vcov_type` for downstream introspection. `SurveyDesign` (any kind — analytical weights, stratified, PSU, or replicate-weight) combined with `vcov_type in {"classical","hc2","hc2_bm"}` raises `NotImplementedError`: the survey-design TSL (or replicate-weight refit) variance overrides the analytical sandwich family, and the auto-cluster guard for one-way families would silently downgrade unit-level PSUs to per-observation PSUs. Use `vcov_type="hc1"` (default) for survey designs. `conley` rejected at `__init__` with a deferral message (would require threading 6+ `conley_*` params through the saturated regression call). **Deviation from R:** SA's within-transform HC1 SE differs from `fixest::sunab()` by ~1-2% (~2e-3 absolute) on typical panel sizes due to a different `(n-k)` finite-sample correction (fixest counts absorbed FE in k_total; SA's `solve_ols` counts only within-transformed columns); the IW aggregation step is otherwise identical (pinned at atol=5e-3, tracked in TODO.md). First PR of the Phase 1b standalone-estimator threading initiative (7 PRs to follow: StackedDiD, WooldridgeDiD-OLS, CallawaySantAnna, ImputationDiD, TripleDifference, TwoStageDiD, EfficientDiD).
 - **PreTrendsPower R `pretrends` parity goldens (PR-C closes PR-B's deferred R-parity row).** JSON goldens at `benchmarks/data/r_pretrends_golden.json` generated from the committed `benchmarks/R/generate_pretrends_golden.R` script against `jonathandroth/pretrends` commit `122731d082` (package version 0.1.0, R 4.5.2). 4 fixtures cover regular K=3 grid (`uniform_3_pre_periods_no_anticipation`), irregular K=3 grid `[-5,-3,-1]` (`irregular_pre_periods` — locks the PR-B Step 4 γ-unit linear-weight fix), anticipation-shifted K=4 grid (`anticipation_shifted`), and K=1 closed form (`single_pre_period_closed_form` — Roth Proposition 2 univariate truncated-normal). `TestPretrendsParityR` in `tests/test_methodology_pretrends.py` now active (4 tests): NIS power vs R `pretrends::pretrends()` at `atol=1e-4` across all 4 fixtures × 4 γ values; γ_p MDV vs R `slope_for_power()` at `atol=1e-4` across all 4 fixtures × 2 target_power values; end-to-end `fit()` on irregular grid vs R γ_p at `atol=1e-4` (locks the full `fit() → _extract_pre_period_params → _get_violation_weights → _compute_mdv_nis` chain through the public API); K=1 three-way cross-check (Python ≡ analytical truncated-normal closed form `1 - Φ(z - γ/σ) + Φ(-z - γ/σ)` at `atol=1e-7`; both within `atol=1e-4` of R). Tolerance rationale: R hardcodes `thresholdTstat.Pretest=1.96` while Python uses `scipy.stats.norm.ppf(0.975) = 1.959963984540054` (`dz ≈ 3.6e-5`); R `slope_for_power` uses `uniroot(tol = .Machine$double.eps^0.25 ≈ 1.22e-4)` versus Python `brentq(xtol=2e-12)`; the inverse-solver tolerance gap dominates γ_p, and `mvtnorm::pmvnorm` (R) vs `scipy.stats.multivariate_normal.cdf` (Python) Genz-Bretz randomized-lattice differences bound the K=4 NIS power gap at ~5e-5. `METHODOLOGY_REVIEW.md` PreTrendsPower row promoted `**Complete** (R parity pending)` → `**Complete**`. Roth (2022) paper review's `R \`pretrends\` package version pin (provisional)` Gaps bullet struck. Closes the PR-C TODO row.
 - **`SpilloverDiD(survey_design=...)` integration on HC1 / CR1 paths via Binder TSL (Wave E.1).** Lifts the Wave B/C/D upfront `NotImplementedError` and adds design-based variance for `vcov_type ∈ {"hc1"}` plus `cluster=<col>` (CR1). **Documented synthesis** of Gerber (2026, arXiv:2605.04124) Proposition 1 — Binder Taylor Series Linearization for IF representations of smooth functionals; explicitly derived for TwoStageDiD in the paper's Appendix — composed with the Wave D Gardner GMM first-stage uncertainty correction (Butts 2021 §3.1 + Gardner 2022 §4) applied to SpilloverDiD's ring-indicator stage-2 design. No reference software combines all ingredients. **Mechanical composition:** SpilloverDiD's per-obs Wave D IF `psi_i = gamma_hat' * X_{10,i} * eps_{10,i} - X_{2,i} * eps_{2,i}` (with survey weights threaded through `gamma_hat` solve, eps construction, and bread inversion via Hájek normalization) is aggregated to PSU totals and passed to the audited `_compute_stratified_meat_from_psu_scores` Binder TSL meat helper. Stage-1 FE estimation extends `_iterative_fe_subset` with a `weights=` kwarg implementing WLS-FE via weighted bincount (numerator `bincount(w*resid)` / denominator `bincount(w)`); the `weights is None` path is bit-identical to the Wave B / C / D unweighted bincount. **Degrees of freedom:** t-distribution lookup uses `ResolvedSurveyDesign.df_survey` (4-way branch: PSU+strata → `n_PSU - n_strata`; PSU only → `n_PSU - 1`; strata only → `n_obs - n_strata`; neither → `n_obs - 1`), threaded through all four `safe_inference` call sites (aggregate `tau_total`, per-ring `delta_j`, event-study per-event-time `tau_k` / `delta_jk`, scalar `att` lincom). **Survey-array subsetting:** when `finite_mask` drops baseline-treated rows, `survey_weights` and `ResolvedSurveyDesign.{weights, strata, psu, fpc, replicate_weights}` are subsetted in parallel; `n_psu`, `n_strata`, and `survey_metadata` are recomputed (mirrors `TwoStageDiD.fit:567-601`). **Cluster + survey resolution:** when `cluster=<col>` and `survey_design.psu` are both supplied with different groupings, a `UserWarning` fires and PSU wins (mirrors `_resolve_effective_cluster` at `survey.py:1253-1275`; TwoStageDiD parity). When `cluster=<col>` is supplied without `survey_design.psu`, the cluster column is injected as the effective PSU via `_inject_cluster_as_psu`, which now honors `SurveyDesign.nest`: under `nest=False`, cluster labels must be globally unique across strata (raises if they repeat, matching the explicit-PSU resolver's contract). **Saturated `df_survey = 0` NaN-fail:** when `lonely_psu="remove"` removes all strata (singleton PSUs), the meat helper returns `(_, var_computed=False, legit_zero=0)` and SpilloverDiD's Wave E.1 path returns NaN meat with a `UserWarning` matching `"df_survey"` so callers can `pytest.warns(UserWarning, match="df_survey")`. This is a **departure from TwoStageDiD** (`two_stage.py:2003-2005`) which currently NaN-fails SILENTLY; Wave E.1 surfaces the diagnostic per `feedback_no_silent_failures`. **Subpopulation limitation (Wave E.3 follow-up):** `SurveyDesign.subpopulation()`-derived designs with zero-weight padding rows that lose stage-1 FE support have those rows physically removed by `finite_mask`, so `n_psu` / `df_survey` / Binder centering reflect the reduced fit sample rather than the full domain design (documented in REGISTRY; Wave E.3 will preserve full-design bookkeeping). **Public surface restrictions:** `vcov_type="conley" + survey_design=` raises `NotImplementedError` pointing at planned Wave E.2 (Conley × survey product-kernel synthesis with within-stratum Conley sandwich on PSU totals); replicate-weight variance (BRR / Fay / JK1 / JKn / SDR) raises `NotImplementedError` — per Gerber (2026) Appendix A, the IF-reweighting shortcut does not apply to TwoStageDiD-class estimators because `gamma_hat` is weight-sensitive; correct support requires per-replicate full re-fit and is queued as a follow-up; non-pweight (`weight_type ∈ {"fweight", "aweight"}`) raises `ValueError` (the Binder TSL assumes probability weights). **Implementation:** `_compute_gmm_corrected_meat` extended with `survey_weights` + `resolved_survey` kwargs at `diff_diff/two_stage.py:56` (TYPE_CHECKING forward reference for `ResolvedSurveyDesign` to avoid circular import); new module-level helper `_compute_binder_tsl_meat` at `diff_diff/two_stage.py` wraps `_compute_stratified_meat_from_psu_scores` with implicit per-obs PSU synthesis for no-PSU survey designs + the Wave E.1 NaN-fail + warning; `_iterative_fe_subset` weighted path at `diff_diff/spillover.py:1382` (in-place extension, bit-identical fallback, positive-weight identification gate); `_inject_cluster_as_psu` honors `nest` (shared survey-helper fix that also benefits TwoStageDiD); `ResolvedSurveyDesign` gains a `nest` field propagated through all 5 construction sites. `SpilloverDiDResults` extended with `survey_metadata`, `n_psu`, `n_strata` fields at `diff_diff/results.py`. **Tests:** new `TestSpilloverDiDWaveE1SurveyDesignHc1` (17 tests: bit-identity fallback, Binder TSL hand-check uniform + non-uniform weights, lonely_psu modes, FPC degenerate limits ×3, saturated NaN-fail with `pytest.warns(match="df_survey")`, cluster+survey warn-and-use-PSU, no-PSU regressions (weights-only, weights+strata, cluster-without-PSU, cluster overlap with nest=False/True), zero-weight Omega_0 exclusion + all-zero raises, replicate-weight + non-pweight + Conley+survey rejections, fit idempotency, finite_mask subsetting) and `TestSpilloverDiDWaveE1SurveyDesignEventStudy` (7 tests: event-study + survey on both `is_staggered` branches with `df_survey` lincom verification, distinguishability between survey-share and sample-share lincom rules via manual reconstruction with cohort-correlated weights + non-constant tau_k, aggregate-vs-event-study parity, drift goldens, subset-path invariant). Wave B/C/D bullets below are unchanged; this entry replaces the pre-Wave-E.1 `survey_design=` rejection.
 

TODO.md

@@ -98,7 +98,9 @@ Deferred items from PR reviews that were not addressed before merge.
 | PreTrendsPower: CS/SA `anticipation=1` R-parity fixture. The PR-C R-parity goldens cover NIS power + γ_p MDV at `atol=1e-4` on four shifted-grid / regular / irregular / K=1 fixtures, but R `pretrends` has no anticipation parameter so the Python-side `_extract_pre_period_params` anticipation filter (`if t < _pre_cutoff` in `pretrends.py` lines 1138-1150 for CS; mirror in SA branch) is not R-parity-locked. Build a synthetic `CallawaySantAnnaResults` (or `SunAbrahamResults`) with `anticipation=1` and a t=-1 event-study entry that should be filtered before reaching `_compute_power_nis`, then assert the resulting γ_p matches R's `slope_for_power()` on the K=4 shifted-grid fixture. Existing PR-B MC-based tests (`TestPretrendsPropositions`) and full-VCV tests (`TestPretrendsCovarianceSource`) already cover the filter mechanically; this would close the loop against R. | `tests/test_methodology_pretrends.py::TestPretrendsParityR`, `benchmarks/R/generate_pretrends_golden.R` | PR-C follow-up | Low |
 
 
-| Thread `vcov_type` (classical / hc1 / hc2 / hc2_bm) through the 8 standalone estimators that expose `cluster=`: `CallawaySantAnna`, `SunAbraham`, `ImputationDiD`, `TwoStageDiD`, `TripleDifference`, `StackedDiD`, `WooldridgeDiD`, `EfficientDiD`. Phase 1a added `vcov_type` to the `DifferenceInDifferences` inheritance chain only. | multiple | Phase 1a | Medium |
+| Thread `vcov_type` (classical / hc1 / hc2 / hc2_bm) through the 7 standalone estimators that expose `cluster=` but not yet `vcov_type=`: `CallawaySantAnna`, `ImputationDiD`, `TwoStageDiD`, `TripleDifference`, `StackedDiD`, `WooldridgeDiD`, `EfficientDiD`. Phase 1a added the chain to DiD/MPD/TWFE; Phase 1b PR 1/8 added `SunAbraham` (this row tracks the remaining 7). | multiple | Phase 1b | Medium |
+| Extend `SunAbraham` with `vcov_type="conley"` (Conley spatial-HAC) as a first-class feature: thread `conley_coords` / `conley_cutoff_km` / `conley_metric` / `conley_kernel` / `conley_time` / `conley_unit` / `conley_lag_cutoff` through `_fit_saturated_regression`. Phase 1b PR 1/8 deferred this; SA currently rejects `vcov_type="conley"` at `__init__` with a deferral message. | `diff_diff/sun_abraham.py` | follow-up | Medium |
+| Harmonize SunAbraham's HC1 within-transform finite-sample correction with `fixest::sunab()`. SA's `solve_ols` applies `n / (n - k_dm)` (within-transform columns only); fixest applies `n / (n - k_total)` (counts absorbed FE). SE values differ by ~1-2% on typical panel sizes (documented in REGISTRY.md "Deviation from R"; pinned at `atol=5e-3` in `tests/test_methodology_sun_abraham.py`). Either thread `df_adjustment` into the vcov scaling or document as an intentional difference. | `diff_diff/sun_abraham.py`, `diff_diff/linalg.py::compute_robust_vcov` | follow-up | Low |
 | Weighted one-way Bell-McCaffrey (`vcov_type="hc2_bm"` + `weights`, no cluster) currently raises `NotImplementedError`. `_compute_bm_dof_from_contrasts` builds its hat matrix from the unscaled design via `X (X'WX)^{-1} X' W`, but `solve_ols` solves the WLS problem by transforming to `X* = sqrt(w) X`, so the correct symmetric idempotent residual-maker is `M* = I - sqrt(W) X (X'WX)^{-1} X' sqrt(W)`. Rederive the Satterthwaite `(tr G)^2 / tr(G^2)` ratio on the transformed design and add weighted parity tests before lifting the guard. | `linalg.py::_compute_bm_dof_from_contrasts`, `linalg.py::_validate_vcov_args` | Phase 1a | Medium |
 | Weighted CR2 Bell-McCaffrey cluster-robust (`vcov_type="hc2_bm"` + `cluster_ids` + `weights`) currently raises `NotImplementedError`. Weighted hat matrix and residual rebalancing need threading per clubSandwich WLS handling. | `linalg.py::_compute_cr2_bm` | Phase 1a | Medium |
 | `TwoWayFixedEffects(vcov_type in {"hc2","hc2_bm"})` with replicate-weight survey designs raises `NotImplementedError` (`twfe.py:~233`). The replicate path re-demeans per replicate (re-demeaning depends on the per-replicate weight vector), which doesn't compose with the full-dummy HC2/HC2-BM build — a correct implementation would need per-replicate full-dummy refit. Workaround: use `vcov_type="hc1"` for replicate-weight CR1. | `twfe.py::fit` | follow-up | Low |

benchmarks/R/generate_clubsandwich_golden.R

@@ -283,6 +283,123 @@ output$twfe_two_period <- list(
   dof_bm_unit = as.numeric(ct_twfe_unit$df_Satt)
 )
 
+# --- SunAbraham saturated regression HC2 / HC2-BM scenario (Phase 1b PR 1/8) -
+# Mirrors SunAbraham(vcov_type in {"classical","hc2","hc2_bm"}) on a
+# 5-cohort × 8-period balanced panel. SA's Part G auto-route builds a
+# full-dummy saturated design when vcov_type needs the hat matrix —
+# matches lm(y ~ D_ge interactions + factor(unit) + factor(period)).
+# Targets the (g=4, e=0) cohort × event-time interaction (the canonical
+# at-treatment effect of the earliest cohort).
+
+set.seed(42)
+n_units_per_cohort <- 8
+n_sa_periods <- 8
+sa_cohorts <- c(0, 4, 5, 6, 7)  # 0 = never-treated
+
+d_sa <- expand.grid(u_in_cohort = seq_len(n_units_per_cohort),
+                    period = seq_len(n_sa_periods),
+                    cohort_idx = seq_along(sa_cohorts))
+d_sa <- d_sa[order(d_sa$cohort_idx, d_sa$u_in_cohort, d_sa$period), ]
+d_sa$unit <- (d_sa$cohort_idx - 1) * n_units_per_cohort + d_sa$u_in_cohort - 1
+d_sa$first_treat <- sa_cohorts[d_sa$cohort_idx]
+d_sa$time <- d_sa$period
+d_sa$rel_time <- ifelse(d_sa$first_treat > 0,
+                         d_sa$time - d_sa$first_treat, -999L)
+sa_unit_fe <- rnorm(max(d_sa$unit) + 1, mean = 0, sd = 1)
+d_sa$treated <- as.integer(d_sa$first_treat > 0 & d_sa$time >= d_sa$first_treat)
+d_sa$y <- sa_unit_fe[d_sa$unit + 1] + 0.3 * d_sa$time +
+          1.0 * d_sa$treated + rnorm(nrow(d_sa), sd = 0.5)
+
+# Build cohort × event-time interaction columns (excluding ref period -1).
+# Sanitize negative event times for R formula compatibility (e=-3 → "n3").
+sa_treatment_groups <- sort(unique(d_sa$first_treat[d_sa$first_treat > 0]))
+sa_all_rel_times <- sort(unique(d_sa$rel_time[d_sa$first_treat > 0]))
+sa_all_rel_times <- sa_all_rel_times[sa_all_rel_times != -1]
+sa_interaction_cols <- c()
+sa_col_map <- list()
+for (g in sa_treatment_groups) {
+  for (e in sa_all_rel_times) {
+    e_safe <- if (e < 0) paste0("n", abs(e)) else as.character(e)
+    col_name <- paste0("D_", g, "_", e_safe)
+    original_name <- paste0("D_", g, "_", e)
+    ind <- as.integer(d_sa$first_treat == g & d_sa$rel_time == e)
+    if (sum(ind) > 0) {
+      d_sa[[col_name]] <- ind
+      sa_interaction_cols <- c(sa_interaction_cols, col_name)
+      sa_col_map[[original_name]] <- col_name
+    }
+  }
+}
+
+sa_target_orig <- "D_4_0"  # the (g=4, e=0) interaction
+sa_target_safe <- sa_col_map[[sa_target_orig]]
+stopifnot(!is.null(sa_target_safe))
+
+sa_rhs <- paste(c(sa_interaction_cols, "factor(unit)", "factor(time)"),
+                collapse = " + ")
+fit_sa <- lm(as.formula(paste("y ~", sa_rhs)), data = d_sa)
+sa_coef_names <- names(coef(fit_sa))
+sa_target_idx <- which(sa_coef_names == sa_target_safe)
+stopifnot(length(sa_target_idx) == 1L)
+
+# Extract SE/DOF for the target only (atol=1e-10 pin in Python tests).
+sa_classical_se <- summary(fit_sa)$coefficients[sa_target_safe, "Std. Error"]
+sa_vcov_hc2 <- sandwich::vcovHC(fit_sa, type = "HC2")
+sa_hc2_se <- sqrt(sa_vcov_hc2[sa_target_safe, sa_target_safe])
+# Singleton-cluster CR2 reduces to one-way HC2-BM.
+sa_vcov_cr2_singleton <- vcovCR(fit_sa, cluster = seq_len(nrow(d_sa)),
+                                type = "CR2")
+sa_cr2_singleton_se <- sqrt(sa_vcov_cr2_singleton[sa_target_safe,
+                                                   sa_target_safe])
+sa_ct_singleton <- coef_test(fit_sa, vcov = sa_vcov_cr2_singleton)
+sa_dof_bm_singleton <- sa_ct_singleton[sa_target_safe, "df_Satt"]
+# CR2-BM clustered at unit (the SA auto-cluster default for hc2_bm).
+sa_vcov_cr2_unit <- vcovCR(fit_sa, cluster = d_sa$unit, type = "CR2")
+sa_cr2_unit_se <- sqrt(sa_vcov_cr2_unit[sa_target_safe, sa_target_safe])
+sa_ct_unit <- coef_test(fit_sa, vcov = sa_vcov_cr2_unit)
+sa_dof_bm_unit <- sa_ct_unit[sa_target_safe, "df_Satt"]
+# fixest::sunab() parity for SA's HC1 cluster-at-unit default path.
+# SA HC1 uses within-transform; fixest also uses within-transform.
+# Note: fixest::sunab requires a specific encoding — first_treat=0 means
+# never-treated. fixest auto-handles that.
+suppressPackageStartupMessages(library(fixest, quietly = TRUE))
+fit_sunab <- fixest::feols(
+  y ~ sunab(first_treat, time) | unit + time,
+  data = d_sa,
+  cluster = ~unit
+)
+# fixest::sunab aggregates to event-study coefficients (IW-aggregated
+# across cohorts). The coefficient labels are "time::<event_time>".
+# Compare SA's event_study_effects[0] (overall e=0 ATT) against fixest's
+# "time::0" event-study SE.
+sunab_coef_table <- as.data.frame(summary(fit_sunab)$coeftable)
+sunab_target_label <- "time::0"
+sunab_hc1_es0_se <- if (sunab_target_label %in% rownames(sunab_coef_table)) {
+  sunab_coef_table[sunab_target_label, "Std. Error"]
+} else {
+  warning("Could not locate fixest sunab event-study target ",
+          sunab_target_label)
+  NA_real_
+}
+
+output$sun_abraham_two_cohort <- list(
+  unit = d_sa$unit,
+  time = d_sa$time,
+  first_treat = d_sa$first_treat,
+  y = d_sa$y,
+  target_cohort_g = 4L,
+  target_event_time_e = 0L,
+  target_col_safe = sa_target_safe,
+  classical_se = unname(sa_classical_se),
+  hc2_se = unname(sa_hc2_se),
+  cr2_bm_singleton_se = unname(sa_cr2_singleton_se),
+  dof_bm_singleton = unname(sa_dof_bm_singleton),
+  cr2_bm_unit_se = unname(sa_cr2_unit_se),
+  dof_bm_unit = unname(sa_dof_bm_unit),
+  sunab_hc1_event_study_e0_se = unname(sunab_hc1_es0_se),
+  sunab_event_study_target_label = sunab_target_label
+)
+
 output$meta <- list(
   source = "clubSandwich",
   clubSandwich_version = as.character(packageVersion("clubSandwich")),