Skip to contents

This vignette documents how fetwfe(), betwfe(), etwfe(), and twfeCovs() compute their standard errors: the tight Gaussian variance used by default under (Ψ-IF), the two variance components V1+V2V_1 + V_2 the package exposes, the three values se_type can take (including an experimental cluster-robust option), and how the bridge penalty lambda is selected. The four estimators share the same inferential machinery, so the discussion applies to all of them; we use fetwfe() in the running example.

Two companion vignettes cover the rest of the package’s inference toolkit:

1. The headline: tight Gaussian CIs by default under (Ψ-IF)

By default, the package’s standard errors — the att_se slot on the returned object and the entries of catt_ses — are the tight Gaussian variance from the paper’s Theorem te.asym.norm.thm(c'), valid under the influence-function condition (Ψ-IF) (Faletto 2025, arXiv:2312.05985, Assumption (Ψ-IF) at paper Eq. psi.if.assum). The combined ATT variance has two pieces,

V1=σε2𝛂𝚺1𝛂,V2=vψ(𝛃0,𝒮), V_1 \;=\; \sigma_\varepsilon^2 \, \boldsymbol{\alpha}^\top \boldsymbol{\Sigma}^{-1} \boldsymbol{\alpha}, \qquad V_2 \;=\; v_\psi(\boldsymbol{\beta}_0, \mathcal{S}),

where V1V_1 is the Kock-2013 regression-coefficient variance contribution and V2V_2 is the propensity-score variance contribution; the names V1V_1 and V2V_2 are the paper’s, and the package exposes both of them by name on the fitted object under result$internal$variance_components. Theorem (c') shows that under (Ψ-IF) the same-data confidence interval is asymptotically Gaussian with the tight variance V1+V2V_1 + V_2, identical to the independent-data variance in part (b) of the theorem. Concretely, the package returns

SÊ(τ̂ATT)=Var̂1+Var̂2=ṽN/N, \widehat{\text{SE}}(\hat{\tau}_{\text{ATT}}) \;=\; \sqrt{\widehat{\text{Var}}_1 + \widehat{\text{Var}}_2} \;=\; \sqrt{\tilde{v}_N \,/\, N},

where ṽN:=v̂N/T\tilde{v}_N := \hat{v}_N / T is the unit-scaled paper-notation variance estimator (paper line 2004), and the Wald confidence interval is

CIN(ṽN;α)=[T̂N±z1α/2ṽN/N] \text{CI}_N(\tilde{v}_N;\,\alpha) \;=\; \left[\, \hat{T}_N \;\pm\; z_{1 - \alpha/2}\,\sqrt{\tilde{v}_N \,/\, N} \,\right]

per paper Eq. conf.int.form (line 1259). The same combined variance backs the per-cohort entries of catt_ses and the catt_df confidence-interval columns.

1.1 Which propensity-score estimators satisfy (Ψ-IF)?

Assumption (Ψ-IF) requires the cohort-weight estimator 𝛙̂N\widehat{\boldsymbol{\psi}}_N to admit an i.i.d.-across-ii unit-level influence-function representation. It is satisfied by

  • Cohort sample proportions π̂g=Ng/N\widehat{\pi}_g = N_g / N — the package’s default propensity-score estimator. So out of the box the tight Gaussian CI is asymptotically exact.
  • Multinomial-logit MLE for P(W=g𝐗)P(W = g \mid \boldsymbol{X}).
  • Any GLM fit by maximum likelihood on W𝐗W \mid \boldsymbol{X}.
  • Kernel or series regression of 𝟏{W=g}\mathbf{1}\{W = g\} on 𝐗\boldsymbol{X}.

(See the discussion at paper line 1268: “When (Ψ-IF) holds — which is the case for every standard parametric or semi-parametric propensity-score estimator listed in Remark psi.if.coverage (cohort sample proportions, multinomial logit, any GLM, kernel/series regression) — Theorem te.asym.norm.thm(c') gives asymptotic Gaussianity with the tight variance V1+V2V_1 + V_2.”)

The (Ψ-IF) condition fails for Robins-Rotnitzky-style augmented doubly-robust estimators that augment the propensity score with outcome residuals. The package does not currently implement these; when it does, users will be able to opt out of the tight CI via se_type = "conservative" (see the conservative-fallback note in Section 2.3).

1.2 What changed in v1.12.0

Versions \le 1.11.7 of the package returned the conservative Cauchy-Schwarz upper bound

SÊ(τ̂ATT)=Var̂1+Var̂2+2Var̂1Var̂2 \widehat{\text{SE}}(\hat{\tau}_{\text{ATT}}) \;=\; \sqrt{\widehat{\text{Var}}_1 + \widehat{\text{Var}}_2 + 2\sqrt{\widehat{\text{Var}}_1\,\widehat{\text{Var}}_2}}

as the same-data default (from paper Theorem (c)). Theorem (c') was added in the paper’s same-data Gaussianity upgrade, and v1.12.0 switches the package default to match: the tight Gaussian variance is now the headline number. The conservative formula remains available as an opt-in fallback (Section 2.3) for any users whose propensity-score estimator violates (Ψ-IF).

2. The two variance pieces and how the package exposes them

The combined ATT variance decomposes as V1+V2V_1 + V_2, with both pieces estimated by the package and exposed as named slots.

2.1 V1V_1: regression-coefficient variance

The regression-coefficient piece is

Var̂1=σε2NTψattĜ1ψatt, \widehat{\text{Var}}_1 \;=\; \frac{\sigma_\varepsilon^2}{NT}\;\psi_{\text{att}}^\top \widehat G^{-1} \psi_{\text{att}},

where Ĝ1\widehat G^{-1} is the Gram inverse on the bridge-selected support and ψatt\psi_{\text{att}} encodes the cohort-weighted ATT contrast. The per-cohort SEs (catt_ses) have the same form with a cohort selector ψg\psi_g in place of ψatt\psi_{\text{att}}. This piece is stored on the fitted object as result$internal$variance_components$att_var_1; the paper-notation per-unit-scaled version is result$internal$variance_components$V_1 (with V_1 = N * att_var_1).

2.2 V2V_2: propensity-score variance

The propensity-score piece accounts for variability from estimating the cohort-membership probabilities π̂g\widehat{\pi}_g. It scales like 1/N1/N and is unrelated to the regression residuals. This piece is stored as result$internal$variance_components$att_var_2, with paper-notation per-unit-scaled V_2 = N * att_var_2.

2.3 The combined variance: tilde_v_N and the paper’s catalogue

The unit-scaled total variance is

ṽN=NVar̂(T̂N), \tilde{v}_N \;=\; N \cdot \widehat{\text{Var}}(\hat{T}_N),

stored as result$internal$variance_components$tilde_v_N. The original (un-scaled) paper notation hat_v_N = T * tilde_v_N is also exposed.

The paper catalogues six variance estimators at line 2006; the package exposes the four that apply to the ATT (the other two apply to the CATT(x), which connects to the on-hold predict() work in issue #33):

Slot Paper Eq. Meaning
tilde_v_N_C v.n.r.t.att.const Fixed-π exact (no propensity-score noise)
tilde_v_N_C_pi_hat v.n.r.t.att.rand Random-π exact (tight Gaussian under (Ψ-IF))
tilde_v_N_C_pi_hat_cons v.n.r.t.att.rand.cons Random-π conservative (Cauchy-Schwarz)
tilde_v_N_cons var.est.kock.wooldridge.subgauss.cons Subgaussian conservative bound

The default att_se value is computed from tilde_v_N_C_pi_hat (via sqrt(tilde_v_N_C_pi_hat / N)); se_type = "conservative" switches the headline to tilde_v_N_C_pi_hat_cons. When indep_counts is supplied (two-sample regime), the two-sample-exact formula tilde_v_N_C_pi_hat is in force regardless of se_type.

Conservative fallback (se_type = "conservative"). If a cohort-weight (propensity-score) estimator violates (Ψ-IF) — the canonical case being a Robins-Rotnitzky-style augmented doubly-robust estimator that injects outcome-residual information into the weights — Theorem (c') no longer applies, and se_type = "conservative" switches the headline SE to the Cauchy-Schwarz upper bound Var̂1+Var̂2+2Var̂1Var̂2\sqrt{\widehat{\text{Var}}_1 + \widehat{\text{Var}}_2 + 2\sqrt{\widehat{\text{Var}}_1\,\widehat{\text{Var}}_2}} (tilde_v_N_C_pi_hat_cons, paper Theorem (c)). The package ships no non-(Ψ-IF) estimator — its default cohort-sample-proportions weights π̂g=Ng/N\widehat{\pi}_g = N_g/N satisfy (Ψ-IF) — so under the default estimator there is no practical reason to set it, and default users never need it. It exists for users who supply their own propensity-score estimates (issue #32) or who want to reproduce the pre-v1.12.0 conservative default; see ?fetwfe.

2.4 A running example

The rest of this vignette uses the following small simulated panel and its default fit; the cluster-robust comparison in Section 3.3 reuses them.

set.seed(2026)
sim_coefs <- genCoefs(G = 3, T = 6, d = 2, density = 0.5, eff_size = 2)
# sim_coefs was built without a seed, so draw from the ambient generator
# (seeded just above); seed = NA selects the ambient generator silently.
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
)

3. Experimental: cluster-robust standard errors via se_type = "cluster"

The package’s default standard errors rely on the model-based covariance structure (paper Assumption F1: i.i.d. units, idiosyncratic noise plus a unit random effect). In applied DiD work it is not unusual to suspect one or more violations of F1:

  • Within-unit serial correlation beyond the random effect. The unit random effect cic_i in F1 absorbs a single, time-constant deviation per unit. It does not model serial correlation in the idiosyncratic shocks. Bertrand, Duflo, and Mullainathan (2004) is the classic warning that ignoring within-unit serial correlation can drastically understate DiD standard errors.
  • Heteroskedasticity across units. F1 imposes a single σε2\sigma_\varepsilon^2 across all units.
  • Higher-level clustering. F1 treats units as i.i.d.; if observations within a state-year or industry-year share unobserved shocks, the unit-level F1 variance will understate the true sampling variability.

Starting in version 1.6.0, all four estimators accept an experimental se_type argument with value "cluster":

fetwfe(..., se_type = "cluster")

Setting se_type = "cluster" swaps the model-based regression-coefficient variance V1V_1 for a unit-clustered Liang-Zeger CR1 sandwich computed on the bridge-selected support. The V2V_2 piece and the (Ψ-IF)-tight combination V1+V2V_1 + V_2 are unchanged — only the V1V_1 formula switches.

3.1 The formula

Let Ŝ\widehat{S} be the support selected by the bridge regression, XŜX_{\widehat{S}} the corresponding design matrix in the coordinate system the regression was solved in (GLS-transformed for ETWFE/twfeCovs, fusion-then-GLS-transformed for FETWFE/BETWFE), and ε̂\widehat\varepsilon the residuals from OLS on that selected support. The cluster-robust variance is

VCR=NN1(XŜXŜ)1(i=1NXiŜε̂iε̂iXiŜ)(XŜXŜ)1, V_{\text{CR}} \;=\; \frac{N}{N-1}\; (X_{\widehat{S}}^\top X_{\widehat{S}})^{-1}\; \left(\sum_{i=1}^N X_{i\cdot\widehat{S}}^\top \widehat\varepsilon_{i\cdot} \widehat\varepsilon_{i\cdot}^\top X_{i\cdot\widehat{S}}\right)\; (X_{\widehat{S}}^\top X_{\widehat{S}})^{-1},

with units i=1,,Ni = 1, \dots, N as clusters and an N/(N1)N/(N-1) small-sample adjustment (matching sandwich::vcovCL(cadjust = TRUE, type = "HC0")). The CATT SE for cohort gg is ψgVCRψg\sqrt{\psi_g^\top V_{\text{CR}} \psi_g} (using a zero-padded ψg\psi_g on the full selected support); the ATT regression-coefficient variance is ψattVCRψatt\psi_{\text{att}}^\top V_{\text{CR}} \psi_{\text{att}}, replacing V1V_1 above. The propensity-score piece V2V_2 is unchanged because it depends on empirical cohort proportions, not regression residuals.

For FETWFE and BETWFE, se_type = "cluster" is only meaningful when q < 1 (the bridge oracle property is required); for q >= 1 the cluster path returns NA just like the default. ETWFE and twfeCovs have no q argument, so the cluster path always runs when the Gram matrix is invertible.

3.2 Why we call this experimental

The CR1 sandwich is a textbook estimator and the package’s implementation matches sandwich::vcovCL() to numerical precision on a clean panel without selection. What is not yet covered by the paper’s theory is:

  • Verification that the bridge oracle property of Theorem 6.2 still holds under the relaxed covariance structure that motivates cluster-robust SEs in the first place.
  • Sandwich consistency after model selection, i.e., that the CR1 sandwich evaluated at the bridge-selected support is a consistent estimator of the true asymptotic variance.

These extensions are mechanically routine but conceptually non-trivial, and they are explicitly outside the package’s current scope. Until they land, se_type = "cluster" is exposed as an opt-in, clearly-labelled experimental feature.

Recommendation. Until the theory lands, treat se_type = "cluster" as a sensitivity check: report both se_type = "default" and se_type = "cluster" in applied work, comment on the gap, and lean on the default for headline numbers.

3.3 Worked cluster comparison

res_cluster <- fetwfeWithSimulatedData(sim_data, se_type = "cluster")

c(
  default = res_default$att_se,
  cluster = res_cluster$att_se
)
#>   default   cluster 
#> 0.1050007 0.1000234

On this F1-conforming simulated panel the two SEs are similar by construction: the data-generating process satisfies F1, so the model-based SE is already valid and the cluster-robust SE estimates the same underlying variance. Under a deliberately serially-correlated DGP (or under heteroskedasticity, or higher-level clustering) the cluster-robust SE will typically be larger.

The print() and summary() methods label the SE so it is clear which one was used:

print(res_cluster)
#> Fused Extended Two-Way Fixed Effects Results
#> ===========================================
#> 
#> Overall Average Treatment Effect (ATT):
#>   Estimate:   -0.7798
#>   Std. Error (cluster-robust): 0.1000
#>   P-value:    6.38e-15
#>   Selected:   TRUE
#>   95% CI:    [-0.9758, -0.5838]
#> 
#> Cohort Average Treatment Effects (CATT) [simultaneous 95% CI]:
#>  cohort   estimate         se     ci_low    ci_high     p_value selected
#>       2 -0.4137526 0.06057451 -0.5489811 -0.2785242 1.69269e-11     TRUE
#>       3 -1.6817852 0.13059071 -1.9733200 -1.3902505 0.00000e+00     TRUE
#>       4  0.0000000 0.00000000  0.0000000  0.0000000          NA    FALSE
#> 
#> Event-Study Average Treatment Effects (per event time) [simultaneous 95% CI]:
#>  event_time n_cohorts   estimate        se     ci_low    ci_high      p_value
#>           0         3  0.0000000 0.0000000  0.0000000  0.0000000           NA
#>           1         3 -0.6359733 0.1104477 -0.8975561 -0.3743905 2.162509e-08
#>           2         3 -0.6359733 0.1104477 -0.8975561 -0.3743905 2.529819e-08
#>           3         2 -3.0326244 0.1919425 -3.4872181 -2.5780307 0.000000e+00
#>           4         1  0.1684927 0.2513507 -0.4268026  0.7637881 8.581460e-01
#> 
#> Model Details:
#>   Units (N)           : 120
#>   Time periods (T)    : 6
#>   Treated cohorts (G) : 3
#>   Covariates (d)      : 2
#>   Features (p)        : 62
#>   Selected size       : 26
#>   Lambda*             : 0.0114

The CATT SEs and confidence intervals in catt_df are recomputed from the same cluster-robust sandwich; the CATT p-values follow accordingly.

4. Choosing the bridge penalty parameter lambda

As of v1.13.0, fetwfe() and betwfe() select the bridge penalty lambda via 10-fold cross-validation by default (using grpreg::cv.grpreg). The prior default was BIC over the same grpreg lambda grid. The change addresses a finite-sample bias issue documented in simulation studies (see issue #164 in the package repository): under the prior BIC default, the overall ATT estimator was biased toward zero at moderate sample sizes, producing 95% confidence intervals whose empirical coverage was as low as 0.00 in some regimes. Cross-validation restores near-nominal coverage in every regime tested and additionally has the lowest MSE in the high-dim regime FETWFE was originally designed for.

The lambda_star and lambda_star_model_size result slots are populated identically under both selection paths; only the underlying selection criterion differs. The CV path stores its provenance in three new top-level slots: lambda_selection, cv_folds, and cv_seed. The seed defaults to as.integer(N * T) so consecutive calls on the same dataset are reproducible without the user having to specify a seed.

To recover the prior BIC behavior — for reproducing results from analyses run against v1.12.0 or earlier, for example — pass lambda_selection = "bic":

fit_bic <- fetwfe(..., lambda_selection = "bic")

To control the CV fold assignment explicitly, pass cv_folds and/or cv_seed:

fit_cv <- fetwfe(..., cv_folds = 20L, cv_seed = 42L)

References

Bertrand, M., Duflo, E., & Mullainathan, S. (2004). “How much should we trust differences-in-differences estimates?” Quarterly Journal of Economics 119(1), 249–275.

Faletto, G. (2025). “Fused Extended Two-Way Fixed Effects for Difference-in-Differences with Staggered Adoptions.” arXiv preprint arXiv:2312.05985.

Kock, A. B. (2013). “Oracle Efficient Variable Selection in Random and Fixed Effects Panel Data Models.” Econometric Theory 29(1), 115–152.

Liang, K.-Y., & Zeger, S. L. (1986). “Longitudinal data analysis using generalized linear models.” Biometrika 73(1), 13–22.