From 3f856ccaa6a554adee106861cb527faee5aa4cbd Mon Sep 17 00:00:00 2001 From: Yiqing Xu Date: Mon, 8 Jun 2026 21:31:40 -0700 Subject: [PATCH 01/14] feat(conformal): estimator-agnostic core (scores + rank p-value + interval inversion) Phase 1 (part 1) of vartype='conformal'. R/conformal.R: .conformal_score (studentized default; ratio = Abadie post/pre RMSPE; meanabs; rmse), .conformal_pval (weighted-capable), .conformal_ci_meanabs (closed-form att.avg), .conformal_ci_grid (refit-free grid inversion for the other scores), .conformal_neff. Isolated check: meanabs coverage 0.902 at target 0.90 (Nco=39); grid studentized CI brackets truth; resolution guard unbounds when alpha < 1/(Nco+1). Next: leave-one-control-out calibration via impute_Y0 (staggered-aware via valid_controls) + vartype dispatch in boot.R/default.R; not yet wired. Spec: statsclaw-workspace/fect/runs/REQ-conformal/spec.md Co-Authored-By: Claude Opus 4.8 --- R/conformal.R | 110 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 110 insertions(+) create mode 100644 R/conformal.R diff --git a/R/conformal.R b/R/conformal.R new file mode 100644 index 00000000..4717d217 --- /dev/null +++ b/R/conformal.R @@ -0,0 +1,110 @@ +## ----------------------------------------------------------------------------- +## conformal.R +## +## Cross-sectional conformal inference for separated, symmetric counterfactual +## estimators (vartype = "conformal"). The treated counterfactual prediction +## error is calibrated against the leave-one-control-out prediction errors of the +## donors and read off as a RANK, rather than resampled and normal-wrapped +## (parametric bootstrap) or quantile-resampled (nonparametric bootstrap). +## +## Design: statsclaw-workspace/fect/runs/REQ-conformal/spec.md +## Decisions (2026-06-08): jackknife+ default (>= 1 - 2 alpha); scores +## "studentized" (default), "ratio" (Abadie post/pre RMSPE), "meanabs", "rmse"; +## staggered adoption via valid_controls(predictive = "notyettreated"). +## +## This file holds the estimator-AGNOSTIC core: scores, the conformal p-value, +## and interval inversion. The calibration loop (leave-one-control-out via +## impute_Y0) and the vartype dispatch live in boot.R / default.R and feed this +## core a control-by-period residual matrix plus the treated gap path. +## ----------------------------------------------------------------------------- + +## ---- 1. nonconformity score ------------------------------------------------- +## e.post : numeric, a unit's post-treatment residual (prediction-error) path +## e.pre : numeric, the same unit's pre-treatment residual path (for normaliz.) +## type : "studentized" | "ratio" | "meanabs" | "rmse" +## Returns a scalar. Validity (exchangeability) is identical across types; the +## choice affects interval width / adaptivity only. +.conformal_score <- function(e.post, e.pre, type = "studentized") { + e.post <- e.post[is.finite(e.post)] + if (length(e.post) == 0L) return(NA_real_) + rms <- function(x) sqrt(mean(x^2)) + if (type == "meanabs") { + return(abs(mean(e.post))) + } else if (type == "rmse") { + return(rms(e.post)) + } else if (type == "ratio") { + ## Abadie post/pre RMSPE ratio. + ep <- e.pre[is.finite(e.pre)] + denom <- if (length(ep) > 0L) rms(ep) else NA_real_ + if (is.na(denom) || denom <= 0) return(NA_real_) + return(rms(e.post) / denom) + } else if (type == "studentized") { + ep <- e.pre[is.finite(e.pre)] + denom <- if (length(ep) > 1L) stats::sd(ep) else NA_real_ + if (is.na(denom) || denom <= 0) return(NA_real_) + return(mean(abs(e.post)) / denom) + } + stop("conformal: unknown score type '", type, "'.") +} + +## ---- 2. conformal p-value (weighted-capable) -------------------------------- +## s.tr : scalar treated score; s.co : numeric vector of control scores. +## w.co : optional non-negative control weights (weighted conformal); NULL = 1. +## w.tr : treated self-weight (default 1). Returns p in (0, 1]. +.conformal_pval <- function(s.tr, s.co, w.co = NULL, w.tr = 1) { + ok <- is.finite(s.co) + s.co <- s.co[ok] + if (is.null(w.co)) w.co <- rep(1, length(s.co)) else w.co <- w.co[ok] + num <- w.tr + sum(w.co[s.co >= s.tr]) + den <- w.tr + sum(w.co) + num / den +} + +## effective sample size for weighted conformal +.conformal_neff <- function(w) { + w <- w[is.finite(w) & w > 0] + if (length(w) == 0L) return(0) + (sum(w)^2) / sum(w^2) +} + +## ---- 3. interval inversion -------------------------------------------------- +## Closed form for the average effect under the "meanabs" score. +## g.tr : scalar treated mean post-period gap; g.co : control mean gaps. +## Returns c(lower, upper) for att.avg. Width = +/- the rank-corrected +## (1 - alpha) quantile of |g.co|. +.conformal_ci_meanabs <- function(g.tr, g.co, alpha = 0.05, w.co = NULL) { + ag <- abs(g.co[is.finite(g.co)]) + n <- length(ag) + k <- ceiling((1 - alpha) * (n + 1)) # conformal quantile rank + if (k > n) { # alpha < 1/(n+1): unbounded + return(c(-Inf, Inf)) + } + q <- sort(ag)[k] + c(g.tr - q, g.tr + q) +} + +## Grid inversion for ratio / studentized / rmse (refit-free: control scores are +## tau-independent; only the treated score is recomputed on the grid). +## gap.tr.post / gap.tr.pre : treated gap paths. s.co : control scores. +## Returns c(lower, upper) = range of accepted constant effects tau. +.conformal_ci_grid <- function(gap.tr.post, gap.tr.pre, s.co, type, alpha = 0.05, + w.co = NULL, w.tr = 1, n.grid = 401L, span = NULL) { + center <- mean(gap.tr.post[is.finite(gap.tr.post)]) + if (is.null(span)) { + sd.co <- stats::sd(s.co[is.finite(s.co)]) + span <- max(abs(gap.tr.post), na.rm = TRUE) + 6 * (if (is.finite(sd.co)) sd.co else 1) + } + grid <- seq(center - span, center + span, length.out = n.grid) + accept <- vapply(grid, function(tau) { + s.tr <- .conformal_score(gap.tr.post - tau, gap.tr.pre, type) + if (!is.finite(s.tr)) return(FALSE) + .conformal_pval(s.tr, s.co, w.co = w.co, w.tr = w.tr) > alpha + }, logical(1)) + if (!any(accept)) return(c(NA_real_, NA_real_)) + rng <- range(grid[accept]) + ## flag if the acceptance set hit a grid edge (interval may be unbounded) + if (accept[1] || accept[n.grid]) { + attr(rng, "edge") <- TRUE + } + rng +} From 08b84c381f027bb8de8e5b52bc9907a26be40e93 Mon Sep 17 00:00:00 2001 From: Yiqing Xu Date: Mon, 8 Jun 2026 22:09:44 -0700 Subject: [PATCH 02/14] feat(conformal): leave-one-control-out calibration + vartype dispatch (WIP) Phase 1 part 2 of vartype='conformal'. R/conformal.R: conformal_calibrate() does the deterministic leave-one-control-out calibration, reusing valid_controls() + impute_Y0() (the held-out, non-resampled analog of draw.error). It runs end to end on fect_boot's internal matrices (N_calib = full donor pool); the fix for the impute_Y0 effect-aggregation error was to pass the real hasRevs rather than 1. R/default.R accepts vartype='conformal' (+ mc guard). R/boot.R adds a conformal branch after the point fit that computes the calibration and currently stops with a diagnostic. tests/testthat/test-conformal.R covers the estimator-agnostic core (scores, meanabs coverage ~0.90, p-value, resolution guard, n_eff) and passes. Not done (next): wire the result through fect's output slots (eff.calendar/est.att/est.avg) by mirroring the parametric slot assembly (an early return skips the calendar-effect computation fect.default expects); fix the studentized grid CI; validate full-pipeline coverage (the LOO-gap assembly currently yields frequent unbounded intervals to debug; over-coverage is plausibly jackknife+ conservatism); full testthat + coverage-study gates. Existing vartypes are unaffected (conformal branch is gated and early-returns; validation change is additive); bootstrap regression smoke passes. Spec: statsclaw-workspace/fect/runs/REQ-conformal/spec.md Co-Authored-By: Claude Opus 4.8 --- R/boot.R | 28 +++++++ R/conformal.R | 133 ++++++++++++++++++++++++++++++++ R/default.R | 8 +- tests/testthat/test-conformal.R | 51 ++++++++++++ 4 files changed, 216 insertions(+), 4 deletions(-) create mode 100644 tests/testthat/test-conformal.R diff --git a/R/boot.R b/R/boot.R index 45149d85..28a468a4 100644 --- a/R/boot.R +++ b/R/boot.R @@ -543,6 +543,34 @@ fect_boot <- function( fit.out <- out$Y.ct N_unit <- dim(out$res)[2] + ## ---- vartype = "conformal": cross-sectional conformal interval ----------- + ## Rank the treated gap against the leave-one-control-out donor scores + ## (see conformal.R), bypassing the resampling / jackknife SE machinery. + ## Uses fect_boot's internal preprocessed matrices (Y, D, I, II, T.on). + ## NOTE Phase 1: score / alpha hardcoded; est.att per-period bands + arg + ## threading from fect() come next. + if (vartype == "conformal") { + cc <- conformal_calibrate( + Y = Y, D = D, X = X, I = I, II = II, T.on = T.on, + r.cv = out$r.cv, eff = out$eff, + method = method, predictive = time.component.from, + force = force, hasRevs = hasRevs, tol = tol, max.iteration = max.iteration, + norm.para = norm.para, + score = "studentized", alpha = 0.05 + ) + ## WIP (Phase 1 part 2): the calibration runs on the internal matrices and + ## yields the att.avg interval below. Wiring the result through fect()'s + ## output slots (eff.calendar, est.att, est.avg) must mirror the parametric + ## path's slot assembly (boot.R est.* construction), not an early return, + ## which skips the calendar-effect computation that fect.default expects. + stop(sprintf(paste0( + "vartype = 'conformal': calibration is implemented but output integration ", + "is in progress. Result: att = %.4f, 95%% CI = [%.4f, %.4f], p = %.4f, ", + "N_calib = %d, score = %s, form = %s."), + cc$att, cc$ci[1], cc$ci[2], cc$p.value, cc$n.calib, cc$score, cc$form), + call. = FALSE) + } + if (!is.null(group)) { group.output.origin <- out$group.output group.output.name <- names(out$group.output) diff --git a/R/conformal.R b/R/conformal.R index 4717d217..57729560 100644 --- a/R/conformal.R +++ b/R/conformal.R @@ -108,3 +108,136 @@ } rng } + +## ---- 4. calibration: deterministic leave-one-control-out -------------------- +## Mirror of boot.R::draw.error but HELD-OUT (not resampled): each valid control +## is predicted once from the other controls, giving its out-of-fold residual +## path. Estimator-agnostic — built-ins via impute_Y0(method); a custom separated +## learner via `conformal.fit`. Inputs are the preprocessed matrices available +## inside fect_boot(): Y, D, X, I, II, T.on (all TT x N; X is TT x N x p or NULL), +## r.cv (selected rank), and eff (TT x N point-fit gap matrix). +## +## Returns: list(att, ci, p.value, score.tr, score.co, n.calib, n.eff, score, form). +## Scope (verified): block design (common onset), method gsynth/ife, nevertreated +## (separated) calibration. Staggered uses the union post-window (per-cohort +## windows are a TODO); pooled (notyettreated) fits fall back to nevertreated with +## a warning. MC is not yet routed through impute_Y0 (upstream stop()). +conformal_calibrate <- function(Y, D, X = NULL, I, II, T.on, r.cv, eff, + method = "gsynth", predictive = "nevertreated", + force = 3L, hasRevs = 0L, tol = 1e-5, + max.iteration = 1000L, norm.para = NULL, + score = "studentized", alpha = 0.05, + conformal.full = FALSE, conformal.weight = "none", + conformal.fit = NULL) { + + TT <- nrow(Y); N <- ncol(Y) + sum.D <- colSums(D) + id.tr <- which(sum.D > 0); id.co <- which(sum.D == 0) + Ntr <- length(id.tr); Nco <- length(id.co) + if (Nco < 2L) stop("conformal: need at least 2 controls.") + + ## --- separation guard: conformal calibration is controls-only (nevertreated) + if (!identical(predictive, "nevertreated")) { + warning("vartype = \"conformal\" requires a separated fit; ", + "calibrating on controls only (predictive = \"nevertreated\").") + predictive <- "nevertreated" + } + + ## --- valid controls (enough pre/post obs); reuse fect's screen + valid.co <- valid_controls(list(D = D, I = I, r.cv = r.cv), method, predictive, force) + valid.co <- intersect(valid.co, id.co) + if (length(valid.co) < 2L) { + stop("conformal: fewer than 2 valid controls after screening.") + } + + ## --- post / pre period windows (block: common onset; staggered: union) + post.idx <- which(rowSums(D[, id.tr, drop = FALSE]) > 0) # any treated on + pre.idx <- setdiff(seq_len(TT), post.idx) + if (length(post.idx) == 0L || length(pre.idx) == 0L) { + stop("conformal: could not identify pre/post windows from D.") + } + if (Ntr > 1L) { + onsets <- vapply(id.tr, function(i) min(which(D[, i] == 1)), integer(1)) + if (length(unique(onsets)) > 1L) { + warning("conformal: staggered onsets detected; using the union post-window. ", + "Per-cohort windows are not yet implemented.") + } + } + + ## --- one held-out fit per valid control, reusing the parametric refit path + d.pattern <- D[, id.tr[1]] # a treated D column (defines onset) + sub3 <- function(A, idx) if (is.null(A)) NULL else A[, idx, , drop = FALSE] + + ## fake-treated column = held-out control j's DATA carrying a treated unit's + ## TIMING (D, T.on, II must agree with the assigned onset; only Y/I are j's). + loo_gap <- function(j) { + co.rest <- setdiff(valid.co, j) + if (!is.null(conformal.fit)) { + y0 <- conformal.fit(Y = Y, X = X, time = seq_len(TT), + control.ids = co.rest, target.id = j, T0 = max(pre.idx)) + return(Y[, j] - y0) + } + Y.ps <- cbind(Y[, j], Y[, co.rest, drop = FALSE]) + D.ps <- cbind(d.pattern, D[, co.rest, drop = FALSE]) + Ton.ps <- cbind(T.on[, id.tr[1]], T.on[, co.rest, drop = FALSE]) + I.ps <- cbind(I[, j], I[, co.rest, drop = FALSE]) + II.ps <- cbind(I[, j] * (d.pattern == 0), II[, co.rest, drop = FALSE]) + X.ps <- if (is.null(X)) NULL else + array(c(X[, j, , drop = FALSE], X[, co.rest, , drop = FALSE]), + dim = c(TT, 1L + length(co.rest), dim(X)[3])) + synth <- try(impute_Y0( + method = method, predictive = predictive, + Y = Y.ps, X = X.ps, D = D.ps, W = NULL, I = I.ps, II = II.ps, + T.on = Ton.ps, tuning = r.cv, boot = 1, + force = force, hasRevs = hasRevs, tol = tol, + max.iteration = max.iteration, norm.para = norm.para + ), silent = TRUE) + if (inherits(synth, "try-error") || !("eff" %in% names(synth))) return(rep(NA_real_, TT)) + g <- as.matrix(synth$eff.tr)[, 1] + if (!is.null(norm.para)) g <- g / norm.para[1] + g + } + + ## --- control scores (tau-independent) and mean post gaps + obs <- (I == 1) + S.co <- rep(NA_real_, length(valid.co)) + g.co <- rep(NA_real_, length(valid.co)) + for (m in seq_along(valid.co)) { + j <- valid.co[m] + gj <- loo_gap(j) + ep <- gj[post.idx][obs[post.idx, j]] + eq <- gj[pre.idx ][obs[pre.idx, j]] + S.co[m] <- .conformal_score(ep, eq, score) + g.co[m] <- mean(ep, na.rm = TRUE) + } + keep <- is.finite(S.co) + S.co <- S.co[keep]; g.co <- g.co[keep] + Ncal <- length(S.co) + + ## --- treated statistic (averaged over treated units, post window) + eff.tr.mat <- eff[, id.tr, drop = FALSE] + tr.post.path <- rowMeans(eff.tr.mat[post.idx, , drop = FALSE], na.rm = TRUE) + tr.pre.path <- rowMeans(eff.tr.mat[pre.idx, , drop = FALSE], na.rm = TRUE) + g.tr <- mean(tr.post.path, na.rm = TRUE) + + ## --- weights (overlap) — placeholder until loading_bound wiring (Phase 3) + w.co <- NULL + if (identical(conformal.weight, "overlap")) { + warning("conformal.weight = \"overlap\" not yet wired; using unweighted.") + } + + ## --- interval + if (score == "meanabs") { + ci <- .conformal_ci_meanabs(g.tr, g.co, alpha = alpha, w.co = w.co) + } else { + ci <- .conformal_ci_grid(tr.post.path, tr.pre.path, S.co, score, + alpha = alpha, w.co = w.co) + } + S.tr <- .conformal_score(tr.post.path, tr.pre.path, score) + p.value <- .conformal_pval(S.tr, S.co, w.co = w.co) + + list(att = g.tr, ci = ci, p.value = p.value, + score.tr = S.tr, score.co = S.co, + n.calib = Ncal, n.eff = if (is.null(w.co)) Ncal else .conformal_neff(w.co), + score = score, form = if (conformal.full) "full" else "jackknife+") +} diff --git a/R/default.R b/R/default.R index ed89bae0..5e13c876 100644 --- a/R/default.R +++ b/R/default.R @@ -560,15 +560,15 @@ fect.default <- function( method_arg <- method if (se == 1) { - if (!vartype %in% c("bootstrap", "jackknife", "parametric")) { + if (!vartype %in% c("bootstrap", "jackknife", "parametric", "conformal")) { stop( - "\"vartype\" must be one of \"bootstrap\", \"jackknife\", or \"parametric\".", + "\"vartype\" must be one of \"bootstrap\", \"jackknife\", \"parametric\", or \"conformal\".", call. = FALSE ) } - if (vartype == "parametric" && method %in% c("mc", "both")) { + if (vartype %in% c("parametric", "conformal") && method %in% c("mc", "both")) { stop( - "The \"parametric\" option is not available for the \"mc\" or \"both\" methods." + "The \"", vartype, "\" option is not available for the \"mc\" or \"both\" methods." ) } if (vartype == "jackknife" && !is.null(cl)) { diff --git a/tests/testthat/test-conformal.R b/tests/testthat/test-conformal.R new file mode 100644 index 00000000..12846ff7 --- /dev/null +++ b/tests/testthat/test-conformal.R @@ -0,0 +1,51 @@ +## Tests for the estimator-agnostic conformal core (R/conformal.R). +## The full vartype = "conformal" pipeline (leave-one-control-out calibration via +## impute_Y0) runs end to end but its output-slot integration and coverage +## validation are in progress; those tests are added with that wiring. + +test_that("conformal scores are finite and ordered as expected", { + set.seed(1) + e.post <- rnorm(12); e.pre <- rnorm(24) + for (s in c("meanabs", "rmse", "ratio", "studentized")) { + expect_true(is.finite(fect:::.conformal_score(e.post, e.pre, s))) + } + ## a larger post deviation gives a larger score (meanabs) + expect_gt(fect:::.conformal_score(e.post + 5, e.pre, "meanabs"), + fect:::.conformal_score(e.post, e.pre, "meanabs")) + ## ratio with zero pre-variation is NA, not Inf-by-accident + expect_true(is.na(fect:::.conformal_score(e.post, rep(0, 24), "ratio"))) +}) + +test_that("meanabs closed-form interval covers near nominal", { + set.seed(20260608) + reps <- 3000; alpha <- 0.10; Nco <- 39; att <- 3 + cov <- replicate(reps, { + g.co <- rnorm(Nco) # no-effect donor mean gaps + g.tr <- att + rnorm(1) # treated mean gap = att + exchangeable noise + ci <- fect:::.conformal_ci_meanabs(g.tr, g.co, alpha = alpha) + (att >= ci[1]) && (att <= ci[2]) + }) + expect_gt(mean(cov), 0.86) + expect_lt(mean(cov), 0.94) # ~ 1 - alpha +}) + +test_that("conformal p-value lies in (0,1] and respects weights", { + s.co <- rnorm(20) + p <- fect:::.conformal_pval(0, s.co) + expect_gt(p, 0); expect_lte(p, 1) + ## a very large treated score is the most extreme -> smallest p = 1/(n+1) + expect_equal(fect:::.conformal_pval(100, s.co), 1 / (length(s.co) + 1)) + ## zero weight on a donor drops it from the count + w <- rep(1, 20); w[which.max(s.co)] <- 0 + expect_true(is.finite(fect:::.conformal_pval(0, s.co, w.co = w))) +}) + +test_that("resolution guard: alpha below 1/(Nco+1) yields an unbounded interval", { + ci <- fect:::.conformal_ci_meanabs(0, rnorm(20), alpha = 1 / 25) # 1/25 < 1/21 + expect_true(all(is.infinite(ci))) +}) + +test_that("effective sample size is sensible", { + expect_equal(fect:::.conformal_neff(rep(1, 10)), 10) + expect_lt(fect:::.conformal_neff(c(rep(1, 9), 100)), 10) # one big weight -> n_eff < n +}) From 1d16ea47360a07c1ed1ff5c7178c6e6a11ae5e6d Mon Sep 17 00:00:00 2001 From: Yiqing Xu Date: Mon, 8 Jun 2026 22:45:39 -0700 Subject: [PATCH 03/14] fix(conformal): resolve coverage items 1 & 2 (unbounded + studentized empties) Validates and corrects the leave-one-control-out calibration. 200-rep block DGP (true effect 0, alpha=0.10): total coverage 0.92-0.95, 0 unbounded across all four scores. ## Item 1 -- unbounded intervals (already fixed; now locked) The "~29/50 unbounded" note predated the hasRevs plumbing fix. Verified 0 unbounded over 800 calibrate calls. Mild over-coverage (0.915 vs 0.90) is the treated-vs-control LOO asymmetry (treated gap from the full fit, each control from an Nco-1 fit) plus jackknife+ conservatism -- safe, not a bug. ## Item 2 -- studentized "flaky NAs" were misdiagnosed They are correct conformal empties: with per-unit pre-period studentization, a treated unit that draws a small sd(pre) has a score exceeding every control at every tau, so the data legitimately reject all constant effects. A wide-grid (2001-pt) diagnostic put the best-tau p-value at the 1/(Nco+1) floor -> genuine, not a grid-span miss. Counted as misses, studentized total coverage is nominal. ## Fixes - .conformal_ci_grid: dimensional span bug fixed. Span was max|gap| + 6*sd(score), mixing the outcome-unit tau grid with score units. Now max|gap| + (max score + 1)*denom via new .conformal_denom(e.pre, type) (sd for studentized, rms for ratio, 1 for rmse). Added adaptive edge-expansion and explicit +/-Inf on a genuinely unbounded side. - Honest empty handling: empty acceptance returns c(NA,NA) with attr "empty"="rejected_all_tau"; conformal_calibrate surfaces status in {ok, empty, unbounded} and warns (no more silent NA). - Default score -> "meanabs" in the boot.R conformal branch (was "studentized"). meanabs is the closed-form level interval for the average effect: conservative, never empty/unbounded above the resolution floor, matches fect's headline att.avg. Dispersion scores stay available (calibrated, can be empty). ## Verification - tests/testthat/test-conformal.R: 33 pass (added grid bracket/empty/unbounded/ expansion, .conformal_denom, + 2 full-pipeline coverage guards: meanabs never empties/unbounds & covers ~nominal; studentized total coverage ~nominal with every empty flagged, no silent NA). - bootstrap & jackknife vartypes smoke-tested unaffected (changes are gated inside the vartype=="conformal" branch; conformal.R is all-new functions). Refs statsclaw runs/2026-06-08-conformal-inference.md (Phase 1 part 3). Co-Authored-By: Claude Opus 4.8 --- R/boot.R | 7 +- R/conformal.R | 87 +++++++++++++++++++----- tests/testthat/test-conformal.R | 117 ++++++++++++++++++++++++++++++-- 3 files changed, 188 insertions(+), 23 deletions(-) diff --git a/R/boot.R b/R/boot.R index 28a468a4..9c894fe5 100644 --- a/R/boot.R +++ b/R/boot.R @@ -556,7 +556,12 @@ fect_boot <- function( method = method, predictive = time.component.from, force = force, hasRevs = hasRevs, tol = tol, max.iteration = max.iteration, norm.para = norm.para, - score = "studentized", alpha = 0.05 + ## meanabs = closed-form level interval for the average effect: conservative, + ## never empty/unbounded (for Ncal above the resolution floor), and matches + ## fect's headline att.avg. The dispersion scores (studentized/ratio/rmse) + ## are constant-effect interval inversions -- correctly calibrated but can + ## return an empty set ~alpha of the time -- exposed later via a score arg. + score = "meanabs", alpha = 0.05 ) ## WIP (Phase 1 part 2): the calibration runs on the internal matrices and ## yields the att.avg interval below. Wiring the result through fect()'s diff --git a/R/conformal.R b/R/conformal.R index 57729560..c4f1056f 100644 --- a/R/conformal.R +++ b/R/conformal.R @@ -83,29 +83,66 @@ c(g.tr - q, g.tr + q) } +## Outcome-scale denominator a given score divides by, used to size the grid in +## the SAME units as tau (the grid is over outcome-scale constant effects, so the +## span must be outcome-scaled, not score-scaled). +.conformal_denom <- function(e.pre, type) { + ep <- e.pre[is.finite(e.pre)] + if (type == "ratio") return(if (length(ep) > 0L) sqrt(mean(ep^2)) else NA_real_) + if (type == "studentized") return(if (length(ep) > 1L) stats::sd(ep) else NA_real_) + 1 # rmse / meanabs: numerator already in outcome units +} + ## Grid inversion for ratio / studentized / rmse (refit-free: control scores are ## tau-independent; only the treated score is recomputed on the grid). ## gap.tr.post / gap.tr.pre : treated gap paths. s.co : control scores. -## Returns c(lower, upper) = range of accepted constant effects tau. +## Returns c(lower, upper) = range of accepted constant effects tau. Two special +## returns, both carried as attributes so the caller can react: +## attr "empty" = "rejected_all_tau" : acceptance set is genuinely empty (no +## constant effect is consistent at level alpha; common with per-unit +## studentization when the treated pre-period sd is small) -> c(NA, NA). +## attr "edge" = TRUE : acceptance reached the (expanded) grid +## edge; the interval is effectively unbounded on that side. .conformal_ci_grid <- function(gap.tr.post, gap.tr.pre, s.co, type, alpha = 0.05, w.co = NULL, w.tr = 1, n.grid = 401L, span = NULL) { - center <- mean(gap.tr.post[is.finite(gap.tr.post)]) - if (is.null(span)) { - sd.co <- stats::sd(s.co[is.finite(s.co)]) - span <- max(abs(gap.tr.post), na.rm = TRUE) + 6 * (if (is.finite(sd.co)) sd.co else 1) + gp <- gap.tr.post[is.finite(gap.tr.post)] + center <- mean(gp) + ## span in OUTCOME units: beyond |tau - center| ~ (max control score) * denom + ## the treated score must exceed every control score, so acceptance is + ## impossible -- this bounds where the grid needs to look. + denom <- .conformal_denom(gap.tr.pre, type) + if (!is.finite(denom) || denom <= 0) denom <- 1 + s.hi <- suppressWarnings(max(s.co[is.finite(s.co)])) + if (!is.finite(s.hi)) s.hi <- 1 + if (is.null(span)) span <- max(abs(gp)) + (s.hi + 1) * denom + + accept_on <- function(span) { + grid <- seq(center - span, center + span, length.out = n.grid) + keep <- vapply(grid, function(tau) { + s.tr <- .conformal_score(gap.tr.post - tau, gap.tr.pre, type) + if (!is.finite(s.tr)) return(FALSE) + .conformal_pval(s.tr, s.co, w.co = w.co, w.tr = w.tr) > alpha + }, logical(1)) + list(grid = grid, keep = keep) } - grid <- seq(center - span, center + span, length.out = n.grid) - accept <- vapply(grid, function(tau) { - s.tr <- .conformal_score(gap.tr.post - tau, gap.tr.pre, type) - if (!is.finite(s.tr)) return(FALSE) - .conformal_pval(s.tr, s.co, w.co = w.co, w.tr = w.tr) > alpha - }, logical(1)) - if (!any(accept)) return(c(NA_real_, NA_real_)) - rng <- range(grid[accept]) - ## flag if the acceptance set hit a grid edge (interval may be unbounded) - if (accept[1] || accept[n.grid]) { - attr(rng, "edge") <- TRUE + a <- accept_on(span) + ## adaptive expansion: if acceptance touches an edge the true interval may run + ## further (guards the span heuristic against unusual gap/score scales). + tries <- 0L + while (any(a$keep) && (a$keep[1L] || a$keep[length(a$keep)]) && tries < 3L) { + span <- span * 2; a <- accept_on(span); tries <- tries + 1L } + if (!any(a$keep)) { + out <- c(NA_real_, NA_real_) + attr(out, "empty") <- "rejected_all_tau" + return(out) + } + ## an accepted grid endpoint means the interval runs past it -> unbounded side. + lo <- min(a$grid[a$keep]); hi <- max(a$grid[a$keep]) + if (a$keep[1L]) lo <- -Inf + if (a$keep[length(a$keep)]) hi <- Inf + rng <- c(lo, hi) + if (is.infinite(lo) || is.infinite(hi)) attr(rng, "edge") <- TRUE rng } @@ -233,11 +270,25 @@ conformal_calibrate <- function(Y, D, X = NULL, I, II, T.on, r.cv, eff, ci <- .conformal_ci_grid(tr.post.path, tr.pre.path, S.co, score, alpha = alpha, w.co = w.co) } + ## status: "ok" | "unbounded" (alpha too small for Ncal) | "empty" (a + ## constant-effect score that rejects every tau -- correct conformal output, + ## not a failure; warn so callers/users know to read it as "no constant effect + ## consistent at level alpha", or switch to score = "meanabs"). + status <- "ok" + if (identical(attr(ci, "empty"), "rejected_all_tau")) { + status <- "empty" + warning("conformal: no constant effect is consistent at level ", alpha, + " under score = \"", score, "\" (the treated path is rejected at ", + "every tau). Reported interval is empty; score = \"meanabs\" gives ", + "a level interval for the average effect that cannot be empty.") + } else if (any(is.infinite(ci))) { + status <- "unbounded" + } S.tr <- .conformal_score(tr.post.path, tr.pre.path, score) p.value <- .conformal_pval(S.tr, S.co, w.co = w.co) - list(att = g.tr, ci = ci, p.value = p.value, - score.tr = S.tr, score.co = S.co, + list(att = g.tr, ci = c(ci[1], ci[2]), p.value = p.value, + score.tr = S.tr, score.co = S.co, status = status, n.calib = Ncal, n.eff = if (is.null(w.co)) Ncal else .conformal_neff(w.co), score = score, form = if (conformal.full) "full" else "jackknife+") } diff --git a/tests/testthat/test-conformal.R b/tests/testthat/test-conformal.R index 12846ff7..a93edeaa 100644 --- a/tests/testthat/test-conformal.R +++ b/tests/testthat/test-conformal.R @@ -1,7 +1,8 @@ -## Tests for the estimator-agnostic conformal core (R/conformal.R). -## The full vartype = "conformal" pipeline (leave-one-control-out calibration via -## impute_Y0) runs end to end but its output-slot integration and coverage -## validation are in progress; those tests are added with that wiring. +## Tests for the estimator-agnostic conformal core (R/conformal.R) AND the full +## leave-one-control-out calibration (conformal_calibrate, via impute_Y0). +## Coverage is validated here through the calibration itself; the fect() output-slot +## wiring (eff.calendar / est.* from vartype = "conformal") is still in progress and +## its print/plot tests are added with that wiring. test_that("conformal scores are finite and ordered as expected", { set.seed(1) @@ -49,3 +50,111 @@ test_that("effective sample size is sensible", { expect_equal(fect:::.conformal_neff(rep(1, 10)), 10) expect_lt(fect:::.conformal_neff(c(rep(1, 9), 100)), 10) # one big weight -> n_eff < n }) + +## ---- grid inversion (dispersion scores) ------------------------------------- + +test_that("grid denominator is outcome-scaled per score", { + e.pre <- rnorm(30) + expect_equal(fect:::.conformal_denom(e.pre, "studentized"), stats::sd(e.pre)) + expect_equal(fect:::.conformal_denom(e.pre, "ratio"), sqrt(mean(e.pre^2))) + expect_equal(fect:::.conformal_denom(e.pre, "rmse"), 1) +}) + +test_that("grid interval brackets a constant effect", { + set.seed(2) + s.co <- abs(rnorm(40, mean = 1, sd = 0.3)) # control scores + gap.post <- 3 + rnorm(8, sd = 0.3) # true constant effect = 3 + gap.pre <- rnorm(15, sd = 1) + ci <- fect:::.conformal_ci_grid(gap.post, gap.pre, s.co, "rmse", alpha = 0.10) + expect_true(all(is.finite(ci))) + expect_lt(ci[1], 3); expect_gt(ci[2], 3) +}) + +test_that("grid returns an explicit empty set when every tau is rejected", { + set.seed(3) + gap.post <- c(-10, 10, -10, 10, -10, 10) # huge dispersion: no tau fits + gap.pre <- rnorm(15, sd = 1) + s.co <- abs(rnorm(40, sd = 0.5)) # tight control scores + ci <- fect:::.conformal_ci_grid(gap.post, gap.pre, s.co, "rmse", alpha = 0.10) + expect_true(all(is.na(ci))) + expect_identical(attr(ci, "empty"), "rejected_all_tau") +}) + +test_that("grid is unbounded when alpha is below the resolution floor", { + set.seed(4) + gap.post <- rnorm(5, sd = 1); gap.pre <- rnorm(15, sd = 1) + s.co <- abs(rnorm(20)) + ci <- fect:::.conformal_ci_grid(gap.post, gap.pre, s.co, "rmse", alpha = 1 / 25) # < 1/21 + expect_true(any(is.infinite(ci))) +}) + +test_that("grid expands when the initial span is too tight", { + set.seed(7) + gap.post <- 5 + rnorm(6, sd = 0.3); gap.pre <- rnorm(15, sd = 1) + s.co <- abs(rnorm(40, mean = 1, sd = 0.3)) + ## start far inside the true center (~5); adaptive doubling must still recover it + ci <- fect:::.conformal_ci_grid(gap.post, gap.pre, s.co, "rmse", alpha = 0.10, span = 0.5) + expect_true(all(is.finite(ci))) + expect_lt(ci[1], 5); expect_gt(ci[2], 5) +}) + +## ---- full leave-one-control-out calibration (coverage) ---------------------- +## A block DGP with a known factor structure and zero true effect. Each rep fits +## gsynth, calibrates against the held-out donor scores, and records coverage. +## These are regression guards (no unbounded intervals; empties always flagged; +## coverage near nominal) -- precise calibration is the coverage study's job. + +.conf_dgp <- function(seed, N = 25, T = 20, T0 = 15, r = 2) { + set.seed(seed) + Fm <- matrix(rnorm(T * r), T, r); L <- matrix(rnorm(N * r), N, r) + D <- matrix(0, T, N); D[(T0 + 1):T, 1] <- 1 + list(mu = Fm %*% t(L), D = D, N = N, T = T) +} + +test_that("meanabs calibration covers near nominal and never empties/unbounds", { + skip_on_cran() + g <- .conf_dgp(20260608); reps <- 50; alpha <- 0.10 + cov <- logical(reps); bad <- 0L + for (b in seq_len(reps)) { + Y <- g$mu + matrix(rnorm(g$N * g$T), g$T, g$N) + dat <- data.frame(id = rep(1:g$N, each = g$T), time = rep(1:g$T, g$N), + Y = as.vector(Y), D = as.vector(g$D)) + fit <- suppressMessages(fect(Y ~ D, data = dat, index = c("id", "time"), + method = "gsynth", force = 3, CV = FALSE, r = 2, se = FALSE)) + cc <- conformal_calibrate(Y = Y, D = g$D, I = fit$I, II = fit$II, T.on = fit$T.on, + r.cv = fit$r.cv, eff = fit$eff, method = "gsynth", + score = "meanabs", alpha = alpha) + if (cc$status != "ok") bad <- bad + 1L + cov[b] <- (0 >= cc$ci[1]) && (0 <= cc$ci[2]) + } + expect_equal(bad, 0L) # meanabs: no empty, no unbounded + expect_gt(mean(cov), 0.82) # >= nominal, allowing MC noise over 50 reps + expect_lte(mean(cov), 1.0) +}) + +test_that("studentized calibration is calibrated and flags every empty", { + skip_on_cran() + ## 80 reps so the coverage bound is robust to Monte-Carlo noise (SE ~ 0.033 at + ## nominal 0.90); the bound is a gross-miscalibration guard, not a precise check. + g <- .conf_dgp(424242); reps <- 80; alpha <- 0.10 + cov.tot <- logical(reps); na.unflagged <- 0L; empties <- 0L + for (b in seq_len(reps)) { + Y <- g$mu + matrix(rnorm(g$N * g$T), g$T, g$N) + dat <- data.frame(id = rep(1:g$N, each = g$T), time = rep(1:g$T, g$N), + Y = as.vector(Y), D = as.vector(g$D)) + fit <- suppressMessages(fect(Y ~ D, data = dat, index = c("id", "time"), + method = "gsynth", force = 3, CV = FALSE, r = 2, se = FALSE)) + cc <- suppressWarnings(conformal_calibrate(Y = Y, D = g$D, I = fit$I, II = fit$II, + T.on = fit$T.on, r.cv = fit$r.cv, eff = fit$eff, method = "gsynth", + score = "studentized", alpha = alpha)) + if (any(is.na(cc$ci))) { + empties <- empties + 1L + if (cc$status != "empty") na.unflagged <- na.unflagged + 1L # never silent + } + ## total coverage: an empty set does NOT contain the truth -> FALSE + cov.tot[b] <- if (any(is.na(cc$ci))) FALSE else (0 >= cc$ci[1]) && (0 <= cc$ci[2]) + } + expect_equal(na.unflagged, 0L) # every empty interval is flagged + expect_lt(empties / reps, 0.30) # empties are the exception, not the rule + expect_gt(mean(cov.tot), 0.78) # total coverage (empties = miss) ~ nominal +}) From 695d1a24f07e35bf1508f5762d34b099ee3f907d Mon Sep 17 00:00:00 2001 From: Yiqing Xu Date: Mon, 8 Jun 2026 23:49:37 -0700 Subject: [PATCH 04/14] =?UTF-8?q?refactor(conformal):=20Family=20A=20only?= =?UTF-8?q?=20=E2=80=94=20level-statistic=20intervals,=20never=20empty?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Per design decision, remove the path / constant-effect scores (the only ones that could return an empty acceptance set). Every interval is now a level statistic S_i(tau) = |m_i - tau| / scale_i, inverted in closed form to m_tr +/- scale_tr*Q. ## Core (R/conformal.R, rewritten) - .conformal_center (mean | median), .conformal_scale (none | sd | rmspe | mad | diff; model-se supplied externally), .conformal_quantile (weighted (1-alpha) quantile with the treated +Inf atom; reduces to the ceil((1-alpha)(n+1))-th order statistic), .conformal_pval, .conformal_neff, .conformal_ci_level (closed form, never empty; unbounded only below the resolution floor). - Removed .conformal_score, .conformal_ci_meanabs, .conformal_denom, .conformal_ci_grid and the empty/status="empty" machinery. The grid-span and empty-handling fixes from earlier are retired with them: net simpler. - conformal_calibrate takes scale/center/weight (was score); returns status in {ok, unbounded}. LOO calibration plumbing unchanged. ## Dispatch (R/boot.R) - conformal branch calls scale = "none" (meanabs, the safe provisional default); still stop()s with a diagnostic, replaced by output-slot integration in Phase 2. ## Verification - tests/testthat/test-conformal.R rewritten for the Family-A API: 30 pass. - Coverage (200 reps, block DGP, effect 0, alpha=0.10), bad=0 throughout: scale=none 0.915/0.895; scale=sd 0.930/0.875 (two seeds). Both calibrated. - bootstrap & jackknife vartypes smoke-tested unaffected. Refs statsclaw runs/2026-06-08-conformal-inference.md (overnight build, Phase 1). Co-Authored-By: Claude Opus 4.8 --- R/boot.R | 26 ++- R/conformal.R | 297 +++++++++++++------------------- tests/testthat/test-conformal.R | 184 +++++++------------- 3 files changed, 196 insertions(+), 311 deletions(-) diff --git a/R/boot.R b/R/boot.R index 9c894fe5..682724b8 100644 --- a/R/boot.R +++ b/R/boot.R @@ -556,23 +556,21 @@ fect_boot <- function( method = method, predictive = time.component.from, force = force, hasRevs = hasRevs, tol = tol, max.iteration = max.iteration, norm.para = norm.para, - ## meanabs = closed-form level interval for the average effect: conservative, - ## never empty/unbounded (for Ncal above the resolution floor), and matches - ## fect's headline att.avg. The dispersion scores (studentized/ratio/rmse) - ## are constant-effect interval inversions -- correctly calibrated but can - ## return an empty set ~alpha of the time -- exposed later via a score arg. - score = "meanabs", alpha = 0.05 + ## Family A level interval (never empty/unbounded above the resolution + ## floor). scale = "none" is meanabs, the safe provisional default; the full + ## scale/center/weight knobs are threaded from fect() in Phase 2. + scale = "none", alpha = 0.05 ) - ## WIP (Phase 1 part 2): the calibration runs on the internal matrices and - ## yields the att.avg interval below. Wiring the result through fect()'s - ## output slots (eff.calendar, est.att, est.avg) must mirror the parametric - ## path's slot assembly (boot.R est.* construction), not an early return, - ## which skips the calendar-effect computation that fect.default expects. + ## WIP (Phase 1): the calibration runs on the internal matrices and yields the + ## att.avg interval below. Wiring the result through fect()'s output slots + ## (eff.calendar, est.att, est.avg) mirrors the parametric path's slot + ## assembly and replaces this stop() in Phase 2. stop(sprintf(paste0( "vartype = 'conformal': calibration is implemented but output integration ", - "is in progress. Result: att = %.4f, 95%% CI = [%.4f, %.4f], p = %.4f, ", - "N_calib = %d, score = %s, form = %s."), - cc$att, cc$ci[1], cc$ci[2], cc$p.value, cc$n.calib, cc$score, cc$form), + "is in progress. Result: att = %.4f, %.0f%% CI = [%.4f, %.4f], p = %.4f, ", + "N_calib = %d, scale = %s, status = %s."), + cc$att, 100 * (1 - 0.05), cc$ci[1], cc$ci[2], cc$p.value, cc$n.calib, + cc$scale, cc$status), call. = FALSE) } diff --git a/R/conformal.R b/R/conformal.R index c4f1056f..294f7e2e 100644 --- a/R/conformal.R +++ b/R/conformal.R @@ -2,62 +2,79 @@ ## conformal.R ## ## Cross-sectional conformal inference for separated, symmetric counterfactual -## estimators (vartype = "conformal"). The treated counterfactual prediction -## error is calibrated against the leave-one-control-out prediction errors of the -## donors and read off as a RANK, rather than resampled and normal-wrapped -## (parametric bootstrap) or quantile-resampled (nonparametric bootstrap). +## estimators (vartype = "conformal"). The treated unit's post-treatment +## prediction error is calibrated against the leave-one-control-out prediction +## errors of the donors and read off as a RANK, rather than resampled and +## normal-wrapped (parametric bootstrap) or quantile-resampled (nonparametric +## bootstrap). ## ## Design: statsclaw-workspace/fect/runs/REQ-conformal/spec.md -## Decisions (2026-06-08): jackknife+ default (>= 1 - 2 alpha); scores -## "studentized" (default), "ratio" (Abadie post/pre RMSPE), "meanabs", "rmse"; -## staggered adoption via valid_controls(predictive = "notyettreated"). +## Decision (2026-06-08): FAMILY A ONLY. Every interval is a level statistic ## -## This file holds the estimator-AGNOSTIC core: scores, the conformal p-value, -## and interval inversion. The calibration loop (leave-one-control-out via -## impute_Y0) and the vartype dispatch live in boot.R / default.R and feed this -## core a control-by-period residual matrix plus the treated gap path. +## S_i(tau) = | m_i - tau | / scale_i , interval m_tr +/- scale_tr * Q , +## +## where m_i is a location of unit i's post-period gaps (center: mean | median | +## per-horizon) and scale_i is a per-unit scale from its pre-period gaps +## (scale: none | sd | rmspe | mad | diff | model-se). Q is the conformal +## (1 - alpha) quantile of the control scores |m_j| / scale_j. This family is +## never empty (the numerator vanishes at tau = m_tr) and needs no grid inversion. +## The path / constant-effect scores (rmse / studentized-path / ratio), which +## could return an empty acceptance set, were removed because the empty case is +## confusing to report. +## +## This file is the estimator-AGNOSTIC core: centers, scales, the conformal +## quantile / p-value, and the closed-form interval. The calibration loop +## (leave-one-control-out via impute_Y0) and the vartype dispatch live in +## boot.R / default.R and feed this core the per-unit gap paths. ## ----------------------------------------------------------------------------- -## ---- 1. nonconformity score ------------------------------------------------- -## e.post : numeric, a unit's post-treatment residual (prediction-error) path -## e.pre : numeric, the same unit's pre-treatment residual path (for normaliz.) -## type : "studentized" | "ratio" | "meanabs" | "rmse" -## Returns a scalar. Validity (exchangeability) is identical across types; the -## choice affects interval width / adaptivity only. -.conformal_score <- function(e.post, e.pre, type = "studentized") { - e.post <- e.post[is.finite(e.post)] - if (length(e.post) == 0L) return(NA_real_) - rms <- function(x) sqrt(mean(x^2)) - if (type == "meanabs") { - return(abs(mean(e.post))) - } else if (type == "rmse") { - return(rms(e.post)) - } else if (type == "ratio") { - ## Abadie post/pre RMSPE ratio. - ep <- e.pre[is.finite(e.pre)] - denom <- if (length(ep) > 0L) rms(ep) else NA_real_ - if (is.na(denom) || denom <= 0) return(NA_real_) - return(rms(e.post) / denom) - } else if (type == "studentized") { - ep <- e.pre[is.finite(e.pre)] - denom <- if (length(ep) > 1L) stats::sd(ep) else NA_real_ - if (is.na(denom) || denom <= 0) return(NA_real_) - return(mean(abs(e.post)) / denom) - } - stop("conformal: unknown score type '", type, "'.") +## ---- 1. location (center) of a unit's post-period gap path ------------------- +## "per-horizon" is not a scalar; conformal_calibrate handles it by calling this +## per post period, so here we cover the scalar centers only. +.conformal_center <- function(e.post, type = "mean") { + e <- e.post[is.finite(e.post)] + if (length(e) == 0L) return(NA_real_) + if (type == "median") return(stats::median(e)) + mean(e) +} + +## ---- 2. per-unit scale from a unit's pre-period gap path --------------------- +## A common (unit-invariant) scale cancels in the ranking, so only per-unit +## variation matters; that is why "none" (meanabs) and a per-unit scale differ. +## "model-se" is not computable from e.pre alone -- conformal_calibrate supplies +## it per unit -- so this returns NA to signal the caller must override. +.conformal_scale <- function(e.pre, type = "none") { + if (type == "none") return(1) + if (type == "model-se") return(NA_real_) # supplied externally + e <- e.pre[is.finite(e.pre)] + if (type == "sd") return(if (length(e) > 1L) stats::sd(e) else NA_real_) + if (type == "rmspe") return(if (length(e) > 0L) sqrt(mean(e^2)) else NA_real_) + if (type == "mad") return(if (length(e) > 1L) stats::mad(e) else NA_real_) + if (type == "diff") return(if (length(e) > 1L) stats::sd(diff(e)) else NA_real_) + stop("conformal: unknown scale type '", type, "'.") +} + +## ---- 3. conformal quantile and p-value (weighted-capable) -------------------- +## Weighted (1 - alpha) quantile of the control scores, with the treated unit held +## as a +Inf atom of weight w.tr (Tibshirani et al. 2019). Unweighted, this is +## the ceil((1 - alpha)(n + 1))-th smallest control score. Returns Inf when the +## level is unreachable (alpha below the resolution floor 1 / (n + 1)). +.conformal_quantile <- function(s.co, alpha, w.co = NULL, w.tr = 1) { + ok <- is.finite(s.co); s <- s.co[ok] + n <- length(s); if (n == 0L) return(Inf) + w <- if (is.null(w.co)) rep(1, n) else w.co[ok] + ord <- order(s); s <- s[ord]; w <- w[ord] + cum <- cumsum(w) / (sum(w) + w.tr) + idx <- which(cum >= 1 - alpha) + if (length(idx) == 0L) return(Inf) # mass reached only at the +Inf atom + s[idx[1L]] } -## ---- 2. conformal p-value (weighted-capable) -------------------------------- -## s.tr : scalar treated score; s.co : numeric vector of control scores. -## w.co : optional non-negative control weights (weighted conformal); NULL = 1. -## w.tr : treated self-weight (default 1). Returns p in (0, 1]. +## p in (0, 1]: rank of the treated score among the controls. .conformal_pval <- function(s.tr, s.co, w.co = NULL, w.tr = 1) { - ok <- is.finite(s.co) - s.co <- s.co[ok] + ok <- is.finite(s.co); s.co <- s.co[ok] if (is.null(w.co)) w.co <- rep(1, length(s.co)) else w.co <- w.co[ok] - num <- w.tr + sum(w.co[s.co >= s.tr]) - den <- w.tr + sum(w.co) - num / den + (w.tr + sum(w.co[s.co >= s.tr])) / (w.tr + sum(w.co)) } ## effective sample size for weighted conformal @@ -67,105 +84,43 @@ (sum(w)^2) / sum(w^2) } -## ---- 3. interval inversion -------------------------------------------------- -## Closed form for the average effect under the "meanabs" score. -## g.tr : scalar treated mean post-period gap; g.co : control mean gaps. -## Returns c(lower, upper) for att.avg. Width = +/- the rank-corrected -## (1 - alpha) quantile of |g.co|. -.conformal_ci_meanabs <- function(g.tr, g.co, alpha = 0.05, w.co = NULL) { - ag <- abs(g.co[is.finite(g.co)]) - n <- length(ag) - k <- ceiling((1 - alpha) * (n + 1)) # conformal quantile rank - if (k > n) { # alpha < 1/(n+1): unbounded - return(c(-Inf, Inf)) - } - q <- sort(ag)[k] - c(g.tr - q, g.tr + q) -} - -## Outcome-scale denominator a given score divides by, used to size the grid in -## the SAME units as tau (the grid is over outcome-scale constant effects, so the -## span must be outcome-scaled, not score-scaled). -.conformal_denom <- function(e.pre, type) { - ep <- e.pre[is.finite(e.pre)] - if (type == "ratio") return(if (length(ep) > 0L) sqrt(mean(ep^2)) else NA_real_) - if (type == "studentized") return(if (length(ep) > 1L) stats::sd(ep) else NA_real_) - 1 # rmse / meanabs: numerator already in outcome units -} - -## Grid inversion for ratio / studentized / rmse (refit-free: control scores are -## tau-independent; only the treated score is recomputed on the grid). -## gap.tr.post / gap.tr.pre : treated gap paths. s.co : control scores. -## Returns c(lower, upper) = range of accepted constant effects tau. Two special -## returns, both carried as attributes so the caller can react: -## attr "empty" = "rejected_all_tau" : acceptance set is genuinely empty (no -## constant effect is consistent at level alpha; common with per-unit -## studentization when the treated pre-period sd is small) -> c(NA, NA). -## attr "edge" = TRUE : acceptance reached the (expanded) grid -## edge; the interval is effectively unbounded on that side. -.conformal_ci_grid <- function(gap.tr.post, gap.tr.pre, s.co, type, alpha = 0.05, - w.co = NULL, w.tr = 1, n.grid = 401L, span = NULL) { - gp <- gap.tr.post[is.finite(gap.tr.post)] - center <- mean(gp) - ## span in OUTCOME units: beyond |tau - center| ~ (max control score) * denom - ## the treated score must exceed every control score, so acceptance is - ## impossible -- this bounds where the grid needs to look. - denom <- .conformal_denom(gap.tr.pre, type) - if (!is.finite(denom) || denom <= 0) denom <- 1 - s.hi <- suppressWarnings(max(s.co[is.finite(s.co)])) - if (!is.finite(s.hi)) s.hi <- 1 - if (is.null(span)) span <- max(abs(gp)) + (s.hi + 1) * denom - - accept_on <- function(span) { - grid <- seq(center - span, center + span, length.out = n.grid) - keep <- vapply(grid, function(tau) { - s.tr <- .conformal_score(gap.tr.post - tau, gap.tr.pre, type) - if (!is.finite(s.tr)) return(FALSE) - .conformal_pval(s.tr, s.co, w.co = w.co, w.tr = w.tr) > alpha - }, logical(1)) - list(grid = grid, keep = keep) - } - a <- accept_on(span) - ## adaptive expansion: if acceptance touches an edge the true interval may run - ## further (guards the span heuristic against unusual gap/score scales). - tries <- 0L - while (any(a$keep) && (a$keep[1L] || a$keep[length(a$keep)]) && tries < 3L) { - span <- span * 2; a <- accept_on(span); tries <- tries + 1L - } - if (!any(a$keep)) { - out <- c(NA_real_, NA_real_) - attr(out, "empty") <- "rejected_all_tau" - return(out) - } - ## an accepted grid endpoint means the interval runs past it -> unbounded side. - lo <- min(a$grid[a$keep]); hi <- max(a$grid[a$keep]) - if (a$keep[1L]) lo <- -Inf - if (a$keep[length(a$keep)]) hi <- Inf - rng <- c(lo, hi) - if (is.infinite(lo) || is.infinite(hi)) attr(rng, "edge") <- TRUE - rng +## ---- 4. closed-form level interval ------------------------------------------ +## m.tr, scale.tr : treated center and per-unit scale. +## m.co, scale.co : control centers and per-unit scales (vectors). +## Control score S_j = |m_j| / scale_j (no effect under H0); treated +## S_tr(tau) = |m.tr - tau| / scale.tr. Acceptance {tau : S_tr(tau) <= Q} = +## m.tr +/- scale.tr * Q. NEVER empty; unbounded c(-Inf, Inf) only when alpha is +## below the resolution floor (too few controls for the requested level). +.conformal_ci_level <- function(m.tr, scale.tr, m.co, scale.co, alpha = 0.05, w.co = NULL) { + ok <- is.finite(m.co) & is.finite(scale.co) & scale.co > 0 + s.co <- abs(m.co[ok]) / scale.co[ok] + w <- if (is.null(w.co)) NULL else w.co[ok] + Q <- .conformal_quantile(s.co, alpha, w.co = w) + if (!is.finite(Q) || !is.finite(scale.tr) || scale.tr <= 0) return(c(-Inf, Inf)) + half <- scale.tr * Q + c(m.tr - half, m.tr + half) } -## ---- 4. calibration: deterministic leave-one-control-out -------------------- +## ---- 5. calibration: deterministic leave-one-control-out -------------------- ## Mirror of boot.R::draw.error but HELD-OUT (not resampled): each valid control ## is predicted once from the other controls, giving its out-of-fold residual -## path. Estimator-agnostic — built-ins via impute_Y0(method); a custom separated +## path. Estimator-agnostic -- built-ins via impute_Y0(method); a custom separated ## learner via `conformal.fit`. Inputs are the preprocessed matrices available ## inside fect_boot(): Y, D, X, I, II, T.on (all TT x N; X is TT x N x p or NULL), ## r.cv (selected rank), and eff (TT x N point-fit gap matrix). ## -## Returns: list(att, ci, p.value, score.tr, score.co, n.calib, n.eff, score, form). +## Returns: list(att, ci, p.value, score.tr, score.co, status, n.calib, n.eff, +## scale, center, weight, form). ## Scope (verified): block design (common onset), method gsynth/ife, nevertreated ## (separated) calibration. Staggered uses the union post-window (per-cohort -## windows are a TODO); pooled (notyettreated) fits fall back to nevertreated with -## a warning. MC is not yet routed through impute_Y0 (upstream stop()). +## windows are a TODO, Phase 4); pooled (notyettreated) fits fall back to +## nevertreated with a warning. MC is not yet routed through impute_Y0. conformal_calibrate <- function(Y, D, X = NULL, I, II, T.on, r.cv, eff, method = "gsynth", predictive = "nevertreated", force = 3L, hasRevs = 0L, tol = 1e-5, max.iteration = 1000L, norm.para = NULL, - score = "studentized", alpha = 0.05, - conformal.full = FALSE, conformal.weight = "none", - conformal.fit = NULL) { + scale = "none", center = "mean", weight = "cell", + alpha = 0.05, conformal.fit = NULL) { TT <- nrow(Y); N <- ncol(Y) sum.D <- colSums(D) @@ -201,12 +156,10 @@ conformal_calibrate <- function(Y, D, X = NULL, I, II, T.on, r.cv, eff, } } - ## --- one held-out fit per valid control, reusing the parametric refit path - d.pattern <- D[, id.tr[1]] # a treated D column (defines onset) - sub3 <- function(A, idx) if (is.null(A)) NULL else A[, idx, , drop = FALSE] - + ## --- one held-out fit per valid control, reusing the parametric refit path. ## fake-treated column = held-out control j's DATA carrying a treated unit's ## TIMING (D, T.on, II must agree with the assigned onset; only Y/I are j's). + d.pattern <- D[, id.tr[1]] # a treated D column (defines onset) loo_gap <- function(j) { co.rest <- setdiff(valid.co, j) if (!is.null(conformal.fit)) { @@ -235,60 +188,44 @@ conformal_calibrate <- function(Y, D, X = NULL, I, II, T.on, r.cv, eff, g } - ## --- control scores (tau-independent) and mean post gaps + ## --- control centers and per-unit scales (tau-independent) obs <- (I == 1) - S.co <- rep(NA_real_, length(valid.co)) - g.co <- rep(NA_real_, length(valid.co)) - for (m in seq_along(valid.co)) { - j <- valid.co[m] + m.co <- rep(NA_real_, length(valid.co)) + sc.co <- rep(NA_real_, length(valid.co)) + for (mi in seq_along(valid.co)) { + j <- valid.co[mi] gj <- loo_gap(j) ep <- gj[post.idx][obs[post.idx, j]] eq <- gj[pre.idx ][obs[pre.idx, j]] - S.co[m] <- .conformal_score(ep, eq, score) - g.co[m] <- mean(ep, na.rm = TRUE) + m.co[mi] <- .conformal_center(ep, center) + sc.co[mi] <- .conformal_scale(eq, scale) } - keep <- is.finite(S.co) - S.co <- S.co[keep]; g.co <- g.co[keep] - Ncal <- length(S.co) - - ## --- treated statistic (averaged over treated units, post window) - eff.tr.mat <- eff[, id.tr, drop = FALSE] + keep <- is.finite(m.co) & is.finite(sc.co) & sc.co > 0 + m.co <- m.co[keep]; sc.co <- sc.co[keep] + Ncal <- length(m.co) + + ## --- treated aggregate (cell-weighted average over treated units, by default). + ## NOTE: unit / precision weighting and the per-period band land in Phase 3/4; + ## for the block design rowMeans over treated units is the cell-weighted center. + eff.tr.mat <- eff[, id.tr, drop = FALSE] tr.post.path <- rowMeans(eff.tr.mat[post.idx, , drop = FALSE], na.rm = TRUE) tr.pre.path <- rowMeans(eff.tr.mat[pre.idx, , drop = FALSE], na.rm = TRUE) - g.tr <- mean(tr.post.path, na.rm = TRUE) + m.tr <- .conformal_center(tr.post.path, center) + sc.tr <- .conformal_scale(tr.pre.path, scale) - ## --- weights (overlap) — placeholder until loading_bound wiring (Phase 3) + ## --- weights (overlap density ratio) deferred to a later phase w.co <- NULL - if (identical(conformal.weight, "overlap")) { - warning("conformal.weight = \"overlap\" not yet wired; using unweighted.") - } - ## --- interval - if (score == "meanabs") { - ci <- .conformal_ci_meanabs(g.tr, g.co, alpha = alpha, w.co = w.co) - } else { - ci <- .conformal_ci_grid(tr.post.path, tr.pre.path, S.co, score, - alpha = alpha, w.co = w.co) - } - ## status: "ok" | "unbounded" (alpha too small for Ncal) | "empty" (a - ## constant-effect score that rejects every tau -- correct conformal output, - ## not a failure; warn so callers/users know to read it as "no constant effect - ## consistent at level alpha", or switch to score = "meanabs"). - status <- "ok" - if (identical(attr(ci, "empty"), "rejected_all_tau")) { - status <- "empty" - warning("conformal: no constant effect is consistent at level ", alpha, - " under score = \"", score, "\" (the treated path is rejected at ", - "every tau). Reported interval is empty; score = \"meanabs\" gives ", - "a level interval for the average effect that cannot be empty.") - } else if (any(is.infinite(ci))) { - status <- "unbounded" - } - S.tr <- .conformal_score(tr.post.path, tr.pre.path, score) - p.value <- .conformal_pval(S.tr, S.co, w.co = w.co) + ## --- interval (closed form; never empty) + ci <- .conformal_ci_level(m.tr, sc.tr, m.co, sc.co, alpha = alpha, w.co = w.co) + status <- if (any(is.infinite(ci))) "unbounded" else "ok" + + s.tr <- if (is.finite(sc.tr) && sc.tr > 0) abs(m.tr) / sc.tr else NA_real_ + s.co.vec <- abs(m.co) / sc.co + p.value <- .conformal_pval(s.tr, s.co.vec, w.co = w.co) - list(att = g.tr, ci = c(ci[1], ci[2]), p.value = p.value, - score.tr = S.tr, score.co = S.co, status = status, + list(att = m.tr, ci = c(ci[1], ci[2]), p.value = p.value, + score.tr = s.tr, score.co = s.co.vec, status = status, n.calib = Ncal, n.eff = if (is.null(w.co)) Ncal else .conformal_neff(w.co), - score = score, form = if (conformal.full) "full" else "jackknife+") + scale = scale, center = center, weight = weight, form = "jackknife+") } diff --git a/tests/testthat/test-conformal.R b/tests/testthat/test-conformal.R index a93edeaa..a516c4cf 100644 --- a/tests/testthat/test-conformal.R +++ b/tests/testthat/test-conformal.R @@ -1,108 +1,68 @@ -## Tests for the estimator-agnostic conformal core (R/conformal.R) AND the full +## Tests for the Family-A conformal core (R/conformal.R) and the full ## leave-one-control-out calibration (conformal_calibrate, via impute_Y0). -## Coverage is validated here through the calibration itself; the fect() output-slot -## wiring (eff.calendar / est.* from vartype = "conformal") is still in progress and -## its print/plot tests are added with that wiring. +## Every Family-A interval is a level statistic |m - tau| / scale: closed-form and +## never empty. Coverage is validated through the calibration itself; the fect() +## output-slot wiring (vartype = "conformal") gets its print/plot tests once that +## lands. -test_that("conformal scores are finite and ordered as expected", { - set.seed(1) - e.post <- rnorm(12); e.pre <- rnorm(24) - for (s in c("meanabs", "rmse", "ratio", "studentized")) { - expect_true(is.finite(fect:::.conformal_score(e.post, e.pre, s))) - } - ## a larger post deviation gives a larger score (meanabs) - expect_gt(fect:::.conformal_score(e.post + 5, e.pre, "meanabs"), - fect:::.conformal_score(e.post, e.pre, "meanabs")) - ## ratio with zero pre-variation is NA, not Inf-by-accident - expect_true(is.na(fect:::.conformal_score(e.post, rep(0, 24), "ratio"))) +test_that("center: mean and median", { + set.seed(1); e <- rnorm(8) + expect_equal(fect:::.conformal_center(e, "mean"), mean(e)) + expect_equal(fect:::.conformal_center(e, "median"), stats::median(e)) + expect_true(is.na(fect:::.conformal_center(numeric(0), "mean"))) +}) + +test_that("scale: none / sd / rmspe / mad / diff (+ model-se is external)", { + set.seed(2); e <- rnorm(20) + expect_equal(fect:::.conformal_scale(e, "none"), 1) + expect_equal(fect:::.conformal_scale(e, "sd"), stats::sd(e)) + expect_equal(fect:::.conformal_scale(e, "rmspe"), sqrt(mean(e^2))) + expect_equal(fect:::.conformal_scale(e, "mad"), stats::mad(e)) + expect_equal(fect:::.conformal_scale(e, "diff"), stats::sd(diff(e))) + expect_true(is.na(fect:::.conformal_scale(e, "model-se"))) # supplied externally }) -test_that("meanabs closed-form interval covers near nominal", { - set.seed(20260608) - reps <- 3000; alpha <- 0.10; Nco <- 39; att <- 3 - cov <- replicate(reps, { - g.co <- rnorm(Nco) # no-effect donor mean gaps - g.tr <- att + rnorm(1) # treated mean gap = att + exchangeable noise - ci <- fect:::.conformal_ci_meanabs(g.tr, g.co, alpha = alpha) - (att >= ci[1]) && (att <= ci[2]) - }) - expect_gt(mean(cov), 0.86) - expect_lt(mean(cov), 0.94) # ~ 1 - alpha +test_that("conformal quantile reduces to the order statistic and guards resolution", { + set.seed(3); s <- abs(rnorm(20)) + for (a in c(0.05, 0.10, 0.20)) { + k <- ceiling((1 - a) * (length(s) + 1)) + expected <- if (k > length(s)) Inf else sort(s)[k] + expect_equal(fect:::.conformal_quantile(s, a), expected) + } + ## alpha below 1/(n+1) -> Inf (cannot reach the level) + expect_true(is.infinite(fect:::.conformal_quantile(s, 1 / 30))) # 1/30 < 1/21 }) test_that("conformal p-value lies in (0,1] and respects weights", { s.co <- rnorm(20) - p <- fect:::.conformal_pval(0, s.co) - expect_gt(p, 0); expect_lte(p, 1) - ## a very large treated score is the most extreme -> smallest p = 1/(n+1) + p <- fect:::.conformal_pval(0, s.co); expect_gt(p, 0); expect_lte(p, 1) expect_equal(fect:::.conformal_pval(100, s.co), 1 / (length(s.co) + 1)) - ## zero weight on a donor drops it from the count w <- rep(1, 20); w[which.max(s.co)] <- 0 expect_true(is.finite(fect:::.conformal_pval(0, s.co, w.co = w))) }) -test_that("resolution guard: alpha below 1/(Nco+1) yields an unbounded interval", { - ci <- fect:::.conformal_ci_meanabs(0, rnorm(20), alpha = 1 / 25) # 1/25 < 1/21 - expect_true(all(is.infinite(ci))) +test_that("level interval brackets the center, scales with the treated scale, never empty", { + set.seed(4) + m.co <- rnorm(40); sc.co <- rep(1, 40) + ci1 <- fect:::.conformal_ci_level(2, 1, m.co, sc.co, alpha = 0.10) + expect_true(all(is.finite(ci1))) + expect_lt(ci1[1], 2); expect_gt(ci1[2], 2) # brackets the center + ci2 <- fect:::.conformal_ci_level(2, 2, m.co, sc.co, alpha = 0.10) + expect_equal(diff(ci2), 2 * diff(ci1)) # half-width = scale.tr * Q + ## below the resolution floor -> unbounded, never empty + ci3 <- fect:::.conformal_ci_level(2, 1, rnorm(15), rep(1, 15), alpha = 1 / 30) + expect_true(all(is.infinite(ci3))) }) test_that("effective sample size is sensible", { expect_equal(fect:::.conformal_neff(rep(1, 10)), 10) - expect_lt(fect:::.conformal_neff(c(rep(1, 9), 100)), 10) # one big weight -> n_eff < n -}) - -## ---- grid inversion (dispersion scores) ------------------------------------- - -test_that("grid denominator is outcome-scaled per score", { - e.pre <- rnorm(30) - expect_equal(fect:::.conformal_denom(e.pre, "studentized"), stats::sd(e.pre)) - expect_equal(fect:::.conformal_denom(e.pre, "ratio"), sqrt(mean(e.pre^2))) - expect_equal(fect:::.conformal_denom(e.pre, "rmse"), 1) -}) - -test_that("grid interval brackets a constant effect", { - set.seed(2) - s.co <- abs(rnorm(40, mean = 1, sd = 0.3)) # control scores - gap.post <- 3 + rnorm(8, sd = 0.3) # true constant effect = 3 - gap.pre <- rnorm(15, sd = 1) - ci <- fect:::.conformal_ci_grid(gap.post, gap.pre, s.co, "rmse", alpha = 0.10) - expect_true(all(is.finite(ci))) - expect_lt(ci[1], 3); expect_gt(ci[2], 3) -}) - -test_that("grid returns an explicit empty set when every tau is rejected", { - set.seed(3) - gap.post <- c(-10, 10, -10, 10, -10, 10) # huge dispersion: no tau fits - gap.pre <- rnorm(15, sd = 1) - s.co <- abs(rnorm(40, sd = 0.5)) # tight control scores - ci <- fect:::.conformal_ci_grid(gap.post, gap.pre, s.co, "rmse", alpha = 0.10) - expect_true(all(is.na(ci))) - expect_identical(attr(ci, "empty"), "rejected_all_tau") -}) - -test_that("grid is unbounded when alpha is below the resolution floor", { - set.seed(4) - gap.post <- rnorm(5, sd = 1); gap.pre <- rnorm(15, sd = 1) - s.co <- abs(rnorm(20)) - ci <- fect:::.conformal_ci_grid(gap.post, gap.pre, s.co, "rmse", alpha = 1 / 25) # < 1/21 - expect_true(any(is.infinite(ci))) -}) - -test_that("grid expands when the initial span is too tight", { - set.seed(7) - gap.post <- 5 + rnorm(6, sd = 0.3); gap.pre <- rnorm(15, sd = 1) - s.co <- abs(rnorm(40, mean = 1, sd = 0.3)) - ## start far inside the true center (~5); adaptive doubling must still recover it - ci <- fect:::.conformal_ci_grid(gap.post, gap.pre, s.co, "rmse", alpha = 0.10, span = 0.5) - expect_true(all(is.finite(ci))) - expect_lt(ci[1], 5); expect_gt(ci[2], 5) + expect_lt(fect:::.conformal_neff(c(rep(1, 9), 100)), 10) }) ## ---- full leave-one-control-out calibration (coverage) ---------------------- -## A block DGP with a known factor structure and zero true effect. Each rep fits -## gsynth, calibrates against the held-out donor scores, and records coverage. -## These are regression guards (no unbounded intervals; empties always flagged; -## coverage near nominal) -- precise calibration is the coverage study's job. +## Block DGP, known factor structure, zero true effect. Regression guards: never +## empty, never unbounded (Ncal above the floor), coverage near nominal. Precise +## calibration is the coverage study's job. .conf_dgp <- function(seed, N = 25, T = 20, T0 = 15, r = 2) { set.seed(seed) @@ -111,10 +71,8 @@ test_that("grid expands when the initial span is too tight", { list(mu = Fm %*% t(L), D = D, N = N, T = T) } -test_that("meanabs calibration covers near nominal and never empties/unbounds", { - skip_on_cran() - g <- .conf_dgp(20260608); reps <- 50; alpha <- 0.10 - cov <- logical(reps); bad <- 0L +.conf_cover <- function(seed, scale, reps = 80, alpha = 0.10) { + g <- .conf_dgp(seed); cov <- logical(reps); bad <- 0L for (b in seq_len(reps)) { Y <- g$mu + matrix(rnorm(g$N * g$T), g$T, g$N) dat <- data.frame(id = rep(1:g$N, each = g$T), time = rep(1:g$T, g$N), @@ -123,38 +81,30 @@ test_that("meanabs calibration covers near nominal and never empties/unbounds", method = "gsynth", force = 3, CV = FALSE, r = 2, se = FALSE)) cc <- conformal_calibrate(Y = Y, D = g$D, I = fit$I, II = fit$II, T.on = fit$T.on, r.cv = fit$r.cv, eff = fit$eff, method = "gsynth", - score = "meanabs", alpha = alpha) + scale = scale, alpha = alpha) if (cc$status != "ok") bad <- bad + 1L cov[b] <- (0 >= cc$ci[1]) && (0 <= cc$ci[2]) } - expect_equal(bad, 0L) # meanabs: no empty, no unbounded - expect_gt(mean(cov), 0.82) # >= nominal, allowing MC noise over 50 reps - expect_lte(mean(cov), 1.0) + list(cov = mean(cov), bad = bad) +} + +test_that("meanabs (scale = none) covers near nominal, never empty/unbounded", { + skip_on_cran() + ## scale = none is stable (~0.90-0.92 at 200 reps); loose bound guards gross + ## miscalibration without flaking on Monte-Carlo noise. Precise calibration is + ## the Phase 5 simulation study, not this unit test. + r <- .conf_cover(20260608, "none") + expect_equal(r$bad, 0L) + expect_gt(r$cov, 0.80); expect_lte(r$cov, 1.0) }) -test_that("studentized calibration is calibrated and flags every empty", { +test_that("studentized-mean (scale = sd) runs, never empty/unbounded, ballpark coverage", { skip_on_cran() - ## 80 reps so the coverage bound is robust to Monte-Carlo noise (SE ~ 0.033 at - ## nominal 0.90); the bound is a gross-miscalibration guard, not a precise check. - g <- .conf_dgp(424242); reps <- 80; alpha <- 0.10 - cov.tot <- logical(reps); na.unflagged <- 0L; empties <- 0L - for (b in seq_len(reps)) { - Y <- g$mu + matrix(rnorm(g$N * g$T), g$T, g$N) - dat <- data.frame(id = rep(1:g$N, each = g$T), time = rep(1:g$T, g$N), - Y = as.vector(Y), D = as.vector(g$D)) - fit <- suppressMessages(fect(Y ~ D, data = dat, index = c("id", "time"), - method = "gsynth", force = 3, CV = FALSE, r = 2, se = FALSE)) - cc <- suppressWarnings(conformal_calibrate(Y = Y, D = g$D, I = fit$I, II = fit$II, - T.on = fit$T.on, r.cv = fit$r.cv, eff = fit$eff, method = "gsynth", - score = "studentized", alpha = alpha)) - if (any(is.na(cc$ci))) { - empties <- empties + 1L - if (cc$status != "empty") na.unflagged <- na.unflagged + 1L # never silent - } - ## total coverage: an empty set does NOT contain the truth -> FALSE - cov.tot[b] <- if (any(is.na(cc$ci))) FALSE else (0 >= cc$ci[1]) && (0 <= cc$ci[2]) - } - expect_equal(na.unflagged, 0L) # every empty interval is flagged - expect_lt(empties / reps, 0.30) # empties are the exception, not the rule - expect_gt(mean(cov.tot), 0.78) # total coverage (empties = miss) ~ nominal + ## scale = sd is calibrated on average (~0.90) but has higher coverage variance + ## than none, because its width is random via the treated pre-period sd. The + ## strict assertion is bad == 0 (Family A can never empty); coverage is a wide + ## gross-guard here, with the precise check deferred to the Phase 5 study. + r <- .conf_cover(424242, "sd") + expect_equal(r$bad, 0L) + expect_gt(r$cov, 0.75); expect_lte(r$cov, 1.0) }) From 5e06189f80c2aab8119a3cb2ff33eea258b51fb3 Mon Sep 17 00:00:00 2001 From: Yiqing Xu Date: Tue, 9 Jun 2026 00:08:26 -0700 Subject: [PATCH 05/14] feat(conformal): make vartype="conformal" usable via fect() (Phase 2) Output-slot integration + argument threading. fect(vartype="conformal") now returns a complete, printable, plottable object; no bootstrap draws. ## conformal_calibrate - Also returns a per-period pointwise band (calendar-indexed) built from the full leave-one-control-out gap matrix: pre rows = placebo/pre-trend band, post rows = effect band. Guards model-se scale and non-cell weight as "not yet implemented". ## boot.R conformal branch - Builds est.avg / est.avg.unit / est.att (event-time via the treated T.on map) / est.att90 / att.bound / est.eff.calendar(.fit) + a `conformal` meta slot, then returns c(out, result). Nominal S.E. back-derived from the symmetric CI (display only). Unsupported options (group/reversals/W/balance/placebo/carryover) error. ## default.R - Thread conformal.scale/center/weight through fect / fect.formula / fect.default and both pass-throughs (+ fect_boot signature); validate the values. - Skip the bootstrap-based diagtest (F/equivalence test) for conformal, which has no draws (consistent with se=FALSE); fixes an apply(!is.na(att.boot)) error. ## Verification - 40 tests pass (test-conformal.R), incl. 2 fect() end-to-end integration tests. - print() and plot(gap|counterfactual) work; 5 scales (none/sd/rmspe/mad/diff) + median center give distinct finite intervals; CI covers the true ATT. - bootstrap / jackknife / parametric unaffected; their test.out still produced. Refs statsclaw runs/2026-06-08-conformal-inference.md (overnight build, Phase 2). Co-Authored-By: Claude Opus 4.8 --- R/boot.R | 100 +++++++++++++++++++++++++------- R/conformal.R | 80 ++++++++++++++++--------- R/default.R | 35 ++++++++++- tests/testthat/test-conformal.R | 40 +++++++++++++ 4 files changed, 206 insertions(+), 49 deletions(-) diff --git a/R/boot.R b/R/boot.R index 682724b8..ebb7e935 100644 --- a/R/boot.R +++ b/R/boot.R @@ -142,6 +142,9 @@ fect_boot <- function( carryover.period = NULL, vartype = "bootstrap", para.error = "auto", + conformal.scale = "none", + conformal.center = "mean", + conformal.weight = "cell", quantile.CI = FALSE, nboots = 200, parallel = TRUE, @@ -544,34 +547,91 @@ fect_boot <- function( N_unit <- dim(out$res)[2] ## ---- vartype = "conformal": cross-sectional conformal interval ----------- - ## Rank the treated gap against the leave-one-control-out donor scores - ## (see conformal.R), bypassing the resampling / jackknife SE machinery. - ## Uses fect_boot's internal preprocessed matrices (Y, D, I, II, T.on). - ## NOTE Phase 1: score / alpha hardcoded; est.att per-period bands + arg - ## threading from fect() come next. + ## Rank the treated gaps against the leave-one-control-out donor gaps (Family A + ## level statistic, see conformal.R), bypassing the resampling / jackknife SE + ## machinery. Populate the est.* slots so print / plot / esplot work unchanged, + ## then return early (no bootstrap draws). if (vartype == "conformal") { + ## Phase 2 scope is the core SC case; unsupported options error clearly (they + ## are added in later phases). Reversals, weights, balance, placebo, carryover + ## and group estimates still go through bootstrap / jackknife. + if (!is.null(group)) { + stop("vartype = 'conformal' does not yet support group-level estimates; ", + "use vartype = 'bootstrap' or 'jackknife'.", call. = FALSE) + } + if (!is.null(balance.period) || !is.null(W) || hasRevs == 1 || + isTRUE(placeboTest) || isTRUE(carryoverTest)) { + stop("vartype = 'conformal' does not yet support reversals, weights, ", + "balanced-panel, placebo or carryover tests; use vartype = ", + "'bootstrap' or 'jackknife' for those.", call. = FALSE) + } + + id.tr.conf <- which(colSums(D) > 0) cc <- conformal_calibrate( Y = Y, D = D, X = X, I = I, II = II, T.on = T.on, r.cv = out$r.cv, eff = out$eff, method = method, predictive = time.component.from, force = force, hasRevs = hasRevs, tol = tol, max.iteration = max.iteration, norm.para = norm.para, - ## Family A level interval (never empty/unbounded above the resolution - ## floor). scale = "none" is meanabs, the safe provisional default; the full - ## scale/center/weight knobs are threaded from fect() in Phase 2. - scale = "none", alpha = 0.05 + scale = conformal.scale, center = conformal.center, + weight = conformal.weight, alpha = alpha + ) + + ## back out a nominal S.E. from the symmetric conformal CI (display only; + ## conformal inference is rank-based, the S.E. column is informational). + z <- stats::qnorm(1 - alpha / 2) + se.from <- function(lo, hi) ifelse(is.finite(hi - lo), (hi - lo) / (2 * z), NA_real_) + + ## --- average effect -> est.avg / est.avg.unit + se.avg <- se.from(cc$ci[1], cc$ci[2]) + est.avg <- t(as.matrix(c(att.avg, se.avg, cc$ci[1], cc$ci[2], cc$p.value))) + colnames(est.avg) <- c("ATT.avg", "S.E.", "CI.lower", "CI.upper", "p.value") + att.avg.unit.val <- if (length(att.avg.unit)) att.avg.unit else att.avg + est.avg.unit <- t(as.matrix(c(att.avg.unit.val, se.avg, cc$ci[1], cc$ci[2], cc$p.value))) + colnames(est.avg.unit) <- c("ATT.avg.unit", "S.E.", "CI.lower", "CI.upper", "p.value") + + ## --- per-period band (calendar-indexed) -> est.eff.calendar(.fit) directly + band <- cc$band + se.cal <- se.from(band[, "CI.lower"], band[, "CI.upper"]) + est.eff.calendar <- cbind(calendar.eff, se.cal, + band[, "CI.lower"], band[, "CI.upper"], band[, "p.value"], calendar.N) + colnames(est.eff.calendar) <- c("ATT-calendar", "S.E.", "CI.lower", "CI.upper", "p.value", "count") + est.eff.calendar.fit <- cbind(calendar.eff.fit, se.cal, + band[, "CI.lower"], band[, "CI.upper"], band[, "p.value"], calendar.N) + colnames(est.eff.calendar.fit) <- colnames(est.eff.calendar) + + ## --- map the calendar band to event time (treated unit's relative period) + ## for est.att (rownames = out$time). Common-onset: a 1:1 map. + rel <- T.on[, id.tr.conf[1]] + est.att <- matrix(NA_real_, length(time.on), 6, + dimnames = list(time.on, + c("ATT", "S.E.", "CI.lower", "CI.upper", "p.value", "count"))) + for (k in seq_along(time.on)) { + rows <- which(rel == time.on[k]) + if (length(rows) == 0L) next + lo <- mean(band[rows, "CI.lower"], na.rm = TRUE) + hi <- mean(band[rows, "CI.upper"], na.rm = TRUE) + pv <- mean(band[rows, "p.value"], na.rm = TRUE) + est.att[k, ] <- c(att[k], se.from(lo, hi), lo, hi, pv, out$count[k]) + } + ## Phase 4 builds a true (1 - 2*alpha) inner band and the simultaneous band; + ## for now the 90%-slot mirrors the main band so plot paths do not break. + est.att90 <- est.att + att.bound <- est.att[, c("CI.lower", "CI.upper"), drop = FALSE] + rownames(att.bound) <- time.on + + result <- list( + est.avg = est.avg, + est.avg.unit = est.avg.unit, + est.att = est.att, + est.att90 = est.att90, + att.bound = att.bound, + est.eff.calendar = est.eff.calendar, + est.eff.calendar.fit = est.eff.calendar.fit, + vartype = "conformal", + conformal = cc[c("scale", "center", "weight", "status", "n.calib", "form")] ) - ## WIP (Phase 1): the calibration runs on the internal matrices and yields the - ## att.avg interval below. Wiring the result through fect()'s output slots - ## (eff.calendar, est.att, est.avg) mirrors the parametric path's slot - ## assembly and replaces this stop() in Phase 2. - stop(sprintf(paste0( - "vartype = 'conformal': calibration is implemented but output integration ", - "is in progress. Result: att = %.4f, %.0f%% CI = [%.4f, %.4f], p = %.4f, ", - "N_calib = %d, scale = %s, status = %s."), - cc$att, 100 * (1 - 0.05), cc$ci[1], cc$ci[2], cc$p.value, cc$n.calib, - cc$scale, cc$status), - call. = FALSE) + return(c(out, result)) } if (!is.null(group)) { diff --git a/R/conformal.R b/R/conformal.R index 294f7e2e..ebfce0c8 100644 --- a/R/conformal.R +++ b/R/conformal.R @@ -128,6 +128,16 @@ conformal_calibrate <- function(Y, D, X = NULL, I, II, T.on, r.cv, eff, Ntr <- length(id.tr); Nco <- length(id.co) if (Nco < 2L) stop("conformal: need at least 2 controls.") + ## not-yet-implemented knobs (land in a later phase): fail loudly, not silently. + if (identical(scale, "model-se")) { + stop("conformal: scale = \"model-se\" is not yet implemented; use one of ", + "\"none\", \"sd\", \"rmspe\", \"mad\", \"diff\".", call. = FALSE) + } + if (!identical(weight, "cell")) { + stop("conformal: weight = \"", weight, "\" is not yet implemented; ", + "use weight = \"cell\" (the per-treated-cell ATT).", call. = FALSE) + } + ## --- separation guard: conformal calibration is controls-only (nevertreated) if (!identical(predictive, "nevertreated")) { warning("vartype = \"conformal\" requires a separated fit; ", @@ -188,44 +198,60 @@ conformal_calibrate <- function(Y, D, X = NULL, I, II, T.on, r.cv, eff, g } - ## --- control centers and per-unit scales (tau-independent) + ## --- full per-control LOO gap paths (calendar-indexed), and a per-unit scale + ## from each control's pre-period gaps. Keeping the whole path lets us serve + ## both the scalar average interval AND the per-period band from one pass. obs <- (I == 1) - m.co <- rep(NA_real_, length(valid.co)) + G.co <- matrix(NA_real_, length(valid.co), TT) # control x calendar gap sc.co <- rep(NA_real_, length(valid.co)) for (mi in seq_along(valid.co)) { j <- valid.co[mi] gj <- loo_gap(j) - ep <- gj[post.idx][obs[post.idx, j]] - eq <- gj[pre.idx ][obs[pre.idx, j]] - m.co[mi] <- .conformal_center(ep, center) - sc.co[mi] <- .conformal_scale(eq, scale) + gj[!obs[, j]] <- NA_real_ + G.co[mi, ] <- gj + sc.co[mi] <- .conformal_scale(gj[pre.idx], scale) } - keep <- is.finite(m.co) & is.finite(sc.co) & sc.co > 0 - m.co <- m.co[keep]; sc.co <- sc.co[keep] - Ncal <- length(m.co) - - ## --- treated aggregate (cell-weighted average over treated units, by default). - ## NOTE: unit / precision weighting and the per-period band land in Phase 3/4; - ## for the block design rowMeans over treated units is the cell-weighted center. - eff.tr.mat <- eff[, id.tr, drop = FALSE] - tr.post.path <- rowMeans(eff.tr.mat[post.idx, , drop = FALSE], na.rm = TRUE) - tr.pre.path <- rowMeans(eff.tr.mat[pre.idx, , drop = FALSE], na.rm = TRUE) - m.tr <- .conformal_center(tr.post.path, center) - sc.tr <- .conformal_scale(tr.pre.path, scale) - - ## --- weights (overlap density ratio) deferred to a later phase - w.co <- NULL - - ## --- interval (closed form; never empty) - ci <- .conformal_ci_level(m.tr, sc.tr, m.co, sc.co, alpha = alpha, w.co = w.co) - status <- if (any(is.infinite(ci))) "unbounded" else "ok" + ## --- treated aggregate path (cell-weighted average over treated units). + ## NOTE: unit / precision weighting land in Phase 3; for the block design + ## rowMeans over treated units is the cell-weighted (per-treated-cell) center. + eff.tr.mat <- eff[, id.tr, drop = FALSE] + tr.path <- rowMeans(eff.tr.mat, na.rm = TRUE) # calendar + sc.tr <- .conformal_scale(tr.path[pre.idx], scale) + + w.co <- NULL # overlap density-ratio weighting deferred to a later phase + + ## --- scalar average effect: center over the post window, rank vs controls + m.co <- apply(G.co[, post.idx, drop = FALSE], 1L, + function(z) .conformal_center(z, center)) + keep <- is.finite(m.co) & is.finite(sc.co) & sc.co > 0 + m.tr <- .conformal_center(tr.path[post.idx], center) + ci <- .conformal_ci_level(m.tr, sc.tr, m.co[keep], sc.co[keep], + alpha = alpha, w.co = w.co) + Ncal <- sum(keep) + status <- if (any(is.infinite(ci))) "unbounded" else "ok" s.tr <- if (is.finite(sc.tr) && sc.tr > 0) abs(m.tr) / sc.tr else NA_real_ - s.co.vec <- abs(m.co) / sc.co + s.co.vec <- abs(m.co[keep]) / sc.co[keep] p.value <- .conformal_pval(s.tr, s.co.vec, w.co = w.co) + ## --- per-period pointwise band (calendar-indexed): at each period the treated + ## gap is ranked against the control gaps at the same period. Pre-period rows + ## give the placebo / pre-trend band; post-period rows give the effect band. + band <- matrix(NA_real_, TT, 4L, + dimnames = list(NULL, c("eff", "CI.lower", "CI.upper", "p.value"))) + for (t in seq_len(TT)) { + gco.t <- G.co[, t] + okt <- is.finite(gco.t) & is.finite(sc.co) & sc.co > 0 + if (!is.finite(tr.path[t]) || sum(okt) < 2L) next + ci.t <- .conformal_ci_level(tr.path[t], sc.tr, gco.t[okt], sc.co[okt], alpha = alpha) + s.tr.t <- if (is.finite(sc.tr) && sc.tr > 0) abs(tr.path[t]) / sc.tr else NA_real_ + band[t, ] <- c(tr.path[t], ci.t[1], ci.t[2], + .conformal_pval(s.tr.t, abs(gco.t[okt]) / sc.co[okt])) + } + list(att = m.tr, ci = c(ci[1], ci[2]), p.value = p.value, score.tr = s.tr, score.co = s.co.vec, status = status, n.calib = Ncal, n.eff = if (is.null(w.co)) Ncal else .conformal_neff(w.co), - scale = scale, center = center, weight = weight, form = "jackknife+") + scale = scale, center = center, weight = weight, form = "jackknife+", + band = band, post.idx = post.idx, pre.idx = pre.idx) } diff --git a/R/default.R b/R/default.R index 5e13c876..1cbdcdf4 100644 --- a/R/default.R +++ b/R/default.R @@ -58,6 +58,9 @@ fect <- function( se = FALSE, # report uncertainties vartype = "bootstrap", # bootstrap or jackknife para.error = "auto", # parametric bootstrap error strategy: "auto", "ar", "empirical", "wild" + conformal.scale = "none", # vartype="conformal": per-unit scale none|sd|rmspe|mad|diff|model-se + conformal.center = "mean", # vartype="conformal": post-period location mean|median + conformal.weight = "cell", # vartype="conformal": multi-treated aggregation cell|unit|precision cl = NULL, ci.method = "normal", # CI method for fect's est.* slots: "normal" (Wald: theta_hat +- z * SE) or "basic" (reflected pivot, Davison-Hinkley 1997 Sec. 5.2.1). For percentile / bc / bca on alternative estimands (att.cumu, aptt, log.att), call estimand(fit, type, ci.method) post-fit quantile.CI = NULL, # DEPRECATED: use ci.method instead. NULL sentinel = "not supplied"; legacy FALSE -> ci.method = "normal", legacy TRUE -> ci.method = "basic" @@ -141,6 +144,9 @@ fect.formula <- function( se = FALSE, # report uncertainties vartype = "bootstrap", # bootstrap or jackknife para.error = "auto", # parametric bootstrap error strategy: "auto", "ar", "empirical", "wild" + conformal.scale = "none", # vartype="conformal": per-unit scale none|sd|rmspe|mad|diff|model-se + conformal.center = "mean", # vartype="conformal": post-period location mean|median + conformal.weight = "cell", # vartype="conformal": multi-treated aggregation cell|unit|precision cl = NULL, ci.method = "normal", # CI method for fect's est.* slots: "normal" or "basic" quantile.CI = NULL, # DEPRECATED: use ci.method instead @@ -256,6 +262,9 @@ fect.formula <- function( se = se, vartype = vartype, para.error = para.error, + conformal.scale = conformal.scale, + conformal.center = conformal.center, + conformal.weight = conformal.weight, cl = cl, ci.method = ci.method, quantile.CI = quantile.CI, @@ -341,6 +350,9 @@ fect.default <- function( se = FALSE, # report uncertainties vartype = "bootstrap", # bootstrap or jackknife para.error = "auto", # parametric bootstrap error strategy: "auto", "ar", "empirical", "wild" + conformal.scale = "none", # vartype="conformal": per-unit scale none|sd|rmspe|mad|diff|model-se + conformal.center = "mean", # vartype="conformal": post-period location mean|median + conformal.weight = "cell", # vartype="conformal": multi-treated aggregation cell|unit|precision cl = NULL, ci.method = "normal", # CI method for fect's est.* slots: "normal" or "basic" quantile.CI = NULL, # DEPRECATED: use ci.method instead @@ -571,6 +583,19 @@ fect.default <- function( "The \"", vartype, "\" option is not available for the \"mc\" or \"both\" methods." ) } + if (vartype == "conformal") { + if (!conformal.scale %in% c("none", "sd", "rmspe", "mad", "diff", "model-se")) { + stop("conformal.scale must be one of \"none\", \"sd\", \"rmspe\", ", + "\"mad\", \"diff\", \"model-se\".", call. = FALSE) + } + if (!conformal.center %in% c("mean", "median")) { + stop("conformal.center must be \"mean\" or \"median\".", call. = FALSE) + } + if (!conformal.weight %in% c("cell", "unit", "precision")) { + stop("conformal.weight must be \"cell\", \"unit\", or \"precision\".", + call. = FALSE) + } + } if (vartype == "jackknife" && !is.null(cl)) { warning( "vartype = \"jackknife\" with cl = ... : the cl argument is ignored. ", @@ -2750,6 +2775,9 @@ fect.default <- function( carryover.period = carryover.period, vartype = vartype, para.error = para.error, + conformal.scale = conformal.scale, + conformal.center = conformal.center, + conformal.weight = conformal.weight, quantile.CI = .quantile.CI.bool, nboots = nboots, parallel = parallel, @@ -3349,7 +3377,10 @@ fect.default <- function( output$loo <- FALSE } - if (se == 1) { + ## conformal inference is rank-based with no bootstrap draws, so the + ## bootstrap F / equivalence diagnostics (diagtest) do not apply; skip them + ## (as with se = FALSE). The pre-period conformal band still shows pre-trends. + if (se == 1 && vartype != "conformal") { suppressWarnings( test.out <- diagtest( output, @@ -3366,7 +3397,7 @@ fect.default <- function( if (loo == TRUE) { output$loo <- TRUE } - if (loo == TRUE && se == 1) { + if (loo == TRUE && se == 1 && vartype != "conformal") { suppressWarnings( test.out <- diagtest( output, diff --git a/tests/testthat/test-conformal.R b/tests/testthat/test-conformal.R index a516c4cf..0d9ef189 100644 --- a/tests/testthat/test-conformal.R +++ b/tests/testthat/test-conformal.R @@ -108,3 +108,43 @@ test_that("studentized-mean (scale = sd) runs, never empty/unbounded, ballpark c expect_equal(r$bad, 0L) expect_gt(r$cov, 0.75); expect_lte(r$cov, 1.0) }) + +## ---- fect(vartype = "conformal") end-to-end integration --------------------- + +test_that("fect(vartype = 'conformal') returns a complete, printable object", { + skip_on_cran() + g <- .conf_dgp(1) + Y <- g$mu + matrix(rnorm(g$N * g$T), g$T, g$N) + 1.5 * g$D # true ATT = 1.5 + dat <- data.frame(id = rep(1:g$N, each = g$T), time = rep(1:g$T, g$N), + Y = as.vector(Y), D = as.vector(g$D)) + f <- suppressMessages(fect(Y ~ D, data = dat, index = c("id", "time"), + method = "gsynth", force = 3, CV = FALSE, r = 2, se = TRUE, + vartype = "conformal")) + expect_equal(f$vartype, "conformal") + expect_identical(f$conformal$status, "ok") + ## est.avg: finite, ordered CI with the right columns + expect_true(all(c("ATT.avg", "CI.lower", "CI.upper", "p.value") %in% colnames(f$est.avg))) + expect_lt(f$est.avg[1, "CI.lower"], f$est.avg[1, "CI.upper"]) + ## est.att: event-time rownames, post-period CIs finite + expect_equal(rownames(f$est.att), as.character(f$time)) + post <- f$est.att[as.numeric(rownames(f$est.att)) >= 0, , drop = FALSE] + expect_true(all(is.finite(post[, c("CI.lower", "CI.upper")]))) + expect_silent(invisible(capture.output(print(f)))) +}) + +test_that("conformal.scale threads through fect() and changes the interval", { + skip_on_cran() + g <- .conf_dgp(7) + Y <- g$mu + matrix(rnorm(g$N * g$T), g$T, g$N) + 1.5 * g$D + dat <- data.frame(id = rep(1:g$N, each = g$T), time = rep(1:g$T, g$N), + Y = as.vector(Y), D = as.vector(g$D)) + run <- function(scale) suppressMessages(fect(Y ~ D, data = dat, index = c("id", "time"), + method = "gsynth", force = 3, CV = FALSE, r = 2, se = TRUE, + vartype = "conformal", conformal.scale = scale)) + w_none <- diff(run("none")$est.avg[1, c("CI.lower", "CI.upper")]) + w_diff <- diff(run("diff")$est.avg[1, c("CI.lower", "CI.upper")]) + expect_true(is.finite(w_none) && is.finite(w_diff)) + expect_false(isTRUE(all.equal(w_none, w_diff))) # the scale actually does something + ## not-yet-implemented knobs error clearly + expect_error(run("model-se"), "not yet implemented") +}) From 6e4a4691076183e1563df1872a2372f831dc3a64 Mon Sep 17 00:00:00 2001 From: Yiqing Xu Date: Tue, 9 Jun 2026 00:17:20 -0700 Subject: [PATCH 06/14] feat(conformal): multi-treated weights cell/unit/precision (Phase 3) Complete the knob set. Weights aggregate across treated units; model-se scale is deferred (needs estimator prediction-SE plumbing) and stops with a clear message. ## conformal_calibrate - Scalar ATT = weighted mean of per-unit post centers; per-period band uses the matching weighted treated trajectory. cell = per-treated-cell (n_i), unit = equal per unit, precision = inverse pre-period variance. Precision is decoupled from the `scale` knob so it down-weights noisy units even under the default scale = "none". ## boot.R - Report the conformal weighted center (cc$att / band eff) as the point estimates, so the point and the symmetric CI share one estimand. For weight = "cell" this equals fect's canonical att.avg (verified to 1e-6). ## Scope - scale: none/sd/rmspe/mad/diff work; model-se deferred (clear stop()). - center: mean/median; per-horizon is the band (always produced). - weight: cell/unit/precision. ## Verification - 48 tests pass: added multi-treated weight + precision-down-weighting tests. - Existing vartypes untouched (all changes gated to the conformal path). Refs statsclaw runs/2026-06-08-conformal-inference.md (overnight build, Phase 3). Co-Authored-By: Claude Opus 4.8 --- R/boot.R | 18 ++++++++----- R/conformal.R | 48 ++++++++++++++++++++++++--------- tests/testthat/test-conformal.R | 39 +++++++++++++++++++++++++++ 3 files changed, 85 insertions(+), 20 deletions(-) diff --git a/R/boot.R b/R/boot.R index ebb7e935..57fc5cf7 100644 --- a/R/boot.R +++ b/R/boot.R @@ -582,18 +582,21 @@ fect_boot <- function( z <- stats::qnorm(1 - alpha / 2) se.from <- function(lo, hi) ifelse(is.finite(hi - lo), (hi - lo) / (2 * z), NA_real_) - ## --- average effect -> est.avg / est.avg.unit + ## --- average effect -> est.avg / est.avg.unit. Point = the conformal + ## weighted center cc$att (so the symmetric CI is centered on the reported + ## estimate); for the default weight = "cell" this equals fect's att.avg. + ## conformal reports a single weighting, so both print rows show it. se.avg <- se.from(cc$ci[1], cc$ci[2]) - est.avg <- t(as.matrix(c(att.avg, se.avg, cc$ci[1], cc$ci[2], cc$p.value))) + est.avg <- t(as.matrix(c(cc$att, se.avg, cc$ci[1], cc$ci[2], cc$p.value))) colnames(est.avg) <- c("ATT.avg", "S.E.", "CI.lower", "CI.upper", "p.value") - att.avg.unit.val <- if (length(att.avg.unit)) att.avg.unit else att.avg - est.avg.unit <- t(as.matrix(c(att.avg.unit.val, se.avg, cc$ci[1], cc$ci[2], cc$p.value))) + est.avg.unit <- t(as.matrix(c(cc$att, se.avg, cc$ci[1], cc$ci[2], cc$p.value))) colnames(est.avg.unit) <- c("ATT.avg.unit", "S.E.", "CI.lower", "CI.upper", "p.value") - ## --- per-period band (calendar-indexed) -> est.eff.calendar(.fit) directly + ## --- per-period band (calendar-indexed) -> est.eff.calendar(.fit). The + ## ATT-calendar point is the conformal weighted per-period effect (band eff). band <- cc$band se.cal <- se.from(band[, "CI.lower"], band[, "CI.upper"]) - est.eff.calendar <- cbind(calendar.eff, se.cal, + est.eff.calendar <- cbind(band[, "eff"], se.cal, band[, "CI.lower"], band[, "CI.upper"], band[, "p.value"], calendar.N) colnames(est.eff.calendar) <- c("ATT-calendar", "S.E.", "CI.lower", "CI.upper", "p.value", "count") est.eff.calendar.fit <- cbind(calendar.eff.fit, se.cal, @@ -609,10 +612,11 @@ fect_boot <- function( for (k in seq_along(time.on)) { rows <- which(rel == time.on[k]) if (length(rows) == 0L) next + ef <- mean(band[rows, "eff"], na.rm = TRUE) lo <- mean(band[rows, "CI.lower"], na.rm = TRUE) hi <- mean(band[rows, "CI.upper"], na.rm = TRUE) pv <- mean(band[rows, "p.value"], na.rm = TRUE) - est.att[k, ] <- c(att[k], se.from(lo, hi), lo, hi, pv, out$count[k]) + est.att[k, ] <- c(ef, se.from(lo, hi), lo, hi, pv, out$count[k]) } ## Phase 4 builds a true (1 - 2*alpha) inner band and the simultaneous band; ## for now the 90%-slot mirrors the main band so plot paths do not break. diff --git a/R/conformal.R b/R/conformal.R index ebfce0c8..9c0abcb0 100644 --- a/R/conformal.R +++ b/R/conformal.R @@ -128,15 +128,11 @@ conformal_calibrate <- function(Y, D, X = NULL, I, II, T.on, r.cv, eff, Ntr <- length(id.tr); Nco <- length(id.co) if (Nco < 2L) stop("conformal: need at least 2 controls.") - ## not-yet-implemented knobs (land in a later phase): fail loudly, not silently. + ## model-se needs the estimator's per-unit prediction SE (not wired yet). if (identical(scale, "model-se")) { stop("conformal: scale = \"model-se\" is not yet implemented; use one of ", "\"none\", \"sd\", \"rmspe\", \"mad\", \"diff\".", call. = FALSE) } - if (!identical(weight, "cell")) { - stop("conformal: weight = \"", weight, "\" is not yet implemented; ", - "use weight = \"cell\" (the per-treated-cell ATT).", call. = FALSE) - } ## --- separation guard: conformal calibration is controls-only (nevertreated) if (!identical(predictive, "nevertreated")) { @@ -212,20 +208,46 @@ conformal_calibrate <- function(Y, D, X = NULL, I, II, T.on, r.cv, eff, sc.co[mi] <- .conformal_scale(gj[pre.idx], scale) } - ## --- treated aggregate path (cell-weighted average over treated units). - ## NOTE: unit / precision weighting land in Phase 3; for the block design - ## rowMeans over treated units is the cell-weighted (per-treated-cell) center. - eff.tr.mat <- eff[, id.tr, drop = FALSE] - tr.path <- rowMeans(eff.tr.mat, na.rm = TRUE) # calendar - sc.tr <- .conformal_scale(tr.path[pre.idx], scale) + ## --- treated aggregate across treated units, combined per `weight`: + ## cell = weight each unit by its post-obs count (per-treated-cell ATT), + ## unit = each treated unit counts equally, + ## precision = inverse pre-period variance (down-weights noisy units). + ## The scalar center is the weighted mean of the per-unit post centers; the + ## per-period band uses the matching weighted treated trajectory. They coincide + ## for a single treated unit and (mean center) for a balanced panel. + ## Multi-treated NOTE: the aggregate is ranked against single-control gaps; a + ## placebo-AVERAGE calibration (ranking against averages of Ntr controls) is a + ## refinement deferred to the simulation phase. + trc <- eff[, id.tr, drop = FALSE] # TT x Ntr + m.i <- apply(trc[post.idx, , drop = FALSE], 2L, function(z) .conformal_center(z, center)) + n.i <- colSums(!is.na(trc[post.idx, , drop = FALSE])) + ## precision weight uses each unit's pre-period VARIANCE directly (decoupled + ## from the `scale` knob, which only normalizes the score), so it down-weights + ## noisy units even under the default scale = "none". + var.i <- apply(trc[pre.idx, , drop = FALSE], 2L, function(z) { + z <- z[is.finite(z)]; if (length(z) > 1L) stats::var(z) else NA_real_ + }) + wi <- switch(weight, + cell = n.i, + unit = rep(1, Ntr), + precision = ifelse(is.finite(var.i) & var.i > 0, 1 / var.i, 0)) + wi[!is.finite(wi)] <- 0 + if (!any(wi > 0)) wi <- rep(1, Ntr) + Wt <- matrix(wi, TT, Ntr, byrow = TRUE); Wt[is.na(trc)] <- 0 + denom <- rowSums(Wt) + tr.path <- ifelse(denom > 0, rowSums(trc * Wt, na.rm = TRUE) / denom, NA_real_) # calendar + sc.tr <- .conformal_scale(tr.path[pre.idx], scale) w.co <- NULL # overlap density-ratio weighting deferred to a later phase - ## --- scalar average effect: center over the post window, rank vs controls + ## --- scalar average effect: weighted mean of per-unit post centers, ranked vs + ## the per-control post centers. m.co <- apply(G.co[, post.idx, drop = FALSE], 1L, function(z) .conformal_center(z, center)) keep <- is.finite(m.co) & is.finite(sc.co) & sc.co > 0 - m.tr <- .conformal_center(tr.path[post.idx], center) + okm <- is.finite(m.i) & wi > 0 + m.tr <- if (any(okm)) sum(wi[okm] * m.i[okm]) / sum(wi[okm]) else + .conformal_center(tr.path[post.idx], center) ci <- .conformal_ci_level(m.tr, sc.tr, m.co[keep], sc.co[keep], alpha = alpha, w.co = w.co) Ncal <- sum(keep) diff --git a/tests/testthat/test-conformal.R b/tests/testthat/test-conformal.R index 0d9ef189..6363934a 100644 --- a/tests/testthat/test-conformal.R +++ b/tests/testthat/test-conformal.R @@ -148,3 +148,42 @@ test_that("conformal.scale threads through fect() and changes the interval", { ## not-yet-implemented knobs error clearly expect_error(run("model-se"), "not yet implemented") }) + +## ---- multi-treated weighting (cell / unit / precision) ---------------------- + +test_that("conformal weights run for multiple treated units; cell matches att.avg", { + skip_on_cran() + set.seed(3); N <- 30; T <- 20; T0 <- 15; r <- 2 + Fm <- matrix(rnorm(T * r), T, r); L <- matrix(rnorm(N * r), N, r) + D <- matrix(0, T, N); D[(T0 + 1):T, 1:4] <- 1 + Y <- Fm %*% t(L) + matrix(rnorm(N * T), T, N) + 1.5 * D + dat <- data.frame(id = rep(1:N, each = T), time = rep(1:T, N), + Y = as.vector(Y), D = as.vector(D)) + pf <- suppressMessages(fect(Y ~ D, data = dat, index = c("id", "time"), + method = "gsynth", force = 3, CV = FALSE, r = 2, se = FALSE)) + for (w in c("cell", "unit", "precision")) { + f <- suppressMessages(fect(Y ~ D, data = dat, index = c("id", "time"), + method = "gsynth", force = 3, CV = FALSE, r = 2, se = TRUE, + vartype = "conformal", conformal.weight = w)) + expect_identical(f$conformal$status, "ok") + expect_true(all(is.finite(f$est.avg[1, c("CI.lower", "CI.upper")]))) + if (w == "cell") { + expect_equal(unname(f$est.avg[1, "ATT.avg"]), pf$att.avg, tolerance = 1e-6) + } + } +}) + +test_that("precision weight down-weights a noisy treated unit", { + skip_on_cran() + set.seed(5); N <- 30; T <- 20; T0 <- 15; r <- 2 + Fm <- matrix(rnorm(T * r), T, r); L <- matrix(rnorm(N * r), N, r) + D <- matrix(0, T, N); D[(T0 + 1):T, 1:4] <- 1 + E <- matrix(rnorm(N * T), T, N); E[, 1] <- E[, 1] * 4 # treated unit 1 noisy + Y <- Fm %*% t(L) + E + 1.5 * D + dat <- data.frame(id = rep(1:N, each = T), time = rep(1:T, N), + Y = as.vector(Y), D = as.vector(D)) + att <- function(w) suppressMessages(fect(Y ~ D, data = dat, index = c("id", "time"), + method = "gsynth", force = 3, CV = FALSE, r = 2, se = TRUE, + vartype = "conformal", conformal.weight = w))$est.avg[1, "ATT.avg"] + expect_false(isTRUE(all.equal(att("unit"), att("precision")))) +}) From fc2a415df6dd8ea32a26f560222b079310d5cdf0 Mon Sep 17 00:00:00 2001 From: Yiqing Xu Date: Tue, 9 Jun 2026 00:29:11 -0700 Subject: [PATCH 07/14] feat(conformal): pointwise / simultaneous bands + pooled cutoff (Phase 4) ## Bands - conformal.band: "pointwise" (default) | "simultaneous"; conformal.cutoff: "per-period" (default) | "pooled". Threaded through all signatures + validated. - pointwise per-period: treated gap ranked vs control gaps at that period. - pooled: one cutoff over all post-window standardized control gaps. - simultaneous: one multiplier c = conformal quantile of the per-control MAX standardized gap over the post window (sup-t uniform band); pre rows keep the per-period placebo band. - est.att carries the selected band; est.att.sim (uniform band, event-time mapped) is always stored. conformal meta gains band.type / cutoff. ## Coverage (150 reps, block, effect 0, alpha=0.10) - pointwise per-period 0.879; pointwise JOINT 0.440 (multiple-comparison failure). - simultaneous JOINT 0.773: materially better jointly, but undercovers nominal in small N (near the jackknife+ floor 1-2*alpha + the treated-vs-control LOO asymmetry). A symmetric treated-gap fix is flagged for the simulation phase. ## Verification - 55 tests pass (added band-width, est.att.sim, pooled, joint-coverage tests). Loopy tests use parallel = FALSE (avoids PSOCK cluster exhaustion). - bootstrap / jackknife unaffected. Refs statsclaw runs/2026-06-08-conformal-inference.md (overnight build, Phase 4). Co-Authored-By: Claude Opus 4.8 --- R/boot.R | 44 +++++++++++++++---------- R/conformal.R | 58 ++++++++++++++++++++++++--------- R/default.R | 18 ++++++++++ tests/testthat/test-conformal.R | 58 ++++++++++++++++++++++++++++++++- 4 files changed, 145 insertions(+), 33 deletions(-) diff --git a/R/boot.R b/R/boot.R index 57fc5cf7..fb428570 100644 --- a/R/boot.R +++ b/R/boot.R @@ -145,6 +145,8 @@ fect_boot <- function( conformal.scale = "none", conformal.center = "mean", conformal.weight = "cell", + conformal.band = "pointwise", + conformal.cutoff = "per-period", quantile.CI = FALSE, nboots = 200, parallel = TRUE, @@ -574,7 +576,8 @@ fect_boot <- function( force = force, hasRevs = hasRevs, tol = tol, max.iteration = max.iteration, norm.para = norm.para, scale = conformal.scale, center = conformal.center, - weight = conformal.weight, alpha = alpha + weight = conformal.weight, band.type = conformal.band, + cutoff = conformal.cutoff, alpha = alpha ) ## back out a nominal S.E. from the symmetric conformal CI (display only; @@ -603,23 +606,28 @@ fect_boot <- function( band[, "CI.lower"], band[, "CI.upper"], band[, "p.value"], calendar.N) colnames(est.eff.calendar.fit) <- colnames(est.eff.calendar) - ## --- map the calendar band to event time (treated unit's relative period) - ## for est.att (rownames = out$time). Common-onset: a 1:1 map. + ## --- map a calendar band to event time (treated unit's relative period); + ## rownames = out$time. Common-onset: a 1:1 map. rel <- T.on[, id.tr.conf[1]] - est.att <- matrix(NA_real_, length(time.on), 6, - dimnames = list(time.on, - c("ATT", "S.E.", "CI.lower", "CI.upper", "p.value", "count"))) - for (k in seq_along(time.on)) { - rows <- which(rel == time.on[k]) - if (length(rows) == 0L) next - ef <- mean(band[rows, "eff"], na.rm = TRUE) - lo <- mean(band[rows, "CI.lower"], na.rm = TRUE) - hi <- mean(band[rows, "CI.upper"], na.rm = TRUE) - pv <- mean(band[rows, "p.value"], na.rm = TRUE) - est.att[k, ] <- c(ef, se.from(lo, hi), lo, hi, pv, out$count[k]) + map_band <- function(bnd) { + m <- matrix(NA_real_, length(time.on), 6, + dimnames = list(time.on, + c("ATT", "S.E.", "CI.lower", "CI.upper", "p.value", "count"))) + for (k in seq_along(time.on)) { + rows <- which(rel == time.on[k]) + if (length(rows) == 0L) next + ef <- mean(bnd[rows, "eff"], na.rm = TRUE) + lo <- mean(bnd[rows, "CI.lower"], na.rm = TRUE) + hi <- mean(bnd[rows, "CI.upper"], na.rm = TRUE) + pv <- mean(bnd[rows, "p.value"], na.rm = TRUE) + m[k, ] <- c(ef, se.from(lo, hi), lo, hi, pv, out$count[k]) + } + m } - ## Phase 4 builds a true (1 - 2*alpha) inner band and the simultaneous band; - ## for now the 90%-slot mirrors the main band so plot paths do not break. + est.att <- map_band(band) # the selected band (conformal.band) + est.att.sim <- map_band(cc$band.sim) # the uniform / simultaneous band + ## Phase 4 note: est.att90 (the inner equivalence band) mirrors the main band + ## for now; a true (1 - 2*alpha) conformal inner band is a later refinement. est.att90 <- est.att att.bound <- est.att[, c("CI.lower", "CI.upper"), drop = FALSE] rownames(att.bound) <- time.on @@ -628,12 +636,14 @@ fect_boot <- function( est.avg = est.avg, est.avg.unit = est.avg.unit, est.att = est.att, + est.att.sim = est.att.sim, est.att90 = est.att90, att.bound = att.bound, est.eff.calendar = est.eff.calendar, est.eff.calendar.fit = est.eff.calendar.fit, vartype = "conformal", - conformal = cc[c("scale", "center", "weight", "status", "n.calib", "form")] + conformal = cc[c("scale", "center", "weight", "band.type", "cutoff", + "status", "n.calib", "form")] ) return(c(out, result)) } diff --git a/R/conformal.R b/R/conformal.R index 9c0abcb0..975fc088 100644 --- a/R/conformal.R +++ b/R/conformal.R @@ -120,6 +120,7 @@ conformal_calibrate <- function(Y, D, X = NULL, I, II, T.on, r.cv, eff, force = 3L, hasRevs = 0L, tol = 1e-5, max.iteration = 1000L, norm.para = NULL, scale = "none", center = "mean", weight = "cell", + band.type = "pointwise", cutoff = "per-period", alpha = 0.05, conformal.fit = NULL) { TT <- nrow(Y); N <- ncol(Y) @@ -256,24 +257,51 @@ conformal_calibrate <- function(Y, D, X = NULL, I, II, T.on, r.cv, eff, s.co.vec <- abs(m.co[keep]) / sc.co[keep] p.value <- .conformal_pval(s.tr, s.co.vec, w.co = w.co) - ## --- per-period pointwise band (calendar-indexed): at each period the treated - ## gap is ranked against the control gaps at the same period. Pre-period rows - ## give the placebo / pre-trend band; post-period rows give the effect band. - band <- matrix(NA_real_, TT, 4L, - dimnames = list(NULL, c("eff", "CI.lower", "CI.upper", "p.value"))) - for (t in seq_len(TT)) { - gco.t <- G.co[, t] - okt <- is.finite(gco.t) & is.finite(sc.co) & sc.co > 0 - if (!is.finite(tr.path[t]) || sum(okt) < 2L) next - ci.t <- .conformal_ci_level(tr.path[t], sc.tr, gco.t[okt], sc.co[okt], alpha = alpha) - s.tr.t <- if (is.finite(sc.tr) && sc.tr > 0) abs(tr.path[t]) / sc.tr else NA_real_ - band[t, ] <- c(tr.path[t], ci.t[1], ci.t[2], - .conformal_pval(s.tr.t, abs(gco.t[okt]) / sc.co[okt])) + ## --- per-period band(s) (calendar-indexed). Standardized control gaps + ## g.tilde[j,t] = |G.co[j,t]| / sc.co[j]; the treated deviation at period t is + ## |tr.path[t] - tau| / sc.tr. Three constructions: + ## pointwise + per-period (default): Q_t = conformal quantile of g.tilde[,t]. + ## pointwise + pooled: one Q over all post-window g.tilde values. + ## simultaneous: one multiplier c = quantile of the per-control MAX over the + ## post window (sup-t / uniform band over the post path); pre-period rows + ## keep the per-period placebo band. + ## The pre-period rows always show the per-period placebo band (for pre-trends). + ok.co <- which(is.finite(sc.co) & sc.co > 0) + Gv <- G.co[ok.co, , drop = FALSE] + scv <- sc.co[ok.co] + Gtil <- abs(Gv) / scv # ncal x TT, row j divided by scv[j] + + Qpool <- NULL + if (identical(cutoff, "pooled")) { + pv <- as.vector(Gtil[, post.idx, drop = FALSE]); pv <- pv[is.finite(pv)] + Qpool <- .conformal_quantile(pv, alpha) + } + Mj <- apply(Gtil[, post.idx, drop = FALSE], 1L, + function(z) { z <- z[is.finite(z)]; if (length(z)) max(z) else NA_real_ }) + c.sim <- .conformal_quantile(Mj[is.finite(Mj)], alpha) + + mk_band <- function(mode) { + b <- matrix(NA_real_, TT, 4L, + dimnames = list(NULL, c("eff", "CI.lower", "CI.upper", "p.value"))) + for (t in seq_len(TT)) { + gt <- Gtil[, t]; okt <- is.finite(gt) + if (!is.finite(tr.path[t]) || !is.finite(sc.tr) || sc.tr <= 0 || sum(okt) < 2L) next + Qt <- if (mode == "simultaneous" && t %in% post.idx) c.sim + else if (identical(cutoff, "pooled")) Qpool + else .conformal_quantile(gt[okt], alpha) + half <- if (is.finite(Qt)) sc.tr * Qt else Inf + b[t, ] <- c(tr.path[t], tr.path[t] - half, tr.path[t] + half, + .conformal_pval(abs(tr.path[t]) / sc.tr, gt[okt])) + } + b } + band <- mk_band(band.type) # populates est.att + band.sim <- if (identical(band.type, "simultaneous")) band else mk_band("simultaneous") list(att = m.tr, ci = c(ci[1], ci[2]), p.value = p.value, score.tr = s.tr, score.co = s.co.vec, status = status, n.calib = Ncal, n.eff = if (is.null(w.co)) Ncal else .conformal_neff(w.co), - scale = scale, center = center, weight = weight, form = "jackknife+", - band = band, post.idx = post.idx, pre.idx = pre.idx) + scale = scale, center = center, weight = weight, + band.type = band.type, cutoff = cutoff, form = "jackknife+", + band = band, band.sim = band.sim, post.idx = post.idx, pre.idx = pre.idx) } diff --git a/R/default.R b/R/default.R index 1cbdcdf4..bf26116e 100644 --- a/R/default.R +++ b/R/default.R @@ -61,6 +61,8 @@ fect <- function( conformal.scale = "none", # vartype="conformal": per-unit scale none|sd|rmspe|mad|diff|model-se conformal.center = "mean", # vartype="conformal": post-period location mean|median conformal.weight = "cell", # vartype="conformal": multi-treated aggregation cell|unit|precision + conformal.band = "pointwise", # vartype="conformal": est.att band pointwise|simultaneous + conformal.cutoff = "per-period", # vartype="conformal": pointwise cutoff per-period|pooled cl = NULL, ci.method = "normal", # CI method for fect's est.* slots: "normal" (Wald: theta_hat +- z * SE) or "basic" (reflected pivot, Davison-Hinkley 1997 Sec. 5.2.1). For percentile / bc / bca on alternative estimands (att.cumu, aptt, log.att), call estimand(fit, type, ci.method) post-fit quantile.CI = NULL, # DEPRECATED: use ci.method instead. NULL sentinel = "not supplied"; legacy FALSE -> ci.method = "normal", legacy TRUE -> ci.method = "basic" @@ -147,6 +149,8 @@ fect.formula <- function( conformal.scale = "none", # vartype="conformal": per-unit scale none|sd|rmspe|mad|diff|model-se conformal.center = "mean", # vartype="conformal": post-period location mean|median conformal.weight = "cell", # vartype="conformal": multi-treated aggregation cell|unit|precision + conformal.band = "pointwise", # vartype="conformal": est.att band pointwise|simultaneous + conformal.cutoff = "per-period", # vartype="conformal": pointwise cutoff per-period|pooled cl = NULL, ci.method = "normal", # CI method for fect's est.* slots: "normal" or "basic" quantile.CI = NULL, # DEPRECATED: use ci.method instead @@ -265,6 +269,8 @@ fect.formula <- function( conformal.scale = conformal.scale, conformal.center = conformal.center, conformal.weight = conformal.weight, + conformal.band = conformal.band, + conformal.cutoff = conformal.cutoff, cl = cl, ci.method = ci.method, quantile.CI = quantile.CI, @@ -353,6 +359,8 @@ fect.default <- function( conformal.scale = "none", # vartype="conformal": per-unit scale none|sd|rmspe|mad|diff|model-se conformal.center = "mean", # vartype="conformal": post-period location mean|median conformal.weight = "cell", # vartype="conformal": multi-treated aggregation cell|unit|precision + conformal.band = "pointwise", # vartype="conformal": est.att band pointwise|simultaneous + conformal.cutoff = "per-period", # vartype="conformal": pointwise cutoff per-period|pooled cl = NULL, ci.method = "normal", # CI method for fect's est.* slots: "normal" or "basic" quantile.CI = NULL, # DEPRECATED: use ci.method instead @@ -595,6 +603,14 @@ fect.default <- function( stop("conformal.weight must be \"cell\", \"unit\", or \"precision\".", call. = FALSE) } + if (!conformal.band %in% c("pointwise", "simultaneous")) { + stop("conformal.band must be \"pointwise\" or \"simultaneous\".", + call. = FALSE) + } + if (!conformal.cutoff %in% c("per-period", "pooled")) { + stop("conformal.cutoff must be \"per-period\" or \"pooled\".", + call. = FALSE) + } } if (vartype == "jackknife" && !is.null(cl)) { warning( @@ -2778,6 +2794,8 @@ fect.default <- function( conformal.scale = conformal.scale, conformal.center = conformal.center, conformal.weight = conformal.weight, + conformal.band = conformal.band, + conformal.cutoff = conformal.cutoff, quantile.CI = .quantile.CI.bool, nboots = nboots, parallel = parallel, diff --git a/tests/testthat/test-conformal.R b/tests/testthat/test-conformal.R index 6363934a..89d1d4b4 100644 --- a/tests/testthat/test-conformal.R +++ b/tests/testthat/test-conformal.R @@ -78,7 +78,8 @@ test_that("effective sample size is sensible", { dat <- data.frame(id = rep(1:g$N, each = g$T), time = rep(1:g$T, g$N), Y = as.vector(Y), D = as.vector(g$D)) fit <- suppressMessages(fect(Y ~ D, data = dat, index = c("id", "time"), - method = "gsynth", force = 3, CV = FALSE, r = 2, se = FALSE)) + method = "gsynth", force = 3, CV = FALSE, r = 2, se = FALSE, + parallel = FALSE)) cc <- conformal_calibrate(Y = Y, D = g$D, I = fit$I, II = fit$II, T.on = fit$T.on, r.cv = fit$r.cv, eff = fit$eff, method = "gsynth", scale = scale, alpha = alpha) @@ -187,3 +188,58 @@ test_that("precision weight down-weights a noisy treated unit", { vartype = "conformal", conformal.weight = w))$est.avg[1, "ATT.avg"] expect_false(isTRUE(all.equal(att("unit"), att("precision")))) }) + +## ---- bands: pointwise / simultaneous / pooled cutoff ------------------------ + +test_that("simultaneous band is wider than pointwise; est.att.sim always present", { + skip_on_cran() + g <- .conf_dgp(1) + Y <- g$mu + matrix(rnorm(g$N * g$T), g$T, g$N) + 1.5 * g$D + dat <- data.frame(id = rep(1:g$N, each = g$T), time = rep(1:g$T, g$N), + Y = as.vector(Y), D = as.vector(g$D)) + run <- function(...) suppressMessages(fect(Y ~ D, data = dat, index = c("id", "time"), + method = "gsynth", force = 3, CV = FALSE, r = 2, se = TRUE, + vartype = "conformal", ...)) + fp <- run(); fs <- run(conformal.band = "simultaneous") + post_w <- function(m) { + pm <- m[as.numeric(rownames(m)) >= 0, c("CI.lower", "CI.upper"), drop = FALSE] + mean(pm[, 2] - pm[, 1]) + } + expect_gte(post_w(fs$est.att), post_w(fp$est.att)) # uniform band is wider + expect_false(is.null(fp$est.att.sim)) # always computed + expect_gte(post_w(fp$est.att.sim), post_w(fp$est.att)) + ## pooled cutoff runs and is finite + fpool <- run(conformal.cutoff = "pooled") + expect_identical(fpool$conformal$status, "ok") + expect_true(all(is.finite(fpool$est.att[as.numeric(rownames(fpool$est.att)) >= 0, + c("CI.lower", "CI.upper")]))) +}) + +test_that("simultaneous band gives materially better joint coverage than pointwise", { + skip_on_cran() + ## Joint (whole-post-path) coverage. Pointwise undercovers jointly (the + ## multiple-comparison problem, ~0.44 at 150 reps); the simultaneous band is + ## much better (~0.77). NOTE: it undercovers nominal 1 - alpha in this small-N + ## regime, sitting near the jackknife+ floor 1 - 2*alpha plus the treated-vs- + ## control LOO asymmetry; the precise characterization (and a possible + ## symmetric-treated-gap fix) is a Phase 5 item. Here we assert the robust + ## qualitative property only. + g <- .conf_dgp(909); reps <- 40; alpha <- 0.10 + joint <- function(band) { + post <- band[as.numeric(rownames(band)) >= 0, c("CI.lower", "CI.upper"), drop = FALSE] + all(post[, 1] <= 0 & 0 <= post[, 2]) + } + jp <- js <- logical(reps) + for (b in seq_len(reps)) { + Y <- g$mu + matrix(rnorm(g$N * g$T), g$T, g$N) # true effect 0 everywhere + dat <- data.frame(id = rep(1:g$N, each = g$T), time = rep(1:g$T, g$N), + Y = as.vector(Y), D = as.vector(g$D)) + f <- suppressMessages(fect(Y ~ D, data = dat, index = c("id", "time"), + method = "gsynth", force = 3, CV = FALSE, r = 2, se = TRUE, + vartype = "conformal", alpha = alpha, parallel = FALSE)) + jp[b] <- joint(f$est.att) # pointwise + js[b] <- joint(f$est.att.sim) # simultaneous + } + expect_gt(mean(js) - mean(jp), 0.15) # simultaneous materially better jointly + expect_gt(mean(js), 0.65) # and well above pointwise's joint coverage +}) From 61401e7f175ad94537da986e56983c19f4f1e5d5 Mon Sep 17 00:00:00 2001 From: Yiqing Xu Date: Tue, 9 Jun 2026 00:38:09 -0700 Subject: [PATCH 08/14] feat(conformal): default scale = "sd" from simulation; docs + NEWS (Phase 5/6) ## Default chosen by simulation - tests/coverage-study/conformal_default_study.R sweeps scale x DGP (iid/ar1/hetero/nonstat, 150 reps, true effect 0) + a weight sweep on a staggered+hetero DGP. Results in results/conformal_default_study.md. - scale = "sd" wins: highest worst-case scalar coverage (0.920) and smallest mean width (3.16). Decisive under heteroskedasticity: none over-covers at width 5.43 vs sd 0.920 at 2.82. Set as the provisional default across all signatures (flagged for the user's final call). - Weight finding: under treated-unit heterogeneity cell/unit undercover (0.873) and precision recovers (0.993); cell stays the default estimand. - bad = 0 throughout (never empty/unbounded, as designed). ## Docs - man/fect.Rd: document vartype = "conformal" + the five conformal.* args; usage block updated. - NEWS.md: 2.4.6 entry. DESCRIPTION: 2.4.5 -> 2.4.6, date 2026-06-09. ## Verification - 55 conformal tests pass with the new sd default. Refs statsclaw runs/2026-06-08-conformal-inference.md (overnight build, Phase 5/6 docs). Co-Authored-By: Claude Opus 4.8 --- DESCRIPTION | 4 +- NEWS.md | 6 + R/boot.R | 2 +- R/default.R | 6 +- man/fect.Rd | 44 ++++- .../coverage-study/conformal_default_study.R | 159 ++++++++++++++++++ .../results/conformal_default_study.md | 47 ++++++ 7 files changed, 257 insertions(+), 11 deletions(-) create mode 100644 tests/coverage-study/conformal_default_study.R create mode 100644 tests/coverage-study/results/conformal_default_study.md diff --git a/DESCRIPTION b/DESCRIPTION index 62405a1c..6e96ba59 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: fect Type: Package Title: Fixed Effects Counterfactual Estimators -Version: 2.4.5 -Date: 2026-05-29 +Version: 2.4.6 +Date: 2026-06-09 Authors@R: c(person("Yiqing", "Xu", , "yiqingxu@stanford.edu", role = c("aut", "cre")), person("Licheng", "Liu", , "lichengl@stanford.edu", role = c("aut")), diff --git a/NEWS.md b/NEWS.md index e1f53b11..d3903397 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,10 @@ +# fect 2.4.6 + +* New `vartype = "conformal"`: a cross-sectional conformal prediction interval for the average treatment effect on the treated. It ranks the treated unit's average post-treatment prediction error against the leave-one-control-out errors of the donors, giving a distribution-free interval with no bootstrap draws. The interval is closed-form and never empty. +* Conformal options: `conformal.scale` (`"sd"` default, the studentized statistic; also `"none"`, `"rmspe"`, `"mad"`, `"diff"`), `conformal.center` (`"mean"`/`"median"`), `conformal.weight` (`"cell"`/`"unit"`/`"precision"` for multiple treated units), `conformal.band` (`"pointwise"`/`"simultaneous"`), and `conformal.cutoff` (`"per-period"`/`"pooled"`). The simultaneous (uniform) band is also stored in `fit$est.att.sim`. +* `vartype = "conformal"` requires a separated (controls-only) fit and `method` not in `c("mc", "both")`. It currently covers the core synthetic-control case; group, reversal, weighted (`W`), balanced-panel, placebo, and carryover options still use `bootstrap`/`jackknife`. + # fect 2.4.5 * Add `group.fe` to `fect()` for absorbing coarser fixed effects, such as state FE with county-level data. Closes #139. Clustered SE defaults to `group.fe[1]`; override with `cl = ""`. diff --git a/R/boot.R b/R/boot.R index fb428570..6e458099 100644 --- a/R/boot.R +++ b/R/boot.R @@ -142,7 +142,7 @@ fect_boot <- function( carryover.period = NULL, vartype = "bootstrap", para.error = "auto", - conformal.scale = "none", + conformal.scale = "sd", conformal.center = "mean", conformal.weight = "cell", conformal.band = "pointwise", diff --git a/R/default.R b/R/default.R index bf26116e..2b3e0f01 100644 --- a/R/default.R +++ b/R/default.R @@ -58,7 +58,7 @@ fect <- function( se = FALSE, # report uncertainties vartype = "bootstrap", # bootstrap or jackknife para.error = "auto", # parametric bootstrap error strategy: "auto", "ar", "empirical", "wild" - conformal.scale = "none", # vartype="conformal": per-unit scale none|sd|rmspe|mad|diff|model-se + conformal.scale = "sd", # vartype="conformal": per-unit scale none|sd|rmspe|mad|diff|model-se conformal.center = "mean", # vartype="conformal": post-period location mean|median conformal.weight = "cell", # vartype="conformal": multi-treated aggregation cell|unit|precision conformal.band = "pointwise", # vartype="conformal": est.att band pointwise|simultaneous @@ -146,7 +146,7 @@ fect.formula <- function( se = FALSE, # report uncertainties vartype = "bootstrap", # bootstrap or jackknife para.error = "auto", # parametric bootstrap error strategy: "auto", "ar", "empirical", "wild" - conformal.scale = "none", # vartype="conformal": per-unit scale none|sd|rmspe|mad|diff|model-se + conformal.scale = "sd", # vartype="conformal": per-unit scale none|sd|rmspe|mad|diff|model-se conformal.center = "mean", # vartype="conformal": post-period location mean|median conformal.weight = "cell", # vartype="conformal": multi-treated aggregation cell|unit|precision conformal.band = "pointwise", # vartype="conformal": est.att band pointwise|simultaneous @@ -356,7 +356,7 @@ fect.default <- function( se = FALSE, # report uncertainties vartype = "bootstrap", # bootstrap or jackknife para.error = "auto", # parametric bootstrap error strategy: "auto", "ar", "empirical", "wild" - conformal.scale = "none", # vartype="conformal": per-unit scale none|sd|rmspe|mad|diff|model-se + conformal.scale = "sd", # vartype="conformal": per-unit scale none|sd|rmspe|mad|diff|model-se conformal.center = "mean", # vartype="conformal": post-period location mean|median conformal.weight = "cell", # vartype="conformal": multi-treated aggregation cell|unit|precision conformal.band = "pointwise", # vartype="conformal": est.att band pointwise|simultaneous diff --git a/man/fect.Rd b/man/fect.Rd index 40b049a2..4fbba516 100644 --- a/man/fect.Rd +++ b/man/fect.Rd @@ -15,7 +15,10 @@ assumptions.} cv.nobs = 3, cv.donut = 1, cv.buffer = 1, criterion = "mspe", binary = FALSE, QR = FALSE, method = "fe", se = FALSE, vartype = "bootstrap", - para.error = "auto", cl = NULL, + para.error = "auto", + conformal.scale = "sd", conformal.center = "mean", + conformal.weight = "cell", conformal.band = "pointwise", + conformal.cutoff = "per-period", cl = NULL, ci.method = "normal", quantile.CI = NULL, nboots = 200, alpha = 0.05, parallel = TRUE, cores = NULL, tol = 1e-5, @@ -83,17 +86,23 @@ In v2.3.1, \code{W.est} and \code{W.agg} (when both supplied) must point to the \code{"gsynth"}, or \code{"cfe"}. Default is \code{"fe"}.} \item{se}{a logical flag indicating whether uncertainty estimates will be produced.} \item{vartype}{a string specifying the type of variance estimator, e.g. - \code{"bootstrap"}. Three values are supported: + \code{"bootstrap"}. Four values are supported: \code{"bootstrap"} (nonparametric cluster-bootstrap; the safe default), - \code{"jackknife"} (leave-one-unit-out), and - \code{"parametric"} (two-stage pseudo-treated parametric bootstrap). + \code{"jackknife"} (leave-one-unit-out), + \code{"parametric"} (two-stage pseudo-treated parametric bootstrap), and + \code{"conformal"} (cross-sectional conformal prediction interval; rank-based + and distribution-free, with no bootstrap draws). The \code{"parametric"} option is restricted to the gsynth-style regime: it requires \code{time.component.from = "nevertreated"}, no treatment reversal, and \code{method} not in \code{c("mc", "both")}. These three conditions correspond to Gates A, B, and C in the three-gate defense system (see \code{ARCHITECTURE.md}). - For all other settings, \code{vartype = "bootstrap"} is recommended.} + \code{"conformal"} likewise requires a separated (controls-only) + fit and \code{method} not in \code{c("mc", "both")}; it currently + supports the core synthetic-control case (no group / reversal / weighted / + balanced / placebo / carryover options). See the \code{conformal.*} + arguments. For all other settings, \code{vartype = "bootstrap"} is recommended.} \item{para.error}{a string specifying the residual-error model used by the parametric bootstrap path; sub-option of \code{vartype = "parametric"} (silently ignored otherwise). One of \code{"auto"}, \code{"ar"}, @@ -105,6 +114,31 @@ In v2.3.1, \code{W.est} and \code{W.agg} (when both supplied) must point to the resamples residuals i.i.d. from the main-fit pool (requires fully- observed panel); \code{"wild"} applies unit-level Rademacher sign-flips over the empirical residual pool (requires fully-observed panel).} +\item{conformal.scale}{for \code{vartype = "conformal"}: the per-unit scale that + normalizes each unit's prediction error before ranking. One of \code{"none"} + (the unscaled mean-absolute statistic), \code{"sd"} (studentize by the unit's + pre-period residual standard deviation), \code{"rmspe"} (Abadie pre-period + root-mean-square error), \code{"mad"} (robust median-absolute-deviation), or + \code{"diff"} (standard deviation of first-differenced pre-period residuals, + robust to nonstationary trends). \code{"model-se"} is reserved and not yet + implemented.} +\item{conformal.center}{for \code{vartype = "conformal"}: the post-period + location summarized into the average effect, \code{"mean"} or \code{"median"} + (the latter is robust to a single aberrant post period). The per-period band + is always reported alongside.} +\item{conformal.weight}{for \code{vartype = "conformal"} with multiple treated + units: how unit effects are aggregated into the average. \code{"cell"} weights + each treated unit-period equally (the conventional ATT; equals + \code{fect}'s \code{att.avg}), \code{"unit"} weights each treated unit equally, + and \code{"precision"} weights by inverse pre-period variance (down-weighting + noisy units).} +\item{conformal.band}{for \code{vartype = "conformal"}: the band reported in + \code{est.att}. \code{"pointwise"} (each period covers on its own) or + \code{"simultaneous"} (a uniform sup-t band over the post window for whole-path + statements). The uniform band is also always stored in \code{fit$est.att.sim}.} +\item{conformal.cutoff}{for \code{vartype = "conformal"} pointwise bands: whether + the rank cutoff is computed \code{"per-period"} (separately at each period) or + \code{"pooled"} (one cutoff over all post-period control gaps).} \item{cl}{a string specifying the cluster column for cluster bootstrapping. When \code{group.fe} is set to a single column and \code{cl} is unset, \code{cl} auto-defaults to \code{group.fe[1]} --- the natural choice when diff --git a/tests/coverage-study/conformal_default_study.R b/tests/coverage-study/conformal_default_study.R new file mode 100644 index 00000000..49b120b8 --- /dev/null +++ b/tests/coverage-study/conformal_default_study.R @@ -0,0 +1,159 @@ +## ----------------------------------------------------------------------------- +## conformal_default_study.R +## +## Picks the default for vartype = "conformal" on evidence. Sweeps the `scale` +## knob (and, separately, the `weight` knob on a staggered DGP) across five data- +## generating processes, recording, for a true effect of zero: +## - scalar att.avg coverage (target 1 - alpha) and median CI width, +## - status != "ok" count (must be 0; Family A never empties), +## - pointwise and simultaneous JOINT (whole-post-path) coverage. +## +## Decision rule: the option that holds scalar coverage across ALL five DGPs at +## the smallest median width is the default. Run: +## Rscript tests/coverage-study/conformal_default_study.R [reps] [alpha] +## Writes a markdown table to tests/coverage-study/results/. +## ----------------------------------------------------------------------------- + +suppressMessages(devtools::load_all("/Users/xyq/GitHub/fect", quiet = TRUE)) + +args <- commandArgs(trailingOnly = TRUE) +reps <- if (length(args) >= 1) as.integer(args[1]) else 150L +alpha <- if (length(args) >= 2) as.numeric(args[2]) else 0.10 +N <- 25L; T <- 20L; T0 <- 15L; r <- 2L +post.terms <- 1:(T - T0) + +## --- DGPs: return a TT x N outcome matrix with true effect 0 (controls only) -- +## Shared latent structure regenerated per rep via the seed stream. +gen <- function(dgp, Fm, L) { + mu <- Fm %*% t(L) + E <- matrix(rnorm(N * T), T, N) + if (dgp == "iid") { + NULL + } else if (dgp == "ar1") { + rho <- 0.5 + for (t in 2:T) E[t, ] <- rho * E[t - 1, ] + sqrt(1 - rho^2) * E[t, ] + } else if (dgp == "hetero") { + s <- rep(c(1, 3), length.out = N) # half the units 3x noisier + E <- E * matrix(s, T, N, byrow = TRUE) + } else if (dgp == "nonstat") { + ramp <- 1 + 1.2 * (seq_len(T) - 1) / (T - 1) # variance grows over time + g <- runif(N, 0.5, 1.5) # unit-specific drift rate + E <- E * outer(ramp, g) + } + mu + E +} + +dgps <- c("iid", "ar1", "hetero", "nonstat") +scales <- c("none", "sd", "rmspe", "mad", "diff") + +joint <- function(m) { + p <- m[as.numeric(rownames(m)) >= 0, c("CI.lower", "CI.upper"), drop = FALSE] + all(p[, 1] <= 0 & 0 <= p[, 2]) +} + +fit_one <- function(Y, D, scale) { + dat <- data.frame(id = rep(1:N, each = T), time = rep(1:T, N), + Y = as.vector(Y), D = as.vector(D)) + suppressWarnings(suppressMessages(fect(Y ~ D, data = dat, index = c("id", "time"), + method = "gsynth", force = 3, CV = FALSE, r = r, se = TRUE, + vartype = "conformal", conformal.scale = scale, alpha = alpha, + parallel = FALSE))) +} + +cat(sprintf("conformal default study: reps=%d alpha=%.2f (target cov %.2f)\n\n", + reps, alpha, 1 - alpha)) + +## ===== Part 1: scale sweep, single treated unit (block) ====================== +D1 <- matrix(0, T, N); D1[(T0 + 1):T, 1] <- 1 +rows <- list() +t0 <- Sys.time() +for (dgp in dgps) { + for (sc in scales) { + set.seed(20260609) # same stream across scales -> paired + cov <- wid <- jp <- js <- numeric(reps); bad <- 0L + for (b in seq_len(reps)) { + Fm <- matrix(rnorm(T * r), T, r); L <- matrix(rnorm(N * r), N, r) + Y <- gen(dgp, Fm, L) + f <- fit_one(Y, D1, sc) + e <- f$est.avg[1, ] + if (f$conformal$status != "ok") bad <- bad + 1L + cov[b] <- (e["CI.lower"] <= 0) && (0 <= e["CI.upper"]) + wid[b] <- e["CI.upper"] - e["CI.lower"] + jp[b] <- joint(f$est.att); js[b] <- joint(f$est.att.sim) + } + rows[[length(rows) + 1]] <- data.frame( + dgp = dgp, scale = sc, cov = mean(cov), med.width = median(wid), + bad = bad, joint.point = mean(jp), joint.sim = mean(js)) + cat(sprintf(" [%-8s %-6s] cov=%.3f width=%.2f bad=%d joint(pt/sim)=%.2f/%.2f\n", + dgp, sc, mean(cov), median(wid), bad, mean(jp), mean(js))) + } +} +tab <- do.call(rbind, rows) +cat(sprintf("\nPart 1 done in %.0fs\n", as.numeric(Sys.time() - t0, units = "secs"))) + +## worst-case scalar coverage per scale across DGPs, and mean width +summ <- do.call(rbind, lapply(scales, function(sc) { + s <- tab[tab$scale == sc, ] + data.frame(scale = sc, worst.cov = min(s$cov), mean.cov = mean(s$cov), + mean.width = mean(s$med.width), max.bad = max(s$bad)) +})) +summ <- summ[order(-summ$worst.cov, summ$mean.width), ] +cat("\n=== scale ranking (worst-case scalar coverage, then width) ===\n") +print(summ, row.names = FALSE, digits = 3) +best <- summ$scale[1] +cat(sprintf("\nProvisional default scale (most robust at smallest width): %s\n", best)) + +## ===== Part 2: weight comparison on a staggered DGP ========================== +set.seed(424242) +cat("\n=== Part 2: weight sweep, staggered + heteroskedastic treated ===\n") +Ds <- matrix(0, T, N); Ds[(T0 - 2 + 1):T, 1] <- 1; Ds[(T0 + 1):T, 2] <- 1; Ds[(T0 + 3):T, 3] <- 1 +wrows <- list() +for (w in c("cell", "unit", "precision")) { + cov <- wid <- numeric(reps); bad <- 0L + set.seed(13) + for (b in seq_len(reps)) { + Fm <- matrix(rnorm(T * r), T, r); L <- matrix(rnorm(N * r), N, r) + E <- matrix(rnorm(N * T), T, N); E[, 1] <- E[, 1] * 3 # treated 1 noisy + Y <- Fm %*% t(L) + E + dat <- data.frame(id = rep(1:N, each = T), time = rep(1:T, N), + Y = as.vector(Y), D = as.vector(Ds)) + f <- tryCatch(suppressWarnings(suppressMessages(fect(Y ~ D, data = dat, + index = c("id", "time"), method = "gsynth", force = 3, CV = FALSE, + r = r, se = TRUE, vartype = "conformal", conformal.weight = w, + alpha = alpha, parallel = FALSE))), error = function(e) NULL) + if (is.null(f)) { bad <- bad + 1L; cov[b] <- NA; wid[b] <- NA; next } + e <- f$est.avg[1, ] + cov[b] <- (e["CI.lower"] <= 0) && (0 <= e["CI.upper"]); wid[b] <- e["CI.upper"] - e["CI.lower"] + } + wrows[[length(wrows) + 1]] <- data.frame(weight = w, cov = mean(cov, na.rm = TRUE), + med.width = median(wid, na.rm = TRUE), bad = bad) + cat(sprintf(" [weight=%-9s] cov=%.3f width=%.2f bad=%d\n", w, + mean(cov, na.rm = TRUE), median(wid, na.rm = TRUE), bad)) +} +wtab <- do.call(rbind, wrows) + +## ===== write results ========================================================= +out <- "/Users/xyq/GitHub/fect/tests/coverage-study/results/conformal_default_study.md" +con <- file(out, "w") +writeLines(c( + "# Conformal default study", + sprintf("reps=%d, alpha=%.2f, target coverage=%.2f, N=%d T=%d T0=%d r=%d, true effect 0.", + reps, alpha, 1 - alpha, N, T, T0, r), + "", "## Part 1: scale sweep (single treated, block)", "", + "| dgp | scale | cov | med.width | bad | joint.point | joint.sim |", + "|---|---|---|---|---|---|---|"), con) +for (i in seq_len(nrow(tab))) writeLines(sprintf("| %s | %s | %.3f | %.2f | %d | %.2f | %.2f |", + tab$dgp[i], tab$scale[i], tab$cov[i], tab$med.width[i], tab$bad[i], + tab$joint.point[i], tab$joint.sim[i]), con) +writeLines(c("", "## Scale ranking (worst-case scalar coverage, then width)", "", + "| scale | worst.cov | mean.cov | mean.width | max.bad |", + "|---|---|---|---|---|"), con) +for (i in seq_len(nrow(summ))) writeLines(sprintf("| %s | %.3f | %.3f | %.2f | %d |", + summ$scale[i], summ$worst.cov[i], summ$mean.cov[i], summ$mean.width[i], summ$max.bad[i]), con) +writeLines(c("", sprintf("**Provisional default scale: %s** (most robust at smallest width).", best), + "", "## Part 2: weight sweep (staggered + heteroskedastic treated)", "", + "| weight | cov | med.width | bad |", "|---|---|---|---|"), con) +for (i in seq_len(nrow(wtab))) writeLines(sprintf("| %s | %.3f | %.2f | %d |", + wtab$weight[i], wtab$cov[i], wtab$med.width[i], wtab$bad[i]), con) +close(con) +cat(sprintf("\nResults written to %s\n", out)) diff --git a/tests/coverage-study/results/conformal_default_study.md b/tests/coverage-study/results/conformal_default_study.md new file mode 100644 index 00000000..baca004e --- /dev/null +++ b/tests/coverage-study/results/conformal_default_study.md @@ -0,0 +1,47 @@ +# Conformal default study +reps=150, alpha=0.10, target coverage=0.90, N=25 T=20 T0=15 r=2, true effect 0. + +## Part 1: scale sweep (single treated, block) + +| dgp | scale | cov | med.width | bad | joint.point | joint.sim | +|---|---|---|---|---|---|---| +| iid | none | 0.940 | 2.15 | 0 | 0.59 | 0.84 | +| iid | sd | 0.927 | 2.41 | 0 | 0.71 | 0.83 | +| iid | rmspe | 0.927 | 2.41 | 0 | 0.71 | 0.83 | +| iid | mad | 0.913 | 2.63 | 0 | 0.69 | 0.83 | +| iid | diff | 0.927 | 2.51 | 0 | 0.67 | 0.83 | +| ar1 | none | 0.913 | 3.07 | 0 | 0.61 | 0.83 | +| ar1 | sd | 0.927 | 3.34 | 0 | 0.67 | 0.87 | +| ar1 | rmspe | 0.927 | 3.34 | 0 | 0.67 | 0.87 | +| ar1 | mad | 0.920 | 3.66 | 0 | 0.73 | 0.89 | +| ar1 | diff | 0.893 | 3.27 | 0 | 0.69 | 0.87 | +| hetero | none | 1.000 | 5.43 | 0 | 0.97 | 0.99 | +| hetero | sd | 0.920 | 2.82 | 0 | 0.68 | 0.86 | +| hetero | rmspe | 0.920 | 2.82 | 0 | 0.68 | 0.86 | +| hetero | mad | 0.920 | 3.16 | 0 | 0.69 | 0.87 | +| hetero | diff | 0.907 | 3.03 | 0 | 0.66 | 0.85 | +| nonstat | none | 0.920 | 4.51 | 0 | 0.63 | 0.85 | +| nonstat | sd | 0.933 | 4.05 | 0 | 0.67 | 0.85 | +| nonstat | rmspe | 0.933 | 4.05 | 0 | 0.67 | 0.85 | +| nonstat | mad | 0.900 | 4.57 | 0 | 0.66 | 0.80 | +| nonstat | diff | 0.940 | 4.41 | 0 | 0.65 | 0.85 | + +## Scale ranking (worst-case scalar coverage, then width) + +| scale | worst.cov | mean.cov | mean.width | max.bad | +|---|---|---|---|---| +| sd | 0.920 | 0.927 | 3.16 | 0 | +| rmspe | 0.920 | 0.927 | 3.16 | 0 | +| none | 0.913 | 0.943 | 3.79 | 0 | +| mad | 0.900 | 0.913 | 3.50 | 0 | +| diff | 0.893 | 0.917 | 3.30 | 0 | + +**Provisional default scale: sd** (most robust at smallest width). + +## Part 2: weight sweep (staggered + heteroskedastic treated) + +| weight | cov | med.width | bad | +|---|---|---|---| +| cell | 0.873 | 1.99 | 0 | +| unit | 0.873 | 1.99 | 0 | +| precision | 0.993 | 1.99 | 0 | From a9e4a8354df1f961ffdf6aa89bb7a1578b2e9f81 Mon Sep 17 00:00:00 2001 From: Yiqing Xu Date: Tue, 9 Jun 2026 08:39:02 -0700 Subject: [PATCH 09/14] feat(conformal): staggered adoption support (per-cohort, event-time bands) The calibration now handles staggered onsets correctly, not just block / common-onset. Block behavior is preserved exactly (verified by the existing tests). ## conformal.R - Treatment-incidence structure (nt[t], union post window, per-unit onsets); control leave-one-out fits mask the UNION post window so donor gaps are held-out wherever any treated unit is post. - Per-treated-unit summaries over each unit's OWN post/pre window (was the union window for every unit, which mixed post-effect and pre-residual cells under staggering). Treated trajectory is the per-period weighted mean over the units that are post (or, where none are, pre) at t. - Control placebo center = the control's gap aggregated over the union post window the same nt-weighted way as the treated ATT (block: a plain mean). - New EVENT-TIME band: aggregates cohorts at each relative period, pooling the control gaps at the matching calendar cells. For block it equals the calendar band re-indexed. Returned as band.et / band.et.sim alongside the calendar band. ## boot.R - est.att / est.att.sim now come from the event-time band (correct for staggered), matched to out$time by relative-time value, replacing the single-unit T.on map. ## Verification - Staggered coverage (2 cohorts, 100 reps, true effect 0, alpha=0.10): scalar 0.900 (iid) / 0.940 (ar1); per-period band 0.911 / 0.906; bad=0. - 58 conformal tests pass (added a staggered run/alignment/coverage test); all block tests unchanged. Refs statsclaw runs/2026-06-08-conformal-inference.md. Co-Authored-By: Claude Opus 4.8 --- R/boot.R | 26 +++-- R/conformal.R | 171 ++++++++++++++++++-------------- tests/testthat/test-conformal.R | 25 +++++ 3 files changed, 134 insertions(+), 88 deletions(-) diff --git a/R/boot.R b/R/boot.R index 6e458099..95299abd 100644 --- a/R/boot.R +++ b/R/boot.R @@ -568,7 +568,6 @@ fect_boot <- function( "'bootstrap' or 'jackknife' for those.", call. = FALSE) } - id.tr.conf <- which(colSums(D) > 0) cc <- conformal_calibrate( Y = Y, D = D, X = X, I = I, II = II, T.on = T.on, r.cv = out$r.cv, eff = out$eff, @@ -606,26 +605,25 @@ fect_boot <- function( band[, "CI.lower"], band[, "CI.upper"], band[, "p.value"], calendar.N) colnames(est.eff.calendar.fit) <- colnames(est.eff.calendar) - ## --- map a calendar band to event time (treated unit's relative period); - ## rownames = out$time. Common-onset: a 1:1 map. - rel <- T.on[, id.tr.conf[1]] - map_band <- function(bnd) { + ## --- event-time band -> est.att (rownames = out$time). conformal_calibrate + ## returns the band already aggregated by relative period (correct for both + ## block and staggered); match its rownames to out$time by value. + map_et <- function(etb) { m <- matrix(NA_real_, length(time.on), 6, dimnames = list(time.on, c("ATT", "S.E.", "CI.lower", "CI.upper", "p.value", "count"))) + etrow <- as.numeric(rownames(etb)) for (k in seq_along(time.on)) { - rows <- which(rel == time.on[k]) - if (length(rows) == 0L) next - ef <- mean(bnd[rows, "eff"], na.rm = TRUE) - lo <- mean(bnd[rows, "CI.lower"], na.rm = TRUE) - hi <- mean(bnd[rows, "CI.upper"], na.rm = TRUE) - pv <- mean(bnd[rows, "p.value"], na.rm = TRUE) - m[k, ] <- c(ef, se.from(lo, hi), lo, hi, pv, out$count[k]) + r <- which(etrow == time.on[k]) + if (!length(r)) next + v <- etb[r[1L], ] + m[k, ] <- c(v["eff"], se.from(v["CI.lower"], v["CI.upper"]), + v["CI.lower"], v["CI.upper"], v["p.value"], out$count[k]) } m } - est.att <- map_band(band) # the selected band (conformal.band) - est.att.sim <- map_band(cc$band.sim) # the uniform / simultaneous band + est.att <- map_et(cc$band.et) # the selected band (conformal.band) + est.att.sim <- map_et(cc$band.et.sim) # the uniform / simultaneous band ## Phase 4 note: est.att90 (the inner equivalence band) mirrors the main band ## for now; a true (1 - 2*alpha) conformal inner band is a later refinement. est.att90 <- est.att diff --git a/R/conformal.R b/R/conformal.R index 975fc088..c86ab65a 100644 --- a/R/conformal.R +++ b/R/conformal.R @@ -149,19 +149,21 @@ conformal_calibrate <- function(Y, D, X = NULL, I, II, T.on, r.cv, eff, stop("conformal: fewer than 2 valid controls after screening.") } - ## --- post / pre period windows (block: common onset; staggered: union) - post.idx <- which(rowSums(D[, id.tr, drop = FALSE]) > 0) # any treated on - pre.idx <- setdiff(seq_len(TT), post.idx) + ## --- treatment-incidence structure (handles block AND staggered uniformly). + ## post.idx is the UNION post window (any treated unit on); nt[t] is the number + ## of treated units in post at calendar t; onsets are per-unit first-treated + ## periods. For a control's leave-one-out fit, the union window is masked so its + ## gaps are held-out everywhere a treated unit is post (block: a single onset). + Dtr <- D[, id.tr, drop = FALSE] # TT x Ntr, 1 = post + nt <- rowSums(Dtr) # treated post count / period + union.post <- nt > 0 + post.idx <- which(union.post) + pre.idx <- setdiff(seq_len(TT), post.idx) if (length(post.idx) == 0L || length(pre.idx) == 0L) { stop("conformal: could not identify pre/post windows from D.") } - if (Ntr > 1L) { - onsets <- vapply(id.tr, function(i) min(which(D[, i] == 1)), integer(1)) - if (length(unique(onsets)) > 1L) { - warning("conformal: staggered onsets detected; using the union post-window. ", - "Per-cohort windows are not yet implemented.") - } - } + onsets <- apply(Dtr, 2L, function(d) { w <- which(d == 1); if (length(w)) min(w) else NA_integer_ }) + staggered <- length(unique(onsets[is.finite(onsets)])) > 1L ## --- one held-out fit per valid control, reusing the parametric refit path. ## fake-treated column = held-out control j's DATA carrying a treated unit's @@ -178,7 +180,9 @@ conformal_calibrate <- function(Y, D, X = NULL, I, II, T.on, r.cv, eff, D.ps <- cbind(d.pattern, D[, co.rest, drop = FALSE]) Ton.ps <- cbind(T.on[, id.tr[1]], T.on[, co.rest, drop = FALSE]) I.ps <- cbind(I[, j], I[, co.rest, drop = FALSE]) - II.ps <- cbind(I[, j] * (d.pattern == 0), II[, co.rest, drop = FALSE]) + ## mask the UNION post window so the control's gaps are held-out wherever any + ## treated unit is post (block: a single onset; staggered: the union). + II.ps <- cbind(I[, j] * as.numeric(!union.post), II[, co.rest, drop = FALSE]) X.ps <- if (is.null(X)) NULL else array(c(X[, j, , drop = FALSE], X[, co.rest, , drop = FALSE]), dim = c(TT, 1L + length(co.rest), dim(X)[3])) @@ -209,85 +213,75 @@ conformal_calibrate <- function(Y, D, X = NULL, I, II, T.on, r.cv, eff, sc.co[mi] <- .conformal_scale(gj[pre.idx], scale) } - ## --- treated aggregate across treated units, combined per `weight`: - ## cell = weight each unit by its post-obs count (per-treated-cell ATT), - ## unit = each treated unit counts equally, - ## precision = inverse pre-period variance (down-weights noisy units). - ## The scalar center is the weighted mean of the per-unit post centers; the - ## per-period band uses the matching weighted treated trajectory. They coincide - ## for a single treated unit and (mean center) for a balanced panel. - ## Multi-treated NOTE: the aggregate is ranked against single-control gaps; a - ## placebo-AVERAGE calibration (ranking against averages of Ntr controls) is a - ## refinement deferred to the simulation phase. - trc <- eff[, id.tr, drop = FALSE] # TT x Ntr - m.i <- apply(trc[post.idx, , drop = FALSE], 2L, function(z) .conformal_center(z, center)) - n.i <- colSums(!is.na(trc[post.idx, , drop = FALSE])) - ## precision weight uses each unit's pre-period VARIANCE directly (decoupled - ## from the `scale` knob, which only normalizes the score), so it down-weights - ## noisy units even under the default scale = "none". - var.i <- apply(trc[pre.idx, , drop = FALSE], 2L, function(z) { - z <- z[is.finite(z)]; if (length(z) > 1L) stats::var(z) else NA_real_ - }) - wi <- switch(weight, - cell = n.i, - unit = rep(1, Ntr), + ## --- per-treated-unit summaries, each over the unit's OWN post / pre window + ## (block: every unit shares the union window; staggered: per-unit windows). + ## `weight`: cell = per-treated-cell ATT (n_i); unit = equal per unit; + ## precision = inverse pre-period variance (decoupled from the score `scale`). + trc <- eff[, id.tr, drop = FALSE] # TT x Ntr + obstr <- (I[, id.tr, drop = FALSE] == 1) + postmask.i <- (Dtr == 1) # own post cells + premask.i <- (Dtr == 0) & obstr # own observed pre cells + m.i <- vapply(seq_len(Ntr), function(i) + .conformal_center(trc[postmask.i[, i], i], center), numeric(1)) + n.i <- colSums(postmask.i) + var.i <- vapply(seq_len(Ntr), function(i) { + z <- trc[premask.i[, i], i]; z <- z[is.finite(z)] + if (length(z) > 1L) stats::var(z) else NA_real_ }, numeric(1)) + wi <- switch(weight, cell = n.i, unit = rep(1, Ntr), precision = ifelse(is.finite(var.i) & var.i > 0, 1 / var.i, 0)) wi[!is.finite(wi)] <- 0 if (!any(wi > 0)) wi <- rep(1, Ntr) - Wt <- matrix(wi, TT, Ntr, byrow = TRUE); Wt[is.na(trc)] <- 0 - denom <- rowSums(Wt) - tr.path <- ifelse(denom > 0, rowSums(trc * Wt, na.rm = TRUE) / denom, NA_real_) # calendar + + ## treated trajectory (calendar): weighted mean over the units that are POST at + ## t (effect window) or, where none are, over the units that are PRE at t + ## (placebo window). For block this is the per-period weighted mean throughout. + wsum <- function(e, w) { ok <- is.finite(e) & is.finite(w) & w > 0 + if (any(ok)) sum(e[ok] * w[ok]) / sum(w[ok]) else NA_real_ } + Wt <- matrix(wi, TT, Ntr, byrow = TRUE) + postW <- Wt * postmask.i; preW <- Wt * premask.i + tr.path <- vapply(seq_len(TT), function(t) + wsum(trc[t, ], if (nt[t] > 0) postW[t, ] else preW[t, ]), numeric(1)) sc.tr <- .conformal_scale(tr.path[pre.idx], scale) w.co <- NULL # overlap density-ratio weighting deferred to a later phase - ## --- scalar average effect: weighted mean of per-unit post centers, ranked vs - ## the per-control post centers. - m.co <- apply(G.co[, post.idx, drop = FALSE], 1L, - function(z) .conformal_center(z, center)) - keep <- is.finite(m.co) & is.finite(sc.co) & sc.co > 0 - okm <- is.finite(m.i) & wi > 0 - m.tr <- if (any(okm)) sum(wi[okm] * m.i[okm]) / sum(wi[okm]) else - .conformal_center(tr.path[post.idx], center) - ci <- .conformal_ci_level(m.tr, sc.tr, m.co[keep], sc.co[keep], - alpha = alpha, w.co = w.co) - Ncal <- sum(keep) + ## --- scalar average effect: weighted mean of per-unit post centers, ranked + ## against per-control placebo centers (each control's gap aggregated over the + ## union post window the SAME way the treated ATT is, i.e. nt-weighted = the + ## per-treated-cell ATT; for block nt is constant so this is a plain mean). + ntp <- nt[post.idx] + m.co <- apply(G.co, 1L, function(g) wsum(g[post.idx], ntp)) + keep <- is.finite(m.co) & is.finite(sc.co) & sc.co > 0 + okm <- is.finite(m.i) & wi > 0 + m.tr <- if (any(okm)) sum(wi[okm] * m.i[okm]) / sum(wi[okm]) else + .conformal_center(tr.path[post.idx], center) + ci <- .conformal_ci_level(m.tr, sc.tr, m.co[keep], sc.co[keep], alpha = alpha, w.co = w.co) + Ncal <- sum(keep) status <- if (any(is.infinite(ci))) "unbounded" else "ok" s.tr <- if (is.finite(sc.tr) && sc.tr > 0) abs(m.tr) / sc.tr else NA_real_ s.co.vec <- abs(m.co[keep]) / sc.co[keep] p.value <- .conformal_pval(s.tr, s.co.vec, w.co = w.co) - ## --- per-period band(s) (calendar-indexed). Standardized control gaps - ## g.tilde[j,t] = |G.co[j,t]| / sc.co[j]; the treated deviation at period t is - ## |tr.path[t] - tau| / sc.tr. Three constructions: - ## pointwise + per-period (default): Q_t = conformal quantile of g.tilde[,t]. - ## pointwise + pooled: one Q over all post-window g.tilde values. - ## simultaneous: one multiplier c = quantile of the per-control MAX over the - ## post window (sup-t / uniform band over the post path); pre-period rows - ## keep the per-period placebo band. - ## The pre-period rows always show the per-period placebo band (for pre-trends). + ## --- standardized control gaps g.tilde[j,t] = |G.co[j,t]| / sc.co[j], and the + ## simultaneous (sup-t) multiplier from the per-control MAX over the post window. ok.co <- which(is.finite(sc.co) & sc.co > 0) - Gv <- G.co[ok.co, , drop = FALSE] - scv <- sc.co[ok.co] - Gtil <- abs(Gv) / scv # ncal x TT, row j divided by scv[j] - - Qpool <- NULL - if (identical(cutoff, "pooled")) { - pv <- as.vector(Gtil[, post.idx, drop = FALSE]); pv <- pv[is.finite(pv)] - Qpool <- .conformal_quantile(pv, alpha) - } + Gv <- G.co[ok.co, , drop = FALSE]; scv <- sc.co[ok.co] + Gtil <- abs(Gv) / scv # ncal x TT Mj <- apply(Gtil[, post.idx, drop = FALSE], 1L, function(z) { z <- z[is.finite(z)]; if (length(z)) max(z) else NA_real_ }) c.sim <- .conformal_quantile(Mj[is.finite(Mj)], alpha) + Qpool <- if (identical(cutoff, "pooled")) { + pv <- as.vector(Gtil[, post.idx, drop = FALSE]); .conformal_quantile(pv[is.finite(pv)], alpha) + } else NULL - mk_band <- function(mode) { - b <- matrix(NA_real_, TT, 4L, - dimnames = list(NULL, c("eff", "CI.lower", "CI.upper", "p.value"))) + ## --- CALENDAR band (per calendar period) -> est.eff.calendar. + mk_cal <- function(mode) { + b <- matrix(NA_real_, TT, 4L, dimnames = list(NULL, c("eff", "CI.lower", "CI.upper", "p.value"))) for (t in seq_len(TT)) { gt <- Gtil[, t]; okt <- is.finite(gt) if (!is.finite(tr.path[t]) || !is.finite(sc.tr) || sc.tr <= 0 || sum(okt) < 2L) next - Qt <- if (mode == "simultaneous" && t %in% post.idx) c.sim - else if (identical(cutoff, "pooled")) Qpool + Qt <- if (mode == "simultaneous" && union.post[t]) c.sim + else if (!is.null(Qpool) && union.post[t]) Qpool else .conformal_quantile(gt[okt], alpha) half <- if (is.finite(Qt)) sc.tr * Qt else Inf b[t, ] <- c(tr.path[t], tr.path[t] - half, tr.path[t] + half, @@ -295,13 +289,42 @@ conformal_calibrate <- function(Y, D, X = NULL, I, II, T.on, r.cv, eff, } b } - band <- mk_band(band.type) # populates est.att - band.sim <- if (identical(band.type, "simultaneous")) band else mk_band("simultaneous") + band <- mk_cal(band.type) + band.sim <- if (identical(band.type, "simultaneous")) band else mk_cal("simultaneous") + + ## --- EVENT-TIME band (per relative period) -> est.att. For block this is the + ## calendar band re-indexed; for staggered it aggregates cohorts at each relative + ## time, pooling the control gaps at the matching calendar cells (cross-sectional + ## exchangeability across units, with stationarity across cohort onsets). + rel.mat <- T.on[, id.tr, drop = FALSE] + etimes <- sort(unique(rel.mat[obstr])) + mk_et <- function(mode) { + b <- matrix(NA_real_, length(etimes), 4L, + dimnames = list(as.character(etimes), c("eff", "CI.lower", "CI.upper", "p.value"))) + for (k in seq_along(etimes)) { + cells <- which(rel.mat == etimes[k] & obstr, arr.ind = TRUE) # (t, unit) pairs + if (!nrow(cells)) next + tg <- wsum(trc[cells], wi[cells[, 2L]]) + ts <- unique(cells[, 1L]) # calendar periods at this rel time + ispost <- any(union.post[ts]) + gp <- as.vector(Gtil[, ts, drop = FALSE]); gp <- gp[is.finite(gp)] + if (!is.finite(tg) || !is.finite(sc.tr) || sc.tr <= 0 || length(gp) < 2L) next + Qk <- if (mode == "simultaneous" && ispost) c.sim + else if (!is.null(Qpool) && ispost) Qpool + else .conformal_quantile(gp, alpha) + half <- if (is.finite(Qk)) sc.tr * Qk else Inf + b[k, ] <- c(tg, tg - half, tg + half, .conformal_pval(abs(tg) / sc.tr, gp)) + } + b + } + band.et <- mk_et(band.type) + band.et.sim <- if (identical(band.type, "simultaneous")) band.et else mk_et("simultaneous") list(att = m.tr, ci = c(ci[1], ci[2]), p.value = p.value, score.tr = s.tr, score.co = s.co.vec, status = status, n.calib = Ncal, n.eff = if (is.null(w.co)) Ncal else .conformal_neff(w.co), scale = scale, center = center, weight = weight, - band.type = band.type, cutoff = cutoff, form = "jackknife+", - band = band, band.sim = band.sim, post.idx = post.idx, pre.idx = pre.idx) + band.type = band.type, cutoff = cutoff, staggered = staggered, form = "jackknife+", + band = band, band.sim = band.sim, band.et = band.et, band.et.sim = band.et.sim, + post.idx = post.idx, pre.idx = pre.idx) } diff --git a/tests/testthat/test-conformal.R b/tests/testthat/test-conformal.R index 89d1d4b4..a6589c36 100644 --- a/tests/testthat/test-conformal.R +++ b/tests/testthat/test-conformal.R @@ -215,6 +215,31 @@ test_that("simultaneous band is wider than pointwise; est.att.sim always present c("CI.lower", "CI.upper")]))) }) +## ---- staggered adoption ----------------------------------------------------- + +test_that("conformal supports staggered adoption: runs, aligns event time, covers", { + skip_on_cran() + N <- 40; T <- 20; r <- 2 + D <- matrix(0, T, N); D[12:T, 1:4] <- 1; D[16:T, 5:8] <- 1 # cohorts at 12 and 16 + set.seed(20260609); reps <- 40; alpha <- 0.10 + cov <- logical(reps); bad <- 0L; aligned <- TRUE + for (b in seq_len(reps)) { + Fm <- matrix(rnorm(T * r), T, r); L <- matrix(rnorm(N * r), N, r) + Y <- Fm %*% t(L) + matrix(rnorm(N * T), T, N) # true effect 0 + dat <- data.frame(id = rep(1:N, each = T), time = rep(1:T, N), + Y = as.vector(Y), D = as.vector(D)) + f <- suppressWarnings(suppressMessages(fect(Y ~ D, data = dat, index = c("id", "time"), + method = "gsynth", force = 3, CV = FALSE, r = r, se = TRUE, + vartype = "conformal", conformal.scale = "sd", alpha = alpha, parallel = FALSE))) + if (f$conformal$status != "ok") bad <- bad + 1L + cov[b] <- (f$est.avg[1, "CI.lower"] <= 0) && (0 <= f$est.avg[1, "CI.upper"]) + if (!identical(rownames(f$est.att), as.character(f$time))) aligned <- FALSE + } + expect_true(aligned) # est.att indexed by event time, matched to fit$time + expect_equal(bad, 0L) # never empty/unbounded + expect_gt(mean(cov), 0.80) # scalar coverage near nominal under staggering +}) + test_that("simultaneous band gives materially better joint coverage than pointwise", { skip_on_cran() ## Joint (whole-post-path) coverage. Pointwise undercovers jointly (the From bdd8543755e88b7f38eb8a621521992dc9e2f631 Mon Sep 17 00:00:00 2001 From: Yiqing Xu Date: Tue, 9 Jun 2026 08:42:36 -0700 Subject: [PATCH 10/14] test(conformal): plotting smoke for block + staggered fits plot.fect already renders conformal fits correctly (gap + counterfactual), and the simultaneous band is plottable via conformal.band = "simultaneous" (it populates est.att) or the est.att.sim slot -- no plot.R change needed. Adds a test that all plot types build a ggplot for block + staggered + simultaneous fits, and that the uniform band is wider than the pointwise band post-treatment. Verification figures (block/staggered gap + counterfactual, pointwise-vs- simultaneous overlay) saved under the sc-conformal paper folder. Co-Authored-By: Claude Opus 4.8 --- tests/testthat/test-conformal.R | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/tests/testthat/test-conformal.R b/tests/testthat/test-conformal.R index a6589c36..a43bdfdd 100644 --- a/tests/testthat/test-conformal.R +++ b/tests/testthat/test-conformal.R @@ -215,6 +215,36 @@ test_that("simultaneous band is wider than pointwise; est.att.sim always present c("CI.lower", "CI.upper")]))) }) +## ---- plotting --------------------------------------------------------------- + +test_that("plot() works on conformal fits (block + staggered, gap + counterfactual)", { + skip_on_cran() + mkfit <- function(D, seed, ...) { + set.seed(seed); N <- ncol(D); T <- nrow(D); r <- 2 + Fm <- matrix(rnorm(T * r), T, r); L <- matrix(rnorm(N * r), N, r) + Y <- Fm %*% t(L) + matrix(rnorm(N * T), T, N) + 1.5 * D + dat <- data.frame(id = rep(1:N, each = T), time = rep(1:T, N), + Y = as.vector(Y), D = as.vector(D)) + suppressWarnings(suppressMessages(fect(Y ~ D, data = dat, index = c("id", "time"), + method = "gsynth", force = 3, CV = FALSE, r = r, se = TRUE, + vartype = "conformal", parallel = FALSE, ...))) + } + Db <- matrix(0, 20, 25); Db[16:20, 1] <- 1 + Ds <- matrix(0, 20, 40); Ds[12:20, 1:4] <- 1; Ds[16:20, 5:8] <- 1 + fb <- mkfit(Db, 1) + fbs <- mkfit(Db, 1, conformal.band = "simultaneous") + fs <- mkfit(Ds, 2) + for (f in list(fb, fbs, fs)) { + expect_s3_class(plot(f, type = "gap"), "gg") + expect_s3_class(plot(f, type = "counterfactual"), "gg") + } + ## the uniform band is stored and wider than the pointwise band in the post window + post <- as.numeric(rownames(fb$est.att)) >= 0 + w_pt <- mean(fb$est.att[post, "CI.upper"] - fb$est.att[post, "CI.lower"]) + w_sim <- mean(fb$est.att.sim[post, "CI.upper"] - fb$est.att.sim[post, "CI.lower"]) + expect_gte(w_sim, w_pt) +}) + ## ---- staggered adoption ----------------------------------------------------- test_that("conformal supports staggered adoption: runs, aligns event time, covers", { From b49bfa3fdff7e1e73b07c79ce6c7dc16142fff36 Mon Sep 17 00:00:00 2001 From: Yiqing Xu Date: Tue, 9 Jun 2026 08:46:51 -0700 Subject: [PATCH 11/14] feat(conformal): true inner (1-2*alpha) band for est.att90 est.att90 (the inner equivalence band shown in plots) was a placeholder that mirrored the main band. It is now a genuine narrower conformal band at level 1 - 2*alpha, computed from the same control scores at no extra fitting cost (the band builders are parameterized by level; the outer alpha and inner 2*alpha bands reuse one calibration pass). att.bound now derives from this inner band. Verification: alpha=0.05 gives a 95% main band and a narrower 90% inner band (widths 5.20 vs 4.02); 61 conformal tests pass; block bands at alpha unchanged. Co-Authored-By: Claude Opus 4.8 --- R/boot.R | 10 +++--- R/conformal.R | 54 +++++++++++++++++++-------------- tests/testthat/test-conformal.R | 3 ++ 3 files changed, 39 insertions(+), 28 deletions(-) diff --git a/R/boot.R b/R/boot.R index 95299abd..14254df6 100644 --- a/R/boot.R +++ b/R/boot.R @@ -622,12 +622,10 @@ fect_boot <- function( } m } - est.att <- map_et(cc$band.et) # the selected band (conformal.band) - est.att.sim <- map_et(cc$band.et.sim) # the uniform / simultaneous band - ## Phase 4 note: est.att90 (the inner equivalence band) mirrors the main band - ## for now; a true (1 - 2*alpha) conformal inner band is a later refinement. - est.att90 <- est.att - att.bound <- est.att[, c("CI.lower", "CI.upper"), drop = FALSE] + est.att <- map_et(cc$band.et) # the (1 - alpha) band (conformal.band) + est.att.sim <- map_et(cc$band.et.sim) # the uniform / simultaneous band + est.att90 <- map_et(cc$band.et.inner) # the inner (1 - 2*alpha) conformal band + att.bound <- est.att90[, c("CI.lower", "CI.upper"), drop = FALSE] rownames(att.bound) <- time.on result <- list( diff --git a/R/conformal.R b/R/conformal.R index c86ab65a..1db44c84 100644 --- a/R/conformal.R +++ b/R/conformal.R @@ -262,35 +262,41 @@ conformal_calibrate <- function(Y, D, X = NULL, I, II, T.on, r.cv, eff, s.co.vec <- abs(m.co[keep]) / sc.co[keep] p.value <- .conformal_pval(s.tr, s.co.vec, w.co = w.co) - ## --- standardized control gaps g.tilde[j,t] = |G.co[j,t]| / sc.co[j], and the - ## simultaneous (sup-t) multiplier from the per-control MAX over the post window. - ok.co <- which(is.finite(sc.co) & sc.co > 0) - Gv <- G.co[ok.co, , drop = FALSE]; scv <- sc.co[ok.co] - Gtil <- abs(Gv) / scv # ncal x TT - Mj <- apply(Gtil[, post.idx, drop = FALSE], 1L, - function(z) { z <- z[is.finite(z)]; if (length(z)) max(z) else NA_real_ }) - c.sim <- .conformal_quantile(Mj[is.finite(Mj)], alpha) - Qpool <- if (identical(cutoff, "pooled")) { - pv <- as.vector(Gtil[, post.idx, drop = FALSE]); .conformal_quantile(pv[is.finite(pv)], alpha) - } else NULL + ## --- standardized control gaps g.tilde[j,t] = |G.co[j,t]| / sc.co[j], the + ## per-control MAX over the post window (for the sup-t multiplier), and the + ## pooled post-window values (for cutoff = "pooled"). The band builders take a + ## level `a` so the outer (alpha) and inner (2*alpha) bands reuse one pass. + ok.co <- which(is.finite(sc.co) & sc.co > 0) + Gv <- G.co[ok.co, , drop = FALSE]; scv <- sc.co[ok.co] + Gtil <- abs(Gv) / scv # ncal x TT + Mj <- apply(Gtil[, post.idx, drop = FALSE], 1L, + function(z) { z <- z[is.finite(z)]; if (length(z)) max(z) else NA_real_ }) + Mj <- Mj[is.finite(Mj)] + pvpool <- if (identical(cutoff, "pooled")) { + v <- as.vector(Gtil[, post.idx, drop = FALSE]); v[is.finite(v)] + } else NULL + qt_post <- function(gp, a) { # post-window cutoff at level a + if (!is.null(pvpool)) .conformal_quantile(pvpool, a) else .conformal_quantile(gp, a) + } ## --- CALENDAR band (per calendar period) -> est.eff.calendar. - mk_cal <- function(mode) { + mk_cal <- function(mode, a) { + c.sim <- .conformal_quantile(Mj, a) b <- matrix(NA_real_, TT, 4L, dimnames = list(NULL, c("eff", "CI.lower", "CI.upper", "p.value"))) for (t in seq_len(TT)) { gt <- Gtil[, t]; okt <- is.finite(gt) if (!is.finite(tr.path[t]) || !is.finite(sc.tr) || sc.tr <= 0 || sum(okt) < 2L) next Qt <- if (mode == "simultaneous" && union.post[t]) c.sim - else if (!is.null(Qpool) && union.post[t]) Qpool - else .conformal_quantile(gt[okt], alpha) + else if (union.post[t]) qt_post(gt[okt], a) + else .conformal_quantile(gt[okt], a) half <- if (is.finite(Qt)) sc.tr * Qt else Inf b[t, ] <- c(tr.path[t], tr.path[t] - half, tr.path[t] + half, .conformal_pval(abs(tr.path[t]) / sc.tr, gt[okt])) } b } - band <- mk_cal(band.type) - band.sim <- if (identical(band.type, "simultaneous")) band else mk_cal("simultaneous") + band <- mk_cal(band.type, alpha) + band.sim <- if (identical(band.type, "simultaneous")) band else mk_cal("simultaneous", alpha) ## --- EVENT-TIME band (per relative period) -> est.att. For block this is the ## calendar band re-indexed; for staggered it aggregates cohorts at each relative @@ -298,7 +304,8 @@ conformal_calibrate <- function(Y, D, X = NULL, I, II, T.on, r.cv, eff, ## exchangeability across units, with stationarity across cohort onsets). rel.mat <- T.on[, id.tr, drop = FALSE] etimes <- sort(unique(rel.mat[obstr])) - mk_et <- function(mode) { + mk_et <- function(mode, a) { + c.sim <- .conformal_quantile(Mj, a) b <- matrix(NA_real_, length(etimes), 4L, dimnames = list(as.character(etimes), c("eff", "CI.lower", "CI.upper", "p.value"))) for (k in seq_along(etimes)) { @@ -310,21 +317,24 @@ conformal_calibrate <- function(Y, D, X = NULL, I, II, T.on, r.cv, eff, gp <- as.vector(Gtil[, ts, drop = FALSE]); gp <- gp[is.finite(gp)] if (!is.finite(tg) || !is.finite(sc.tr) || sc.tr <= 0 || length(gp) < 2L) next Qk <- if (mode == "simultaneous" && ispost) c.sim - else if (!is.null(Qpool) && ispost) Qpool - else .conformal_quantile(gp, alpha) + else if (ispost) qt_post(gp, a) + else .conformal_quantile(gp, a) half <- if (is.finite(Qk)) sc.tr * Qk else Inf b[k, ] <- c(tg, tg - half, tg + half, .conformal_pval(abs(tg) / sc.tr, gp)) } b } - band.et <- mk_et(band.type) - band.et.sim <- if (identical(band.type, "simultaneous")) band.et else mk_et("simultaneous") + a.inner <- min(2 * alpha, 0.5) # inner (1 - 2*alpha) band + band.et <- mk_et(band.type, alpha) + band.et.inner <- mk_et(band.type, a.inner) + band.et.sim <- if (identical(band.type, "simultaneous")) band.et else mk_et("simultaneous", alpha) list(att = m.tr, ci = c(ci[1], ci[2]), p.value = p.value, score.tr = s.tr, score.co = s.co.vec, status = status, n.calib = Ncal, n.eff = if (is.null(w.co)) Ncal else .conformal_neff(w.co), scale = scale, center = center, weight = weight, band.type = band.type, cutoff = cutoff, staggered = staggered, form = "jackknife+", - band = band, band.sim = band.sim, band.et = band.et, band.et.sim = band.et.sim, + band = band, band.sim = band.sim, band.et = band.et, + band.et.inner = band.et.inner, band.et.sim = band.et.sim, post.idx = post.idx, pre.idx = pre.idx) } diff --git a/tests/testthat/test-conformal.R b/tests/testthat/test-conformal.R index a43bdfdd..1c0f9d28 100644 --- a/tests/testthat/test-conformal.R +++ b/tests/testthat/test-conformal.R @@ -243,6 +243,9 @@ test_that("plot() works on conformal fits (block + staggered, gap + counterfactu w_pt <- mean(fb$est.att[post, "CI.upper"] - fb$est.att[post, "CI.lower"]) w_sim <- mean(fb$est.att.sim[post, "CI.upper"] - fb$est.att.sim[post, "CI.lower"]) expect_gte(w_sim, w_pt) + ## the inner (1 - 2*alpha) band (est.att90) is narrower than the main band + w_in <- mean(fb$est.att90[post, "CI.upper"] - fb$est.att90[post, "CI.lower"]) + expect_lt(w_in, w_pt) }) ## ---- staggered adoption ----------------------------------------------------- From f4f03e6b841a0f9d51446791d063d1e4305bd954 Mon Sep 17 00:00:00 2001 From: Yiqing Xu Date: Tue, 9 Jun 2026 08:50:21 -0700 Subject: [PATCH 12/14] docs(conformal): print inference line; NEWS + Rd for staggered/bands/inner - print.fect shows "Inference: conformal (scale = ..., band = ..., N.calib = ...)" for conformal fits. - NEWS 2.4.6 + man/fect.Rd note staggered-adoption support, event-study/ counterfactual plotting, and the true inner (1-2*alpha) est.att90 band. Co-Authored-By: Claude Opus 4.8 --- NEWS.md | 3 ++- R/print.R | 5 +++++ man/fect.Rd | 6 +++--- 3 files changed, 10 insertions(+), 4 deletions(-) diff --git a/NEWS.md b/NEWS.md index d3903397..233eb316 100644 --- a/NEWS.md +++ b/NEWS.md @@ -3,7 +3,8 @@ * New `vartype = "conformal"`: a cross-sectional conformal prediction interval for the average treatment effect on the treated. It ranks the treated unit's average post-treatment prediction error against the leave-one-control-out errors of the donors, giving a distribution-free interval with no bootstrap draws. The interval is closed-form and never empty. * Conformal options: `conformal.scale` (`"sd"` default, the studentized statistic; also `"none"`, `"rmspe"`, `"mad"`, `"diff"`), `conformal.center` (`"mean"`/`"median"`), `conformal.weight` (`"cell"`/`"unit"`/`"precision"` for multiple treated units), `conformal.band` (`"pointwise"`/`"simultaneous"`), and `conformal.cutoff` (`"per-period"`/`"pooled"`). The simultaneous (uniform) band is also stored in `fit$est.att.sim`. -* `vartype = "conformal"` requires a separated (controls-only) fit and `method` not in `c("mc", "both")`. It currently covers the core synthetic-control case; group, reversal, weighted (`W`), balanced-panel, placebo, and carryover options still use `bootstrap`/`jackknife`. +* Conformal supports both block (common-onset) and **staggered adoption** (per-cohort donor pools, union-window calibration, and event-time-aggregated bands). `plot()` renders the event-study and counterfactual for conformal fits, and `est.att90` is a true inner `1 - 2*alpha` band. +* `vartype = "conformal"` requires a separated (controls-only) fit and `method` not in `c("mc", "both")`. It currently covers the core synthetic-control case (block or staggered); group, reversal, weighted (`W`), balanced-panel, placebo, and carryover options still use `bootstrap`/`jackknife`. # fect 2.4.5 diff --git a/R/print.R b/R/print.R index 327f776e..7c6d4247 100644 --- a/R/print.R +++ b/R/print.R @@ -47,6 +47,11 @@ print.fect <- function(x, if (!is.null(x$cl.label)) { cat("Cluster SE: ", x$cl.label, "\n", sep = "") } + if (identical(x$vartype, "conformal") && !is.null(x$conformal)) { + cf <- x$conformal + cat("Inference: conformal (scale = ", cf$scale, ", band = ", cf$band.type, + ", N.calib = ", cf$n.calib, ")\n", sep = "") + } if (switch.on == TRUE) { if (!is.null(time.on.lim)) { diff --git a/man/fect.Rd b/man/fect.Rd index 4fbba516..176ff37e 100644 --- a/man/fect.Rd +++ b/man/fect.Rd @@ -99,9 +99,9 @@ In v2.3.1, \code{W.est} and \code{W.agg} (when both supplied) must point to the These three conditions correspond to Gates A, B, and C in the three-gate defense system (see \code{ARCHITECTURE.md}). \code{"conformal"} likewise requires a separated (controls-only) - fit and \code{method} not in \code{c("mc", "both")}; it currently - supports the core synthetic-control case (no group / reversal / weighted / - balanced / placebo / carryover options). See the \code{conformal.*} + fit and \code{method} not in \code{c("mc", "both")}; it supports both block + (common-onset) and staggered adoption, but not group / reversal / weighted / + balanced / placebo / carryover options. See the \code{conformal.*} arguments. For all other settings, \code{vartype = "bootstrap"} is recommended.} \item{para.error}{a string specifying the residual-error model used by the parametric bootstrap path; sub-option of \code{vartype = "parametric"} From 9b79579df6cf2550870ddf24b5387b4dae417de7 Mon Sep 17 00:00:00 2001 From: Yiqing Xu Date: Wed, 10 Jun 2026 00:52:48 -0700 Subject: [PATCH 13/14] fix(conformal): null-guard est.beta labelling with covariates; staggered+X test With p > 0 and vartype = 'conformal', fect returns no est.beta/est.marginal (no coefficient draws), and the Xname labelling in default.R crashed with 'attempt to set rownames on an object with no dimensions'. Guard on !is.null() so any vartype that omits inference objects skips labelling. Adds a staggered + covariate conformal regression test (asserts finite ATT, ordered CI, point beta present, est.beta absent). Also makes the test file runnable against the installed package: fect::fect in the new test, fect:::conformal_calibrate in the coverage helper. --- R/default.R | 6 ++++-- tests/testthat/test-conformal.R | 27 ++++++++++++++++++++++++++- 2 files changed, 30 insertions(+), 3 deletions(-) diff --git a/R/default.R b/R/default.R index 2b3e0f01..3c0b6101 100644 --- a/R/default.R +++ b/R/default.R @@ -3228,7 +3228,9 @@ fect.default <- function( if (binary == TRUE) { rownames(out$marginal) <- Xname.tmp } - if (se == TRUE) { + if (se == TRUE && !is.null(out$est.beta)) { + ## vartype = "conformal" returns no coefficient draws, hence no + ## est.beta/est.marginal; skip the labelling rather than crash. rownames(out$est.beta) <- Xname.tmp colnames(out$est.beta) <- c( "Coef", @@ -3237,7 +3239,7 @@ fect.default <- function( "CI.upper", "p.value" ) - if (binary == TRUE) { + if (binary == TRUE && !is.null(out$est.marginal)) { rownames(out$est.marginal) <- Xname.tmp } if (placeboTest == TRUE) { diff --git a/tests/testthat/test-conformal.R b/tests/testthat/test-conformal.R index 1c0f9d28..0e83e2b9 100644 --- a/tests/testthat/test-conformal.R +++ b/tests/testthat/test-conformal.R @@ -80,7 +80,7 @@ test_that("effective sample size is sensible", { fit <- suppressMessages(fect(Y ~ D, data = dat, index = c("id", "time"), method = "gsynth", force = 3, CV = FALSE, r = 2, se = FALSE, parallel = FALSE)) - cc <- conformal_calibrate(Y = Y, D = g$D, I = fit$I, II = fit$II, T.on = fit$T.on, + cc <- fect:::conformal_calibrate(Y = Y, D = g$D, I = fit$I, II = fit$II, T.on = fit$T.on, r.cv = fit$r.cv, eff = fit$eff, method = "gsynth", scale = scale, alpha = alpha) if (cc$status != "ok") bad <- bad + 1L @@ -301,3 +301,28 @@ test_that("simultaneous band gives materially better joint coverage than pointwi expect_gt(mean(js) - mean(jp), 0.15) # simultaneous materially better jointly expect_gt(mean(js), 0.65) # and well above pointwise's joint coverage }) + +test_that("conformal runs with covariates under staggered adoption (est.beta absent, no crash)", { + ## Regression guard: with p > 0 and vartype = "conformal", fect returns no + ## est.beta (no coefficient draws); the Xname labelling in default.R must + ## skip rather than crash ("attempt to set 'rownames' on an object with no + ## dimensions", caught 2026-06-10). + set.seed(424) + N <- 20; TT <- 16 + Fm <- matrix(rnorm(TT * 2), TT, 2) + L <- matrix(rnorm(N * 2), N, 2) + X1 <- matrix(rnorm(TT * N), TT, N) + Y <- Fm %*% t(L) + 1.5 * X1 + matrix(rnorm(TT * N), TT, N) + dat <- data.frame(id = rep(1:N, each = TT), time = rep(1:TT, N), + Y = as.vector(Y), X1 = as.vector(X1)) + dat$D <- as.integer((dat$id == 1 & dat$time >= 9) | + (dat$id == 2 & dat$time >= 13)) + f <- suppressMessages(suppressWarnings( + fect::fect(Y ~ D + X1, data = dat, index = c("id", "time"), method = "gsynth", + force = 3, CV = FALSE, r = 2, se = TRUE, vartype = "conformal", + parallel = FALSE))) + expect_true(is.finite(f$est.avg[1])) + expect_true(f$est.avg[3] < f$est.avg[4]) # CI ordered + expect_true(!is.null(f$beta)) # point coefficient still reported + expect_true(is.null(f$est.beta)) # and no inference object invented +}) From 00483a39d03289811618cfb92b4aedb6886c247f Mon Sep 17 00:00:00 2001 From: Yiqing Xu Date: Wed, 10 Jun 2026 07:41:07 -0700 Subject: [PATCH 14/14] feat(conformal): two-sided conformal.fit hook; time.component.from auto-default - conformal.fit is now user-facing (was dead code: defined in conformal_calibrate but never threaded through fect()). Threaded through all three signatures -> fect_boot -> conformal_calibrate, and made two-sided: each treated unit is scored through the same custom learner (own pre-window, T0 = onset - 1), not just the held-out donors. Validation: a quadprog simplex-SC learner reproduces the hand-rolled Proposition 99 conformal interval exactly (ATT -19.5136, CI95 [-53.6576, 14.6304], p .0769). - time.component.from defaults to NULL = auto: resolves to 'nevertreated' under vartype='conformal' with method fe/ife/cfe (strict separation, matching the paper's Algorithm 1; message emitted) and 'notyettreated' otherwise; explicit 'notyettreated' + conformal warns (weak separation). - Tests: custom-learner ATT equality (block), auto-resolution equivalence + warning capture (ife); the conformal test file now runs against the installed package. - Docs: fect.Rd (conformal.fit entry; tcf NULL semantics), NEWS. GATE NOTE: conformal test file fully green; full testthat suite re-run pending (two session-killed attempts, zero failures at both kill points); finish before merging PR #143. --- NEWS.md | 3 ++ R/boot.R | 4 ++- R/conformal.R | 16 +++++++++++ R/default.R | 29 +++++++++++++++++-- man/fect.Rd | 9 ++++-- tests/testthat/test-conformal.R | 51 +++++++++++++++++++++++++++++++++ 6 files changed, 105 insertions(+), 7 deletions(-) diff --git a/NEWS.md b/NEWS.md index 233eb316..b35830e0 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,6 +5,9 @@ * Conformal options: `conformal.scale` (`"sd"` default, the studentized statistic; also `"none"`, `"rmspe"`, `"mad"`, `"diff"`), `conformal.center` (`"mean"`/`"median"`), `conformal.weight` (`"cell"`/`"unit"`/`"precision"` for multiple treated units), `conformal.band` (`"pointwise"`/`"simultaneous"`), and `conformal.cutoff` (`"per-period"`/`"pooled"`). The simultaneous (uniform) band is also stored in `fit$est.att.sim`. * Conformal supports both block (common-onset) and **staggered adoption** (per-cohort donor pools, union-window calibration, and event-time-aggregated bands). `plot()` renders the event-study and counterfactual for conformal fits, and `est.att90` is a true inner `1 - 2*alpha` band. * `vartype = "conformal"` requires a separated (controls-only) fit and `method` not in `c("mc", "both")`. It currently covers the core synthetic-control case (block or staggered); group, reversal, weighted (`W`), balanced-panel, placebo, and carryover options still use `bootstrap`/`jackknife`. +* `conformal.fit`: a user-facing hook for custom separated learners under `vartype = "conformal"`. Supply `f(Y, X, time, control.ids, target.id, T0)` returning the imputed untreated path; both sides of the calibration --- every held-out control and each treated unit (fitted on its own pre-window) --- run through the same learner. With a simplex synthetic-control learner this reproduces the hand-rolled Proposition 99 conformal interval exactly. +* `time.component.from` now defaults to `NULL` (auto): under `vartype = "conformal"` with `method` `"fe"`/`"ife"`/`"cfe"` it resolves to `"nevertreated"` so the main fit is strictly separated like the calibration (a message is emitted); otherwise the legacy `"notyettreated"`. Explicit `"notyettreated"` with conformal is honored with a weak-separation warning. +* Fix: `vartype = "conformal"` with covariates crashed during output labelling (`est.beta` does not exist for conformal --- no coefficient draws); the labelling is now null-guarded. Staggered + covariate conformal fits are covered by a regression test. # fect 2.4.5 diff --git a/R/boot.R b/R/boot.R index 14254df6..85110fa6 100644 --- a/R/boot.R +++ b/R/boot.R @@ -147,6 +147,7 @@ fect_boot <- function( conformal.weight = "cell", conformal.band = "pointwise", conformal.cutoff = "per-period", + conformal.fit = NULL, quantile.CI = FALSE, nboots = 200, parallel = TRUE, @@ -576,7 +577,8 @@ fect_boot <- function( norm.para = norm.para, scale = conformal.scale, center = conformal.center, weight = conformal.weight, band.type = conformal.band, - cutoff = conformal.cutoff, alpha = alpha + cutoff = conformal.cutoff, alpha = alpha, + conformal.fit = conformal.fit ) ## back out a nominal S.E. from the symmetric conformal CI (display only; diff --git a/R/conformal.R b/R/conformal.R index 1db44c84..66ca9229 100644 --- a/R/conformal.R +++ b/R/conformal.R @@ -199,6 +199,22 @@ conformal_calibrate <- function(Y, D, X = NULL, I, II, T.on, r.cv, eff, g } + ## --- custom separated learner: score the TREATED units through the same + ## hook the donors use, so both sides of the calibration come from one + ## estimator. (Previously the treated gap came from the main fect fit while + ## donors went through conformal.fit -- mismatched arms.) Each treated unit + ## is fitted on its OWN pre-window (onset - 1); donors keep the union + ## pseudo-window in loo_gap, mirroring the built-in impute_Y0 path. + if (!is.null(conformal.fit)) { + for (i in seq_len(Ntr)) { + id <- id.tr[i] + T0i <- if (is.finite(onsets[i])) onsets[i] - 1L else max(pre.idx) + y0 <- conformal.fit(Y = Y, X = X, time = seq_len(TT), + control.ids = valid.co, target.id = id, T0 = T0i) + eff[, id] <- Y[, id] - y0 + } + } + ## --- full per-control LOO gap paths (calendar-indexed), and a per-unit scale ## from each control's pre-period gaps. Keeping the whole path lets us serve ## both the scalar average interval AND the per-period band from one pass. diff --git a/R/default.R b/R/default.R index 3c0b6101..992b6f2d 100644 --- a/R/default.R +++ b/R/default.R @@ -39,7 +39,7 @@ fect <- function( na.rm = FALSE, # remove missing values index, # c(unit, time) indicators force = "two-way", # fixed effects demeaning - time.component.from = "notyettreated", # factor estimation sample: "notyettreated" or "nevertreated" + time.component.from = NULL, # factor estimation sample: NULL = auto ("nevertreated" under vartype="conformal" for fe/ife/cfe, else "notyettreated") em = TRUE, # EM algorithm for missing data; FALSE uses direct SVD (requires complete estimation sample) r = 0, # number of factors lambda = NULL, # mc method: regularization parameter @@ -63,6 +63,7 @@ fect <- function( conformal.weight = "cell", # vartype="conformal": multi-treated aggregation cell|unit|precision conformal.band = "pointwise", # vartype="conformal": est.att band pointwise|simultaneous conformal.cutoff = "per-period", # vartype="conformal": pointwise cutoff per-period|pooled + conformal.fit = NULL, # vartype="conformal": custom separated learner f(Y, X, time, control.ids, target.id, T0) -> y0 path cl = NULL, ci.method = "normal", # CI method for fect's est.* slots: "normal" (Wald: theta_hat +- z * SE) or "basic" (reflected pivot, Davison-Hinkley 1997 Sec. 5.2.1). For percentile / bc / bca on alternative estimands (att.cumu, aptt, log.att), call estimand(fit, type, ci.method) post-fit quantile.CI = NULL, # DEPRECATED: use ci.method instead. NULL sentinel = "not supplied"; legacy FALSE -> ci.method = "normal", legacy TRUE -> ci.method = "basic" @@ -127,7 +128,7 @@ fect.formula <- function( na.rm = FALSE, # remove missing values index, # c(unit, time) indicators force = "two-way", # fixed effects demeaning - time.component.from = "notyettreated", # factor estimation sample: "notyettreated" or "nevertreated" + time.component.from = NULL, # factor estimation sample: NULL = auto ("nevertreated" under vartype="conformal" for fe/ife/cfe, else "notyettreated") em = TRUE, # EM algorithm for missing data; FALSE uses direct SVD (requires complete estimation sample) r = 0, # nubmer of factors lambda = NULL, # mc method: regularization parameter @@ -151,6 +152,7 @@ fect.formula <- function( conformal.weight = "cell", # vartype="conformal": multi-treated aggregation cell|unit|precision conformal.band = "pointwise", # vartype="conformal": est.att band pointwise|simultaneous conformal.cutoff = "per-period", # vartype="conformal": pointwise cutoff per-period|pooled + conformal.fit = NULL, # vartype="conformal": custom separated learner f(Y, X, time, control.ids, target.id, T0) -> y0 path cl = NULL, ci.method = "normal", # CI method for fect's est.* slots: "normal" or "basic" quantile.CI = NULL, # DEPRECATED: use ci.method instead @@ -271,6 +273,7 @@ fect.formula <- function( conformal.weight = conformal.weight, conformal.band = conformal.band, conformal.cutoff = conformal.cutoff, + conformal.fit = conformal.fit, cl = cl, ci.method = ci.method, quantile.CI = quantile.CI, @@ -337,7 +340,7 @@ fect.default <- function( na.rm = FALSE, # remove missing values index, # c(unit, time) indicators force = "two-way", # fixed effects demeaning - time.component.from = "notyettreated", # factor estimation sample: "notyettreated" or "nevertreated" + time.component.from = NULL, # factor estimation sample: NULL = auto ("nevertreated" under vartype="conformal" for fe/ife/cfe, else "notyettreated") em = TRUE, # EM algorithm for missing data; FALSE uses direct SVD (requires complete estimation sample) r = 0, # nubmer of factors lambda = NULL, ## mc method: regularization parameter @@ -361,6 +364,7 @@ fect.default <- function( conformal.weight = "cell", # vartype="conformal": multi-treated aggregation cell|unit|precision conformal.band = "pointwise", # vartype="conformal": est.att band pointwise|simultaneous conformal.cutoff = "per-period", # vartype="conformal": pointwise cutoff per-period|pooled + conformal.fit = NULL, # vartype="conformal": custom separated learner f(Y, X, time, control.ids, target.id, T0) -> y0 path cl = NULL, ci.method = "normal", # CI method for fect's est.* slots: "normal" or "basic" quantile.CI = NULL, # DEPRECATED: use ci.method instead @@ -740,6 +744,24 @@ fect.default <- function( ) } + ## resolve time.component.from (NULL = auto). Under vartype = "conformal" + ## the calibration is strictly separated (controls-only); resolving the + ## MAIN fit to "nevertreated" as well makes the treated unit's gap obey + ## strict separation too, matching the paper's Algorithm 1 exactly. + if (is.null(time.component.from)) { + if (vartype == "conformal" && method %in% c("fe", "ife", "cfe")) { + time.component.from <- "nevertreated" + message("vartype = \"conformal\": time.component.from set to \"nevertreated\" (strict separation; pass \"notyettreated\" explicitly to override).") + } else { + time.component.from <- "notyettreated" + } + } else if (vartype == "conformal" && identical(time.component.from, "notyettreated") && + method %in% c("fe", "ife", "cfe")) { + warning("vartype = \"conformal\" with time.component.from = \"notyettreated\": ", + "the main fit pools treated pre-treatment cells (weak separation only); ", + "\"nevertreated\" gives the exact finite-sample guarantee.", call. = FALSE) + } + ## validate time.component.from if (!time.component.from %in% c("notyettreated", "nevertreated")) { stop("\"time.component.from\" must be \"notyettreated\" or \"nevertreated\".") @@ -2796,6 +2818,7 @@ fect.default <- function( conformal.weight = conformal.weight, conformal.band = conformal.band, conformal.cutoff = conformal.cutoff, + conformal.fit = conformal.fit, quantile.CI = .quantile.CI.bool, nboots = nboots, parallel = parallel, diff --git a/man/fect.Rd b/man/fect.Rd index 176ff37e..255a432b 100644 --- a/man/fect.Rd +++ b/man/fect.Rd @@ -9,7 +9,7 @@ assumptions.} group = NULL, na.rm = FALSE, index, force = "two-way", - time.component.from = "notyettreated", em = TRUE, + time.component.from = NULL, em = TRUE, r = 0, lambda = NULL, nlambda = 10, CV = NULL, k = 20, cv.prop = 0.1, cv.method = "rolling", cv.nobs = 3, cv.donut = 1, cv.buffer = 1, criterion = "mspe", @@ -18,7 +18,8 @@ assumptions.} para.error = "auto", conformal.scale = "sd", conformal.center = "mean", conformal.weight = "cell", conformal.band = "pointwise", - conformal.cutoff = "per-period", cl = NULL, + conformal.cutoff = "per-period", + conformal.fit = NULL, cl = NULL, ci.method = "normal", quantile.CI = NULL, nboots = 200, alpha = 0.05, parallel = TRUE, cores = NULL, tol = 1e-5, @@ -62,7 +63,7 @@ In v2.3.1, \code{W.est} and \code{W.agg} (when both supplied) must point to the \item{na.rm}{a logical flag indicating whether to list-wise delete missing observations. Default to FALSE. If \code{na.rm = FALSE}, it allows the situation when Y is missing but D is not missing for some observations. If \code{na.rm = TRUE}, it will list-wise delete observations whose Y, D, or X is missing.} \item{index}{a character vector specifying the unit (first element) and time (second element) indicators. For most methods, must be of length 2. For \code{method = "cfe"}, additional elements (third, fourth, etc.) specify extra fixed-effect grouping variables. Every observation should be uniquely defined by the pair of the unit and time indicator.} \item{force}{a string indicating whether unit or time or both fixed effects will be imposed. Must be one of the following, "none", "unit", "time", or "two-way". The default is "two-way".} -\item{time.component.from}{Controls which units provide the time-varying model components (time fixed effects, factor structure, temporal dynamics). Options are \code{"notyettreated"} (default) --- all units contribute during their pre-treatment periods, or \code{"nevertreated"} --- only never-treated units estimate the time components, which are then projected onto treated units.} +\item{time.component.from}{Controls which units provide the time-varying model components (time fixed effects, factor structure, temporal dynamics). Options are \code{"notyettreated"} --- all units contribute during their pre-treatment periods, or \code{"nevertreated"} --- only never-treated units estimate the time components, which are then projected onto treated units. The default \code{NULL} resolves automatically: \code{"nevertreated"} under \code{vartype = "conformal"} with \code{method} \code{"fe"}, \code{"ife"}, or \code{"cfe"} (strict separation, matching the conformal calibration; a message is emitted), and \code{"notyettreated"} otherwise. Passing \code{"notyettreated"} explicitly together with \code{vartype = "conformal"} is honored with a warning (the main fit then pools treated pre-treatment cells).} \item{em}{a logical flag indicating whether to use the EM algorithm for missing data in the estimation sample. Default is \code{TRUE}. Setting \code{em = FALSE} requires a complete estimation sample and is only compatible with \code{time.component.from = "nevertreated"}.} \item{r}{an integer specifying the number of factors. If \code{CV = TRUE}, the cross validation procedure will select the optimal number of factors from \code{r} to 5.} \item{lambda}{a single or sequence of positive numbers specifying the hyper-parameter sequence for matrix completion method. If \code{lambda} is a sequence and \code{CV = 1}, cross-validation will be performed.} @@ -139,6 +140,8 @@ In v2.3.1, \code{W.est} and \code{W.agg} (when both supplied) must point to the \item{conformal.cutoff}{for \code{vartype = "conformal"} pointwise bands: whether the rank cutoff is computed \code{"per-period"} (separately at each period) or \code{"pooled"} (one cutoff over all post-period control gaps).} + +\item{conformal.fit}{for \code{vartype = "conformal"}: an optional custom separated learner, a function \code{f(Y, X, time, control.ids, target.id, T0)} returning the imputed untreated path for \code{target.id}, fitted on \code{control.ids} (and the target's own cells up to \code{T0} only). When supplied, both sides of the calibration run through it: every held-out control and each treated unit (the latter on its own pre-treatment window). The learner must be separated in the sense of the accompanying paper; built-in methods are used when \code{NULL} (default).} \item{cl}{a string specifying the cluster column for cluster bootstrapping. When \code{group.fe} is set to a single column and \code{cl} is unset, \code{cl} auto-defaults to \code{group.fe[1]} --- the natural choice when diff --git a/tests/testthat/test-conformal.R b/tests/testthat/test-conformal.R index 0e83e2b9..5e9cbe45 100644 --- a/tests/testthat/test-conformal.R +++ b/tests/testthat/test-conformal.R @@ -326,3 +326,54 @@ test_that("conformal runs with covariates under staggered adoption (est.beta abs expect_true(!is.null(f$beta)) # point coefficient still reported expect_true(is.null(f$est.beta)) # and no inference object invented }) + +test_that("conformal.fit scores the treated unit too (custom two-step DiD learner)", { + ## Regression guard (2026-06-10): previously the hook scored donors only and + ## the treated gap came from the main fect fit -- mismatched arms. Now both + ## sides run through the hook: with a deterministic two-step DiD learner the + ## ATT the interval is built on must equal the hand-computed DiD gap. + set.seed(77) + N <- 18; TT <- 14 + Y <- outer(rnorm(TT), rep(1, N)) + outer(rep(1, TT), rnorm(N)) + + matrix(rnorm(TT * N), TT, N) + dat <- data.frame(id = rep(1:N, each = TT), time = rep(1:TT, N), + Y = as.vector(Y)) + dat$D <- as.integer(dat$id == 1 & dat$time >= 10) + did_fit <- function(Y, X, time, control.ids, target.id, T0) { + trend <- rowMeans(Y[, control.ids, drop = FALSE]) + trend - mean(trend[seq_len(T0)]) + mean(Y[seq_len(T0), target.id]) + } + f <- suppressMessages(fect::fect(Y ~ D, data = dat, index = c("id", "time"), + method = "gsynth", force = 3, CV = FALSE, r = 0, se = TRUE, + vartype = "conformal", conformal.fit = did_fit, parallel = FALSE)) + y0 <- did_fit(Y, NULL, 1:TT, 2:N, 1, 9) + att <- mean((Y[, 1] - y0)[10:14]) + expect_equal(unname(f$est.avg[1]), att, tolerance = 1e-8) + expect_true(f$est.avg[3] < f$est.avg[4]) +}) + +test_that("time.component.from auto-resolves under vartype = 'conformal' (ife)", { + set.seed(88) + N <- 16; TT <- 14 + Fm <- matrix(rnorm(TT), TT, 1); L <- matrix(rnorm(N), N, 1) + Y <- Fm %*% t(L) + matrix(rnorm(TT * N), TT, N) + dat <- data.frame(id = rep(1:N, each = TT), time = rep(1:TT, N), + Y = as.vector(Y)) + dat$D <- as.integer(dat$id == 1 & dat$time >= 10) + run <- function(...) fect::fect(Y ~ D, data = dat, index = c("id", "time"), + method = "ife", force = 3, CV = FALSE, r = 1, se = TRUE, + vartype = "conformal", parallel = FALSE, ...) + expect_message(f1 <- suppressWarnings(run()), "nevertreated") + f2 <- suppressMessages(run(time.component.from = "nevertreated")) + expect_equal(f1$est.avg, f2$est.avg) + w <- testthat::capture_warnings( + f3 <- suppressMessages(run(time.component.from = "notyettreated"))) + expect_true(any(grepl("weak separation", w))) + ## non-conformal callers keep the legacy default + g1 <- suppressMessages(fect::fect(Y ~ D, data = dat, index = c("id", "time"), + method = "ife", force = 3, CV = FALSE, r = 1, se = FALSE, parallel = FALSE)) + g2 <- suppressMessages(fect::fect(Y ~ D, data = dat, index = c("id", "time"), + method = "ife", force = 3, CV = FALSE, r = 1, se = FALSE, parallel = FALSE, + time.component.from = "notyettreated")) + expect_equal(g1$att.avg, g2$att.avg) +})