Standard Errors

IID, HC1, and clustered

Overview

This tutorial shows the standard error options available in leanfe.

All standard error types use YOCO compression from Wong et al. (2021) combined with sparse matrices for blazing fast computation. The sparse matrix approach is used for all SE types (IID, HC1, and clustered) — not just clustered SEs.

YOCO + Sparse Matrix Strategy

When YOCO compression is triggered (low-cardinality FEs, no IV), leanfe:

  1. Compresses data: GROUP BY (regressors + FEs) → sufficient statistics per group
  2. Builds sparse design matrix: FE dummies stored as scipy.sparse (Python) or Matrix::sparseMatrix (R)
  3. Solves sparse WLS: Weighted least squares using sparse matrix operations
  4. Computes SEs: Method varies by vcov type, but all use the sparse design matrix

This applies to all 4 backend/language combinations (Python/Polars, Python/DuckDB, R/Polars, R/DuckDB).

Noteleanfe’s Contribution: Sparse Matrices

The original YOCO paper by Wong et al. (2021) describes the compression strategy and mathematical framework, but focuses on the statistical methodology without addressing how to efficiently represent FE dummy variables in the design matrix. leanfe extends YOCO with sparse matrices — a natural choice since FE dummies are ~99% zeros. This provides ~100x memory savings and faster matrix operations while producing mathematically identical results. See Benchmark Methodology for the full explanation.

Options

Type Syntax Description Sparse Matrix
IID vcov="iid" Homoskedastic errors ✅ Design matrix
HC1 vcov="HC1" Heteroskedasticity-robust ✅ Design matrix
Clustered vcov="cluster" Correlation within groups ✅ Design matrix + cluster matrix W̃_C

Example Data

import numpy as np
import polars as pl
from leanfe import leanfe

np.random.seed(42)
n_obs = 50_000

firm_id = np.random.randint(1, 101, n_obs)
year = np.random.randint(2015, 2025, n_obs)
treatment = np.random.binomial(1, 0.3, n_obs).astype(float)
x1 = np.random.normal(0, 1, n_obs)

# Heteroskedastic + clustered errors
error = np.random.normal(0, 0.5 + 0.5 * np.abs(x1)) + np.random.normal(0, 0.5, 101)[firm_id]

y = 2.5 * treatment + 1.5 * x1 + np.random.normal(0, 1, 101)[firm_id] + error

df = pl.DataFrame({
    "y": y, "treatment": treatment, "x1": x1, "firm_id": firm_id, "year": year
})

print(f"Data: {n_obs:,} obs, {df['firm_id'].n_unique()} firms")
Data: 50,000 obs, 100 firms

IID

result_iid = leanfe(data=df, formula="y ~ treatment + x1 | firm_id + year", vcov="iid")

print("IID:")
for var in ["treatment", "x1"]:
    print(f"  {var}: {result_iid['coefficients'][var]:.4f} (SE: {result_iid['std_errors'][var]:.4f})")
IID:
  treatment: 2.4790 (SE: 0.0094)
  x1: 1.5061 (SE: 0.0043)

HC1 (Robust)

result_hc1 = leanfe(data=df, formula="y ~ treatment + x1 | firm_id + year", vcov="HC1")

print("HC1:")
for var in ["treatment", "x1"]:
    print(f"  {var}: {result_hc1['coefficients'][var]:.4f} (SE: {result_hc1['std_errors'][var]:.4f})")
HC1:
  treatment: 2.4790 (SE: 0.0095)
  x1: 1.5061 (SE: 0.0061)

Clustered

Clustered SEs use the same YOCO + sparse design matrix as IID/HC1, plus an additional **sparse cluster matrix W̃_C** (Section 5.3.1 of Wong et al. 2021) for efficient score aggregation:

  1. Within-cluster compression: cluster ID added to GROUP BY (each compressed record → one cluster)
  2. Sparse design matrix: Same as IID/HC1 — FE dummies as sparse matrix
  3. **Sparse cluster matrix W̃_C**: Indicator matrix (G × C) mapping groups to clusters
  4. Vectorized aggregation: cluster_scores = W̃_Cᵀ @ (X × ẽ⁰) — no loops over clusters!
result_cluster = leanfe(
    data=df,
    formula="y ~ treatment + x1 | firm_id + year",
    vcov="cluster",
    cluster_cols=["firm_id"]
)

print("Clustered (firm):")
for var in ["treatment", "x1"]:
    print(f"  {var}: {result_cluster['coefficients'][var]:.4f} (SE: {result_cluster['std_errors'][var]:.4f})")
print(f"Clusters: {result_cluster['n_clusters']}")
print(f"Strategy: {result_cluster.get('strategy', 'N/A')}")
Clustered (firm):
  treatment: 2.4790 (SE: 0.0229)
  x1: 1.5061 (SE: 0.0061)
Clusters: 100
Strategy: compress

Two-Way Clustering

result_twoway = leanfe(
    data=df,
    formula="y ~ treatment + x1 | firm_id + year",
    vcov="cluster",
    cluster_cols=["firm_id", "year"]
)

print("Two-way (firm + year):")
for var in ["treatment", "x1"]:
    print(f"  {var}: {result_twoway['coefficients'][var]:.4f} (SE: {result_twoway['std_errors'][var]:.4f})")
Two-way (firm + year):
  treatment: 2.4790 (SE: 0.0229)
  x1: 1.5061 (SE: 0.0061)

Comparison

print("SE comparison for 'treatment':")
print(f"  IID:       {result_iid['std_errors']['treatment']:.4f}")
print(f"  HC1:       {result_hc1['std_errors']['treatment']:.4f}")
print(f"  Clustered: {result_cluster['std_errors']['treatment']:.4f}")
print(f"  Two-way:   {result_twoway['std_errors']['treatment']:.4f}")
SE comparison for 'treatment':
  IID:       0.0094
  HC1:       0.0095
  Clustered: 0.0229
  Two-way:   0.0229

Technical Notes

When is YOCO + Sparse Matrix Used?

The YOCO compression path with sparse matrices is triggered when:

  • vcov is “iid”, “HC1”, or “cluster” (all three!)
  • No IV/2SLS instruments
  • FE cardinality is below thresholds (max single FE ≤ 10,000, total FE levels ≤ 20,000)

When these conditions are met, all 4 implementations (Python/Polars, Python/DuckDB, R/Polars, R/DuckDB) use the same sparse matrix strategy.

Fallback to FWL Demeaning

If FE cardinality is too high or IV is used, leanfe falls back to iterative FWL demeaning (no sparse matrix). This is faster for high-cardinality FEs where the sparse matrix would be huge.

Implementation Details

See Benchmark Methodology for the full algorithm description.

Next Steps