---
title: "Simultaneous confidence bands and the zero-effect test"
author: "Gregory Faletto"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
    math_method: "mathml"
vignette: >
  %\VignetteIndexEntry{Simultaneous confidence bands and the zero-effect test}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
library(fetwfe)
```

This vignette covers the package's **family-wise inference** tools for `fetwfe()`, `betwfe()`, `etwfe()`, and `twfeCovs()`: `simultaneousCIs()`, which returns confidence bands that hold *jointly* across a family of effects (not merely one effect at a time), and the **selection-consistency zero-effect test** --- the sense in which `fetwfe()`'s `selected` flag is an implicit test of whether a cohort's effect is truly zero.

It is a companion to `vignette("inference_vignette")` (how the package's standard errors are computed) and `vignette("high_dimensional_vignette")` (debiased inference when $p \geq NT$). We reuse that first vignette's small running example:

```{r setup-fit}
set.seed(2026)
sim_coefs <- genCoefs(G = 3, T = 6, d = 2, density = 0.5, eff_size = 2)
sim_data <- simulateData(
  sim_coefs,
  N = 120,
  sig_eps_sq = 1,
  sig_eps_c_sq = 0.5,
  seed = NA
)

res_default <- fetwfe(
  pdata = sim_data$pdata,
  time_var = sim_data$time_var,
  unit_var = sim_data$unit_var,
  treatment = sim_data$treatment,
  covs = sim_data$covs,
  response = sim_data$response,
  sig_eps_sq = sim_data$sig_eps_sq,
  sig_eps_c_sq = sim_data$sig_eps_c_sq
)
```

# 1. Simultaneous confidence intervals via `simultaneousCIs()`

The package's per-point standard errors and Wald intervals deliver `(1 - alpha)` coverage *per effect*: for any single cohort ATT, or any single event-time coefficient, the CI contains the true value with probability approximately `1 - alpha`. But when you make a claim involving multiple effects jointly --- for example, "no event-time coefficient is significantly positive at the 5% level" --- per-effect coverage does not imply family-wise coverage. The probability that *every* coefficient in a family of `K` effects lies in its per-point CI can be far below `1 - alpha` when `K` is large.

`simultaneousCIs()` returns intervals that maintain family-wise coverage at the requested `1 - alpha`, using the exact joint multivariate-normal critical value $c_{1-\alpha}$ defined as the $(1 - \alpha)$ quantile of $\max_k |Z_k|$ where $Z \sim N(0, \rho)$ and $\rho$ is the correlation matrix of the effect family's joint covariance. Under Theorem (c$'$) (`paper_arxiv.tex:1233`) and Assumption $(\Psi\text{-IF})$ (assumption equation `paper_arxiv.tex:2013`; in-prose discussion of when it implies tight Gaussianity at paper line 1268), this joint Gaussianity is exact asymptotically; under the paper's fixed-dim framing (AE 1(d)), no high-dimensional correction is needed.

For an event-study family on the fitted `res_default` object from the setup above:

```{r}
sci <- simultaneousCIs(res_default, family = "event_study", alpha = 0.05)
print(sci)
```

The function returns a `simultaneous_cis` S3 object with `$ci` (the interval data frame), `$critical_value` (the simultaneous $c_{1-\alpha}$), `$pointwise_critical_value` ($z_{1-\alpha/2}$, for reference), and `$bonferroni_critical_value` ($z_{1-\alpha/(2K)}$, the conservative family-wise upper bound). When effects are positively correlated --- as they typically are in DiD, because event-study coefficients share the regression-coefficient variance piece --- the simultaneous critical value sits strictly between the pointwise and Bonferroni values, giving tighter intervals than Bonferroni while preserving family-wise coverage.

The four available family options are `"event_study"` (one effect per post-treatment event time), `"cohort"` (one effect per treated cohort), `"all_post_treatment"` (one effect per `(g, t)` cell), and `"custom"` (a user-supplied `K x num_treats` contrast matrix following the `multcomp::glht()` convention; see `?simultaneousCIs`).

The critical value is computed via `mvtnorm::qmvnorm(..., algorithm = mvtnorm::GenzBretz())` (mvtnorm's default quasi-Monte Carlo integrator), sub-second on a typical fit for any realistic FETWFE family (no `K` cap of practical concern). The function is deterministic in its inputs: the same fit plus the same family plus the same `alpha` plus the same `contrasts` produces the same critical value across calls. This is achieved by save/restoring the caller's `.Random.seed` and using a fixed internal seed immediately before the `qmvnorm()` call; the function does not mutate the caller's RNG state.

As of version 1.16.0 the `mvtnorm` package is an `Imports` dependency (it ships with the package), because simultaneous bands are now computed by default at fit time --- see the `ci_type` subsection below. (In versions 1.15.x `mvtnorm` was in `Suggests` and only loaded when `simultaneousCIs()` was called.) The `plot.simultaneous_cis` method additionally requires `ggplot2`, which is in `Suggests`; install it with `install.packages("ggplot2")` if you have not already.

The default analytic path matches `did::aggte(cband = TRUE)` from the Callaway-Sant'Anna `did` package, which computes simultaneous bands via the CCK (2013) multiplier bootstrap. Asymptotically the two critical values are identical under fixed-dim Gaussianity; the analytic `simultaneousCIs()` path is deterministic and sub-second, whereas a multiplier bootstrap is stochastic and resampling-based.

## The multiplier-bootstrap method: `method = "bootstrap"`

`simultaneousCIs(..., method = "bootstrap")` computes the same sup-t critical value by the CCK (2013) multiplier bootstrap directly --- the procedure `did` uses. It perturbs the per-unit influence functions $\hat f_i^{(k)}$ with random multipliers $\xi_i$ (`multiplier = "rademacher"` or `"mammen"`), forms $T^{(b)} = \max_k |\frac{1}{N}\sum_i \xi_i \hat f_i^{(k)}| / \widehat{\mathrm{sd}}_k$ over $B$ replicates, and takes the $(1-\alpha)$ quantile. The analytic and bootstrap critical values agree up to Monte-Carlo error, and both sit between the pointwise and Bonferroni values:

```{r}
sci_boot <- simultaneousCIs(res_default, family = "all_post_treatment",
                            method = "bootstrap", B = 1000, seed = 1)
c(pointwise  = sci_boot$pointwise_critical_value,
  analytic   = simultaneousCIs(res_default, family = "all_post_treatment")$critical_value,
  bootstrap  = sci_boot$critical_value,
  bonferroni = sci_boot$bonferroni_critical_value)
```

Two reasons to prefer the bootstrap. First, it **scales to large effect families**: `qmvnorm` integration strains at hundreds of effects (e.g. simultaneous bands over the full disaggregated cohort $\times$ event-time grid), whereas the bootstrap stays stable and linear in $B \cdot K$. Second, it is **heteroskedasticity/cluster-robust by construction** --- its per-unit influence matrix $F$ reproduces the cluster-robust joint covariance ($\mathrm{crossprod}(F)/(NT)^2$ equals it up to the $N/(N-1)$ finite-sample factor), so on a `se_type = "cluster"` fit the bootstrap per-effect standard errors equal the analytic ones exactly (and on a default fit the bootstrap reports the cluster-robust band regardless). Pass an integer `seed` for reproducible draws; with `seed = NULL` the draws come from the ambient RNG.

`method = "bootstrap"` supports all four families. For `family = "event_study"` the bootstrap additionally perturbs a **second, independent** multiplier stream over a per-unit cohort-probability (propensity) influence function, alongside the regression-channel stream --- a two-channel bootstrap reproducing the analytic variance decomposition $\Sigma = \Sigma_1 + \Sigma_2$. The event-study family is the one whose propensity term $\Sigma_2$ is non-zero, because each event-time effect pools across the cohorts treated by that event time, weighted by the estimated cohort probabilities $\hat\pi_g = N_g / N$; its per-unit propensity influence function is the cohort-sample-proportion summand $\xi_i = e_{W_i} - \hat\pi$ (paper assumption $(\Psi\text{-IF})$). The two channels are asymptotically uncorrelated (paper Theorem C.1, Step 5), so their variances add with no cross-term --- which is exactly why the multipliers are drawn independently. On a `se_type = "cluster"` fit the event-study bootstrap per-effect standard errors equal the analytic ones to machine precision, and the critical value matches the analytic `qmvnorm` one up to Monte-Carlo error:

```{r}
es_boot <- simultaneousCIs(res_default, family = "event_study",
                           method = "bootstrap", B = 1000, seed = 1)
c(pointwise  = es_boot$pointwise_critical_value,
  analytic   = simultaneousCIs(res_default, family = "event_study")$critical_value,
  bootstrap  = es_boot$critical_value,
  bonferroni = es_boot$bonferroni_critical_value)
```

**High-dimensional $p \geq NT$.** When the (full) design is high-dimensional --- where the analytic Gram inverse need not exist, so there is no analytic sup-t band --- the bootstrap automatically switches to the full-design **desparsified** construction of `debiasedATT()` (the high-dimensional FETWFE theory): each effect's debiasing direction becomes a nodewise (`riesz_lasso`) relaxed inverse over the singular Gram, with the bridge residuals. This is *uniformly valid* (not post-selection) and is the only route to simultaneous bands in that regime. It is validated in simulation (`fetwfe()` fits only): the scalar overall-ATT coverage (Theorem `debiased.highdim.thm`) and the family-wise *band* coverage (Theorem `debiased.highdim.joint.thm`, ≈0.92 over the event-study family at the $p > NT$ anchor of Faletto 2025) are both near-nominal, so the returned object carries a `regime` field and per-effect `feasibility` / `converged` / `lambda_node` diagnostics, and `lambda_c` tunes the nodewise penalty `lambda_node = lambda_c * max(|a|) * sqrt(log(p) / N)` --- too-small values leave the directions infeasible and trigger a warning. `lambda_c` accepts either a fixed positive number (default the theory scale `1.0`) or the string `"cv"`, which selects the constant **per fit** by cross-validating the Riesz loss over the KKT-feasible region (the desparsified-lasso / auto-DML standard), falling back to `1.0` when no grid penalty is feasible; one constant then serves both the `debiasedATT()` point estimate and every bootstrap band effect (issue #295). The default stays the fixed `1.0`; the overall-ATT coverage validated in Faletto (2025) is robust across the penalty values studied (a CV selector below the theory scale and a fixed constant above it, bracketing the default `1.0`). The simulator can generate the `p > NT` panels this needs (see the `simulateData()` documentation). `family = "event_study"` uses this desparsified path too in the high-dimensional regime, additionally carrying the per-unit propensity channel (the two-channel $\Sigma_1 + \Sigma_2$ bootstrap). In the high-dimensional regime the band is **centered on the debiased estimate** (the high-dimensional FETWFE theory correction, matching `debiasedATT()`) rather than the post-selection bridge estimate, which would otherwise carry the selection bias. The desparsified path is `fetwfe()`-only; a non-`fetwfe()` $p \geq NT$ fit (e.g. `betwfe()`) has no debiased band, so it falls back to the fixed-$p$ selected-support band (rather than erroring) but emits a `warning()`. A coverage study (issue #308) shows that fallback band substantially **under-covers** in the $p \geq NT$ regime --- even when the selected support is low-dimensional --- because the bridge shrinks the band center toward zero (so it is biased away from the truth) and the post-selection standard error understates the sampling variability; treat such a band as unreliable and use a `fetwfe()` fit for valid high-dimensional simultaneous bands. See `vignette("high_dimensional_vignette")` for `debiasedATT()` and the accompanying high-dimensional diagnostics.

## Simultaneous bands are the default: the `ci_type` argument

As of version 1.16.0, the confidence-interval bounds that `fetwfe()` / `etwfe()` / `betwfe()` / `twfeCovs()` *return and display* are simultaneous by default --- matching `did`'s `cband = TRUE` default. A new fit-time argument `ci_type = c("simultaneous", "pointwise")` controls this:

* `ci_type = "simultaneous"` (the default) makes the returned `catt_df$ci_low` / `catt_df$ci_high` the **cohort-family** simultaneous bounds, and `eventStudy()`'s `ci_low` / `ci_high` the **event-study-family** simultaneous bounds --- the very bounds `simultaneousCIs(fit, family = "cohort")` / `simultaneousCIs(fit, family = "event_study")` return. `print()`, `summary()`, and `plot()` display these (labeled `[simultaneous 95% CI]`), and `broom::tidy()` (on the fit, on `eventStudy()`, and on `cohortStudy()`) surfaces them too.
* `ci_type = "pointwise"` restores the pre-1.16.0 per-effect Wald intervals everywhere (the bounds you would get from `estimate ± qnorm(1 - alpha/2) * se`).

The interval bounds **and** the per-cohort p-values (`p_value`) follow `ci_type`. Under `"simultaneous"` the `p_value` is the single-step max-T multiplicity-adjusted (family-wise) p-value matching the band; under `"pointwise"` it is the per-effect Wald p-value. The two correspond exactly: a coefficient lies outside its `1 - alpha` simultaneous band iff its adjusted `p_value` is `< alpha` (version 1.18.0, #200). The standard errors (`se`), selection flags (`selected`), and the overall-ATT confidence interval and p-value (a single scalar, `K = 1`, with no family and hence no widening) are **identical** under both settings --- for the bounds the only difference is the critical-value multiplier (`qnorm(1 - alpha/2)` for pointwise vs the joint $c_{1-\alpha}$ for simultaneous).

```{r}
fit_sim <- fetwfeWithSimulatedData(sim_data) # ci_type = "simultaneous" (default)
fit_pw <- fetwfeWithSimulatedData(sim_data, ci_type = "pointwise")

# Same point estimates and SEs; only the bounds differ (simultaneous wider).
cbind(
  estimate = fit_sim$catt_df$estimate,
  se = fit_sim$catt_df$se,
  simul_width = fit_sim$catt_df$ci_high - fit_sim$catt_df$ci_low,
  pointwise_width = fit_pw$catt_df$ci_high - fit_pw$catt_df$ci_low
)
```

Because the tidiers pass the object's bounds straight through, every CI surface in the package agrees under the default: `broom::tidy(fit)`'s cohort rows, `summary(fit)`, `cohortStudy(fit)`, and `plot(fit, type = "catt")` all report the same simultaneous bounds. (The overall-ATT row of `broom::tidy(fit)` stays the scalar pointwise interval.) When standard errors are unavailable (`q >= 1` for the bridge estimators, or a rank-deficient design) the bounds are `NA` under both settings, and the fit still succeeds --- the simultaneous path degrades cleanly to the (`NA`) pointwise bounds rather than throwing an error.

# 2. Selection consistency as an implicit zero-effect test

A question applied users ask constantly is: *is this cohort's treatment effect distinguishable from zero --- and can I conclude that it actually **is** zero?* For `fetwfe()` (and `betwfe()`), the answer is built into the estimator. The bridge penalty's *restriction selection consistency* (Faletto 2025, Theorem 6.2) makes the `selected` flag on `catt_df` an **implicit hypothesis test** of $H_0: \tau_{\text{ATT}}(g) = 0$, for free, with no extra computation. The mapping is:

* `selected = TRUE` (estimate $\neq 0$) is the asymptotic analogue of **reject $H_0$**: the effect is nonzero.
* `selected = FALSE` (estimate exactly $0$) is the asymptotic analogue of **fail to reject $H_0$** --- and, under Theorem 6.2, it is the *stronger* statement that the truth is zero with probability tending to one.

The conclusion for a selected-out cohort is delivered by the `selected` flag, **not** by a p-value. A cohort the penalty zeros out has `estimate == 0`, `se == 0`, and `p_value == NA` by construction --- there is no Wald p-value to read, and the inferential content lives entirely in `selected`. (A condensed, first-time-user-oriented version of this discussion appears under "Testing the zero-effect null" in the main package vignette; the treatment here is the deeper, inference-focused one.)

## Why this is stronger than a confidence interval on an unregularized estimator

For an unregularized estimator --- OLS, or plain `etwfe()` --- the coefficient estimate $\hat\beta$ is a continuous random variable, so $\hat\beta = 0$ *exactly* is a **measure-zero event**: it essentially never happens, and would carry no information about the truth if it did. The only zero-effect inference such an estimator offers is "the confidence interval contains $0$," i.e. *fail to reject* --- never *accept*. Selection consistency changes the picture qualitatively. Under Theorem 6.2, $\hat\beta_j = 0$ exactly with probability $\to 1$ when the truth is zero, and $\hat\beta_j \neq 0$ with probability $\to 1$ when the truth is nonzero, so observing $\hat\beta_j = 0$ *is* informative. The paper states the consequence directly: if $\tau_{\text{ATT}}(g,t) = 0$ then $\lim_{N\to\infty} P(\hat\tau_{\text{ATT}}(g,t) = 0) = 1$, which it calls "stronger than convergence in probability to $0$" (Faletto 2025, around the statement of Theorem 6.2).

It is worth being precise about *which* part of the theory delivers the test, because the paper is careful here. A conventional Wald statistic **cannot** test the zero null in this setting: when the relevant effect is truly zero the statistic is not asymptotically normal (it converges in probability to $0$ even after scaling by $\sqrt{N}$), so there is no valid Wald test of $H_0: \tau = 0$ to construct. That is *precisely why* the right tool is selection consistency and the `selected` flag, rather than a Wald statistic on the point estimate. The complementary half of the theory backs the effects the estimator *retains*: when an effect is truly nonzero, FETWFE is $\sqrt{N}$-consistent and asymptotically Gaussian on the selected support (the oracle property, Theorem 6.3), with the practical Wald confidence intervals and p-values of Theorem 6.4 --- the same machinery documented in `vignette("inference_vignette")`. So `selected = TRUE` cohorts come with the usual valid CIs and finite p-values, while `selected = FALSE` cohorts carry the stronger, selection-based zero-effect conclusion.

## The `q < 1` caveat

This interpretation rests on selection consistency, which the paper proves **only for the bridge exponent $q \in (0, 1)$**. The package default $q = 0.5$ satisfies it. For $q = 1$ (the lasso) selection consistency of FETWFE is *not* proven, so the implicit-test framing does not apply; for $q = 2$ (ridge) the solver never zeros coefficients at all, so `selected` is trivially `TRUE` everywhere and conveys nothing. Read the `selected` flag as a zero-effect test only under the default (or any $q < 1$) bridge penalty.

## A controlled-truth worked example

The cleanest way to see the implicit test is on simulated data where the truth is known. We draw coefficients in which the *middle* treated cohort has a true average treatment effect of exactly zero while its two neighbors are nonzero, then check whether FETWFE selects out exactly the zero cohort. `getTes()` reports the known per-cohort truth for comparison.

```{r implicit-test-demo, message = FALSE}
# Known truth: cohort 3 (the middle treated cohort) has a true effect of
# exactly zero; cohorts 2 and 4 are nonzero.
coefs <- genCoefs(G = 3, T = 6, d = 2, density = 0.2, eff_size = 2, seed = 3)
getTes(coefs)$actual_cohort_tes

dat <- simulateData(coefs, N = 120, sig_eps_sq = 1, sig_eps_c_sq = 0.5, seed = 3)
res <- fetwfeWithSimulatedData(dat, verbose = FALSE)

# The implicit test: the `selected` and `p_value` columns of catt_df.
res$catt_df[, c("cohort", "estimate", "se", "p_value", "selected")]
```

The true effects are $(-1.6, 0, 1.333)$: cohort 3's effect is exactly zero. FETWFE selects out precisely that cohort --- cohort 3 returns `selected = FALSE` with an estimate of exactly $0$ and `p_value = NA` --- while it recovers the two truly-nonzero cohorts with `selected = TRUE` and finite standard errors and p-values. Reading the `selected` column through the test mapping above: cohort 3 is *"fail to reject (and conclude zero),"* cohorts 2 and 4 are *"reject."* The estimator performed the zero-effect test automatically, and its verdict matches the known truth.

## Scope: this is the *post-treatment* zero-effect test

The package currently estimates only post-treatment treatment effects (event times $e \geq 0$), so the demonstration above is specifically the **post-treatment** zero-effect test, surfaced through `catt_df`. The same selection-consistency logic would support two further applications once the corresponding extended regression specifications exist: (a) selecting out pre-treatment placebo coefficients would give a **parallel-trends test**, and (b) selecting out heterogeneous-trend augmenting parameters would give a **homogeneous-trends test**. The package does not estimate those terms today, so this vignette makes no claim to perform either test; they are noted only as the natural extensions the same `selected`-as-implicit-test reasoning will cover when the specifications land.

# References

Faletto, G. (2025). "Fused Extended Two-Way Fixed Effects for Difference-in-Differences with Staggered Adoptions." *arXiv preprint* [arXiv:2312.05985](https://arxiv.org/abs/2312.05985).
