Measuring Uncertainty

The first step in selecting standard errors is to identify the sources of uncertainty. Uncertainty generally takes two forms: sampling uncertainty and design uncertainty.

Sampling uncertainty arises when we observe a sample and wish to make claims about the population. As noted in the previous section, for a sample mean \(\bar{x}\) drawn from a population with variance \(\sigma^2\), the standard error is defined as: \[SE(\bar{x}) = \frac{\sigma}{\sqrt{n}}\]

Here \(\sigma\) is the population standard deviation; \(1/\sqrt{n}\) reflects that larger samples yield more precise estimates, a consequence of \(\text{Var}(\bar{x}) = \sigma^2/n\).

The animation below shows this shrinkage directly: as the sample size increases, the distribution of sample means tightens around the true population mean, with the spread narrowing at the \(1/\sqrt{n}\) rate.

Design uncertainty originates from the choices researchers make about how to collect and structure data. When the unit of sampling or treatment differs from the unit of analysis (for example, households sampled, individuals analyzed; states treated, counties measured), standard errors must account for the within-group correlation this design creates. The right SE is determined not by what the data look like, but by how they were generated (Abadie et al. 2023; see also Abadie et al. 2020 for the sampling-based versus design-based framework underlying this distinction).

Working Example: Nunn and Wantchekon (2011)

Nunn and Wantchekon (2011) ask whether the Atlantic and Indian Ocean slave trades, which forcibly removed tens of millions of people from Africa between the fifteenth and twentieth centuries, left a lasting legacy of distrust. Using Afrobarometer survey data merged with Nunn (2008)’s historical estimates of slave exports by ancestral ethnic group, they find that individuals today whose ancestors experienced more intensive slave raiding report lower trust in neighbors, relatives, and local institutions.

The dataset is a natural fit for a guide on standard errors. The key regressor, historical slave exports, is measured at the ethnic-group level, and both within-country error correlation and substantial cross-country variation in colonial history are key features of the data. Each standard error correction covered in this guide addresses a different aspect of that structure.

For instructional purposes, this guide uses a simplified bivariate model, trust_neighbors ~ ln_exports, to demonstrate standard-error mechanics with real data. This working example omits controls, fixed effects, survey-design details, and other modeling choices that would matter for a faithful replication of Nunn and Wantchekon (2011). Some sections use this model to show implementation mechanics rather than to recommend a final inferential strategy. When a method requires a different data structure, we say so and use a more appropriate example.

Robust (HC) Standard Errors

nunn       <- haven::read_dta("data/Nunn_Wantchekon_AER_2011.dta")
regression <- lm(trust_neighbors ~ ln_exports, data = nunn)

par(mfrow = c(2, 2))
plot(regression, pch = 19, cex = 0.8)

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.nonparametric.api as nparam
import matplotlib.pyplot as plt
import scipy.stats as stats

nunn = pd.read_stata("data/Nunn_Wantchekon_AER_2011.dta")
reg1 = sm.formula.ols("trust_neighbors ~ ln_exports", data=nunn).fit()
residuals = reg1.resid
fitted_values = reg1.fittedvalues
# Internally studentized residuals for the scale-location plot
std_residuals = reg1.get_influence().resid_studentized_internal

fig, axes = plt.subplots(2, 2, figsize=(12, 9))

axes[0, 0].scatter(fitted_values, residuals, alpha=0.6)
axes[0, 0].axhline(0, color="red", linestyle="--")
lowess = nparam.lowess(residuals, fitted_values, frac=0.3)
axes[0, 0].plot(lowess[:, 0], lowess[:, 1], color="blue", linewidth=2)
axes[0, 0].set(title="Residuals vs Fitted", xlabel="Fitted Values", ylabel="Residuals")

stats.probplot(residuals, dist="norm", plot=axes[0, 1])

((array([-3.98542269, -3.76960615, -3.65152731, ...,  3.65152731,
        3.76960615,  3.98542269], shape=(20580,)), array([-1.74997556, -1.74997556, -1.74997556, ...,  1.27013099,
        1.27013099,  1.27013099], shape=(20580,))), (np.float64(0.9424317803632147), np.float64(5.932626781982231e-16), np.float64(0.9331372837619768)))

axes[0, 1].set_title("Normal Q-Q")

axes[1, 0].hist(residuals, bins=20, alpha=0.7, color="steelblue", edgecolor="white")
axes[1, 0].set(title="Distribution of Residuals", xlabel="Residuals", ylabel="Frequency")

sqrt_std = np.sqrt(np.abs(std_residuals))
axes[1, 1].scatter(fitted_values, sqrt_std, alpha=0.6)
lowess2 = nparam.lowess(sqrt_std, fitted_values, frac=0.3)
axes[1, 1].plot(lowess2[:, 0], lowess2[:, 1], color="red", linewidth=2)
axes[1, 1].set(title="Scale-Location", xlabel="Fitted Values", ylabel="√|Std. Residuals|")

plt.tight_layout()
plt.savefig("images/diagnostic_plots.png", dpi=120)
plt.show()

use "data/Nunn_Wantchekon_AER_2011.dta", clear
regress trust_neighbors ln_exports

predict fitted_vals, xb
predict residuals, resid
predict std_residuals, rstandard
gen sqrt_abs_std = sqrt(abs(std_residuals))
predict leverage, leverage

scatter residuals fitted_vals, yline(0) title("Residuals vs. Fitted") name(p1, replace)
qnorm residuals, title("Normal Q-Q") name(p2, replace)
scatter sqrt_abs_std fitted_vals, title("Scale-Location") name(p3, replace)
scatter std_residuals leverage, title("Residuals vs. Leverage") name(p4, replace)

graph combine p1 p2 p3 p4, rows(2)

      Source |       SS           df       MS      Number of obs   =    20,580
-------------+----------------------------------   F(1, 20578)     =      0.97
       Model |  .989159647         1  .989159647   Prob > F        =    0.3247
    Residual |  20984.2025    20,578  1.01973965   R-squared       =    0.0000
-------------+----------------------------------   Adj R-squared   =   -0.0000
       Total |  20985.1916    20,579  1.01973816   Root MSE        =    1.0098

------------------------------------------------------------------------------
trust_neig~s | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
  ln_exports |   .0029778   .0030234     0.98   0.325    -.0029484     .008904
       _cons |   1.729869   .0093682   184.65   0.000     1.711507    1.748231
------------------------------------------------------------------------------

(120 missing values generated)

(1,242 missing values generated)

(1,242 missing values generated)

(1,242 missing values generated)

(120 missing values generated)

using GLM, DataFrames, ReadStataTables, Plots, Distributions

nunn     = DataFrame(readstat("data/Nunn_Wantchekon_AER_2011.dta"))
model    = lm(@formula(trust_neighbors ~ ln_exports), nunn)
resids   = residuals(model)
fitted   = predict(model)
std_res  = resids ./ std(resids)

p1 = scatter(fitted, resids, xlabel="Fitted Values", ylabel="Residuals",
             title="Residuals vs Fitted", legend=false, markeralpha=0.6)
hline!(p1, [0], linestyle=:dash, color=:red)

# NOTE: this panel plots raw residuals against a normal scaled to their SD,
# whereas the R/Python tabs plot standardized residuals; the two are equivalent
# up to the estimated scale and may look slightly different.
p2 = qqplot(Normal(0, std(resids)), resids, title="Normal Q-Q",
            xlabel="Theoretical Quantiles", ylabel="Sample Quantiles")

p3 = scatter(fitted, sqrt.(abs.(std_res)), xlabel="Fitted Values",
             ylabel="√|Std. Residuals|", title="Scale-Location",
             legend=false, markeralpha=0.6)

plot(p1, p2, p3, layout=(1, 3), size=(1100, 380))

I am running an OLS regression in [R / Python / Stata / Julia].

Model:       [dependent variable] ~ [predictor 1, predictor 2, ...]
Sample size: n = [N]
Data:        [cross-section / panel / time series]

From my residual diagnostic plots I observe:
  - Residuals vs Fitted: [fan shape / funnel / random scatter / other]
  - Scale-Location:      [flat LOWESS / upward slope / other]
  - Q-Q Plot:            [on the line / heavy tails / curvature]

Please:
1. Confirm whether this is consistent with heteroskedasticity.
2. Recommend which HC estimator (HC0-HC5) is appropriate given
   my sample size and research design.
3. Provide the code to apply it and compare the original and
   corrected standard errors side by side.

nunn       <- haven::read_dta("data/Nunn_Wantchekon_AER_2011.dta")
regression <- lm(trust_neighbors ~ ln_exports, data = nunn)

regression_robust <- lmtest::coeftest(
  regression, vcov = sandwich::vcovHC(regression, type = "HC2")
)
regression_robust

t test of coefficients:

             Estimate Std. Error  t value Pr(>|t|)    
(Intercept) 1.7298690  0.0093408 185.1945   <2e-16 ***
ln_exports  0.0029778  0.0031236   0.9533   0.3404    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

import statsmodels.api as sm
import pandas as pd

nunn = pd.read_stata("data/Nunn_Wantchekon_AER_2011.dta")
reg = sm.OLS.from_formula("trust_neighbors ~ ln_exports", data=nunn).fit(cov_type="HC2")
print(reg.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:        trust_neighbors   R-squared:                       0.000
Model:                            OLS   Adj. R-squared:                 -0.000
Method:                 Least Squares   F-statistic:                    0.9088
Date:                Tue, 18 Aug 2026   Prob (F-statistic):              0.340
Time:                        14:55:50   Log-Likelihood:                -29402.
No. Observations:               20580   AIC:                         5.881e+04
Df Residuals:                   20578   BIC:                         5.882e+04
Df Model:                           1                                         
Covariance Type:                  HC2                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      1.7299      0.009    185.195      0.000       1.712       1.748
ln_exports     0.0030      0.003      0.953      0.340      -0.003       0.009
==============================================================================
Omnibus:                     7464.971   Durbin-Watson:                   1.464
Prob(Omnibus):                  0.000   Jarque-Bera (JB):             1166.555
Skew:                          -0.225   Prob(JB):                    4.85e-254
Kurtosis:                       1.924   Cond. No.                         4.32
==============================================================================

Notes:
[1] Standard Errors are heteroscedasticity robust (HC2)

use "data/Nunn_Wantchekon_AER_2011.dta", clear
regress trust_neighbors ln_exports, vce(hc2)

Linear regression                               Number of obs     =     20,580
                                                F(1, 20578)       =       0.91
                                                Prob > F          =     0.3404
                                                R-squared         =     0.0000
                                                Root MSE          =     1.0098

------------------------------------------------------------------------------
             |             Robust HC2
trust_neig~s | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
  ln_exports |   .0029778   .0031236     0.95   0.340    -.0031447    .0091003
       _cons |   1.729869   .0093408   185.19   0.000      1.71156    1.748178
------------------------------------------------------------------------------

using GLM, DataFrames, ReadStataTables, CovarianceMatrices

nunn  = DataFrame(readstat("data/Nunn_Wantchekon_AER_2011.dta"))
model = lm(@formula(trust_neighbors ~ ln_exports), nunn)

stderror(model, HC2())

I have confirmed heteroskedasticity in my OLS model.

Model:       [dependent variable] ~ [predictors]
Sample size: n = [N]
Design:      [observational / randomized experiment / survey]

1. I am choosing between HC2 and HC3. Which is more appropriate
   for my context and why?
2. Provide the code in [R / Python / Stata / Julia] to apply it.
3. My original SE on [key coefficient] was [X] and after applying
   HC[n] it is [Y]. Is this change meaningful, and what does it
   suggest about the error structure?

Cluster-Robust Standard Errors

HC standard errors handle heteroskedasticity in errors that are independent across observations. But independence itself breaks down when observations share something in common by design: students sampled from the same classroom, survey respondents drawn from the same village, firms assigned treatment as a block within the same industry, the same individual measured repeatedly over time, voters exposed to the same precinct-level shock. The trigger for clustering is the design that generated the data: sampling, treatment assignment, repeated measures, or shared exposure. It is not simply an observation that errors within a group happen to be correlated. HC standard errors are still wrong under any of these designs. We need cluster-robust standard errors.

The cluster-robust (“Liang-Zeger”) variance estimator keeps the sandwich structure but builds the meat from cluster-level sums rather than observation-level ones:

\[\hat{V}_{CR}(\hat{\boldsymbol{\beta}}) = (\mathbf{X}^\top\mathbf{X})^{-1}\left(\sum_{g=1}^{G}\mathbf{X}_g^\top \hat{\mathbf{e}}_g \hat{\mathbf{e}}_g^\top \mathbf{X}_g\right)(\mathbf{X}^\top\mathbf{X})^{-1}\]

where \(g = 1, \ldots, G\) indexes clusters, \(\mathbf{X}_g\) stacks the rows of \(\mathbf{X}\) for cluster \(g\), and \(\hat{\mathbf{e}}_g\) is that cluster’s residual vector. The meat allows arbitrary within-cluster correlation while assuming independence across clusters (Liang and Zeger 1986). Software typically multiplies by a finite-sample factor: Stata’s cluster default is \(c = \tfrac{G}{G-1}\cdot\tfrac{N-1}{N-k}\), while sandwich::vcovCL applies \(\tfrac{G}{G-1}\) by default (the full Stata factor is available through its arguments). Inference then uses \(t_{G-1}\) rather than the normal. Consistency requires \(G \to \infty\), which is why few-cluster designs need the corrections discussed below.

The key design questions before you cluster:

• What creates the clustering? The sampling design, treatment assignment mechanism, repeated-measures structure, or shared exposure determines the right clustering level, not the presence of within-group correlation in the outcome. If you sampled districts and surveyed individuals within them, cluster at the district. If treatment was assigned at the school level, cluster at the school, regardless of whether you think school-level unobservables exist. If the same unit is measured repeatedly over time, cluster at the unit. If units share a common exposure (a market, a policy, a shock), cluster at the level of that exposure (Abadie et al. 2023).
• Cluster at the right level. Cluster at the level where sampling or assignment introduced dependence. Clustering too low (at the individual when treatment is at the group) understates uncertainty and inflates test statistics. Clustering too high wastes precision.
• Count your clusters. Cluster-robust inference degrades when the number of clusters is small, roughly below 30 to 50 (Cameron and Miller 2015; Jackson 2020). This 30-to-50 figure is a rough heuristic, not a threshold: what actually governs reliability is the effective degrees of freedom for your contrast (see the near-census and skewed-design FAQs), which can be far below \(G - 1\) when treatment is concentrated in a few clusters. With few clusters, consider the wild cluster bootstrap, CR2 corrections via clubSandwich, or the finite-sample adjustments available through fixest::ssc() where appropriate.
• Weights change the variance estimator. If you run WLS or use survey/probability weights, the ordinary cluster sandwich (sandwich::vcovCL, the feols default) is not consistent for the weighted estimator. Use clubSandwich::coef_test(fit, vcov = "CR2", cluster = ~ g) (or a survey-design estimator) instead; the leverage correction is what restores validity under non-constant weights.
• Two-way clustering. When two sources of correlation exist (say, by firm and by year in a panel), two-way clustering is an option. It requires enough clusters in both dimensions to be reliable; with sparse panels, the resulting variance matrix can be non-positive definite. When that happens, fixest (and multiwayvcov) can apply an eigenvalue-flooring fix (setting negative eigenvalues to zero), but a non-positive-definite result usually signals that one dimension has too few clusters: fall back to one-way clustering on the dimension with more clusters, or use the wild cluster bootstrap with \(\min(G_A, G_B)\) as the effective cluster count.

In Nunn and Wantchekon (2011) the key regressor, log slave exports, is defined at the level of a respondent’s ancestral ethnic group: everyone whose ancestry traces to the same group is assigned the same value. That is the design signature that calls for clustering: the right-hand-side variation lives at the group level, so the errors are almost certainly correlated within group as well. We cluster at the ethnic group (murdock_name), of which there are about 187:

nunn       <- haven::read_dta("data/Nunn_Wantchekon_AER_2011.dta")
regression <- lm(trust_neighbors ~ ln_exports, data = nunn)

regression_clustered <- lmtest::coeftest(
  regression, vcov = sandwich::vcovCL(regression, cluster = ~murdock_name)
)
regression_clustered

t test of coefficients:

             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.7298690  0.0526941 32.8285   <2e-16 ***
ln_exports  0.0029778  0.0285699  0.1042    0.917    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

import statsmodels.api as sm
import pandas as pd

nunn = pd.read_stata("data/Nunn_Wantchekon_AER_2011.dta")
nunn = nunn.dropna(subset=["trust_neighbors", "ln_exports", "murdock_name"])
reg = sm.OLS.from_formula("trust_neighbors ~ ln_exports", data=nunn).fit(
    cov_type="cluster",
    cov_kwds={"groups": nunn["murdock_name"]}
)
# NOTE: statsmodels reports normal (z) critical values and a different
# small-sample adjustment than R/Stata. With many clusters (here ~187
# ethnic groups) the difference is negligible, but with few clusters it
# understates uncertainty: apply t(G-1) critical values manually and prefer
# the wild cluster bootstrap (see the Bootstrap section). pyfixest offers
# CRV3 and the wild bootstrap natively.
print(reg.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:        trust_neighbors   R-squared:                       0.000
Model:                            OLS   Adj. R-squared:                 -0.000
Method:                 Least Squares   F-statistic:                   0.01086
Date:                Tue, 18 Aug 2026   Prob (F-statistic):              0.917
Time:                        14:55:51   Log-Likelihood:                -29402.
No. Observations:               20580   AIC:                         5.881e+04
Df Residuals:                   20578   BIC:                         5.882e+04
Df Model:                           1                                         
Covariance Type:              cluster                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      1.7299      0.053     32.828      0.000       1.627       1.833
ln_exports     0.0030      0.029      0.104      0.917      -0.053       0.059
==============================================================================
Omnibus:                     7464.971   Durbin-Watson:                   1.464
Prob(Omnibus):                  0.000   Jarque-Bera (JB):             1166.555
Skew:                          -0.225   Prob(JB):                    4.85e-254
Kurtosis:                       1.924   Cond. No.                         4.32
==============================================================================

Notes:
[1] Standard Errors are robust to cluster correlation (cluster)

use "data/Nunn_Wantchekon_AER_2011.dta", clear
regress trust_neighbors ln_exports, cluster(murdock_name)

Linear regression                               Number of obs     =     20,580
                                                F(1, 184)         =       0.01
                                                Prob > F          =     0.9171
                                                R-squared         =     0.0000
                                                Root MSE          =     1.0098

                         (Std. err. adjusted for 185 clusters in murdock_name)
------------------------------------------------------------------------------
             |               Robust
trust_neig~s | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
  ln_exports |   .0029778   .0285699     0.10   0.917     -.053389    .0593445
       _cons |   1.729869   .0526941    32.83   0.000     1.625907    1.833831
------------------------------------------------------------------------------

using GLM, DataFrames, ReadStataTables, CovarianceMatrices

nunn   = DataFrame(readstat("data/Nunn_Wantchekon_AER_2011.dta"))
model  = lm(@formula(trust_neighbors ~ ln_exports), nunn)
groups = nunn.murdock_name

stderror(model, CRHC2(groups))

I am deciding whether and how to cluster my standard errors.

Data structure:       [cross-section / panel / survey]
Sampling design:      [how were units selected?]
Treatment assignment: [individual level / group level (what group?)]
Cluster variable:     [candidate variable]
Number of clusters:   G = [?]

1. Based on my design, should I cluster? At what level?
2. With G = [N] clusters, are cluster-robust SEs reliable, or
   should I use the wild cluster bootstrap instead?
3. Provide the code in [R / Python / Stata / Julia] to implement
   the appropriate correction.

Spatial Standard Errors

When your units of observation have a geographic location (countries, regions, counties, grid cells), the independence assumption can fail not because of sampling or treatment design, but because nearby units resemble each other. Soil quality, climate, economic conditions, and institutional history all exhibit spatial continuity. Residuals that share this structure are spatially autocorrelated, and standard errors that ignore it are biased downward.

The standard fix is Conley (1999) standard errors, which extend the sandwich estimator by allowing correlation between observations as a function of geographic distance. Observations within a specified distance cutoff are treated as potentially correlated; those farther apart are treated as independent:

\[\hat{V}_{spatial}(\hat{\boldsymbol{\beta}}) = (\mathbf{X}^\top\mathbf{X})^{-1} \left(\sum_i \sum_j \kappa\!\left(\frac{d_{ij}}{h}\right) \mathbf{x}_i \mathbf{x}_j^\top \hat{e}_i \hat{e}_j \right) (\mathbf{X}^\top\mathbf{X})^{-1}\]

where \(d_{ij} = d_{ji}\) is the distance between units \(i\) and \(j\), \(h\) is the bandwidth (distance cutoff), and \(\kappa(\cdot)\) is a kernel with \(\kappa(0) = 1\) and \(\kappa(u) = 0\) for \(u > 1\), so that pairs farther apart than \(h\) receive zero weight. The uniform kernel \(\kappa(u) = \mathbf{1}\{u \le 1\}\) gives Conley’s original hard-cutoff form (arbitrary correlation within \(h\), independence beyond); the Bartlett kernel \(\kappa(u) = (1 - u)_+\) instead lets correlation decay linearly to zero at \(h\) and guarantees a positive-semidefinite \(\hat{V}\). The \(i = j\) terms reproduce the heteroskedasticity-robust (HC0) diagonal.

Conley SEs require coordinates, such as latitude and longitude, or some other distance metric, for every observation.¹ Without them you cannot compute \(d_{ij}\) and the estimator is not defined. The bandwidth \(h\) is the most consequential choice, and there is no universal rule for setting it; it should reflect the plausible range of spatial correlation in your subject matter. Standard practice is to try several values and report sensitivity across them. If your conclusions change with the bandwidth, that is worth flagging rather than resolving by choosing a convenient default. Conley and Kelly (2025) documents that bandwidth choice in historical persistence regressions is frequently consequential and sometimes selected in ways that favor significance; that caution applies broadly. For the placebo-based spatial inference of Conley and Kelly (2025) specifically, spatial-basis de-trending plus a synthetic-treatment placebo test, the authors maintain the spatInfer R package (github.com/morganwkelly/spatInfer), which automates basis selection and the placebo reference distribution and is the most direct way to reproduce their recommended workflow.

Before reaching for Conley SEs, also consider whether simple regional clustering would suffice. If spatial correlation is well captured by a discrete grouping (continent, province, watershed), clustering by that group is simpler and often sufficient. Conley SEs earn their complexity when correlation decays continuously with distance rather than cleanly within administrative borders. For geographic coordinates, use great-circle distance rather than Euclidean distance; most implementations offer a spherical distance option.

library(fixest)

nunn <- haven::read_dta("data/Nunn_Wantchekon_AER_2011.dta")

# Conley SEs: conley(cutoff, distance) ~ lat_col + lon_col
# cutoff 500 is in KILOMETERS because distance = "spherical"; with the default
# Euclidean distance the cutoff is in the coordinates' own units (degrees here).
model_conley <- feols(
  trust_neighbors ~ ln_exports,
  data = nunn,
  vcov = conley(500, distance = "spherical") ~ centroid_lat + centroid_long
)
summary(model_conley)

OLS estimation, Dep. Var.: trust_neighbors
Observations: 20,580
Standard-errors: Conley (500km) 
            Estimate Std. Error   t value  Pr(>|t|)    
(Intercept) 1.729869   0.091664 18.871797 < 2.2e-16 ***
ln_exports  0.002978   0.040882  0.072837   0.94194    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 1.00977   Adj. R2: -1.457e-6

# NOTE: There is no mature drop-in Conley (1999) estimator in Python. The code
# below computes the Kelejian-Prucha spatial HAC (spreg, robust="hac") -- related
# but NOT identical to Conley. For genuine Conley SEs, use fixest::feols(vcov =
# conley()) in R (see the R tab), or compute the distance-kernel sandwich by hand.
import pandas as pd
import numpy as np
from spreg import OLS

ModuleNotFoundError: No module named 'spreg'

from libpysal.weights import Kernel

ModuleNotFoundError: No module named 'libpysal'

nunn = pd.read_stata("data/Nunn_Wantchekon_AER_2011.dta")

# Drop rows missing coordinates or model variables BEFORE building the
# kernel weights, so the weights object and y/X line up row for row
nunn = nunn.dropna(subset=["trust_neighbors", "ln_exports",
                            "centroid_long", "centroid_lat"])
coords = list(zip(nunn["centroid_long"], nunn["centroid_lat"]))

# Note: bandwidth is in decimal degrees. One degree ≈ 111 km at the equator
# but shrinks toward the poles. For data spanning all of Africa, a bandwidth
# of 4.5 corresponds roughly to 350–500 km depending on latitude; treat as
# approximate. Use projected coordinates (meters) for a precise 500 km cutoff.
gwk = Kernel(coords, fixed=True, bandwidth=4.5, function="triangular",
             diagonal=True)

NameError: name 'Kernel' is not defined

y = nunn[["trust_neighbors"]].values
X = nunn[["ln_exports"]].values

# spreg.OLS with robust="hac" implements the Kelejian-Prucha (1999) spatial HAC
# estimator, similar in spirit to Conley (1999) but different in mechanics.
# It requires a KERNEL weights object passed via gwk= (not a binary
# contiguity/distance-band matrix passed via w=).
# For the Conley (1999) distance-kernel sandwich estimator, use fixest::feols()
# with vcov = conley() in R (see R tab).
model = OLS(y, X, gwk=gwk, robust="hac", name_y="trust_neighbors",
            name_x=["ln_exports"])

NameError: name 'OLS' is not defined

print(model.summary)

NameError: name 'model' is not defined

* Requires: ssc install acreg
use "data/Nunn_Wantchekon_AER_2011.dta", clear

acreg trust_neighbors ln_exports, ///
  lat(centroid_lat) lon(centroid_long) dist(500) spherical

# Spatial SEs are less mature in Julia than in R or Stata.
# The closest current option is CovarianceMatrices.jl with a
# user-supplied distance weight matrix.

using GLM, DataFrames, ReadStataTables, LinearAlgebra

nunn  = DataFrame(readstat("data/Nunn_Wantchekon_AER_2011.dta"))
model = lm(@formula(trust_neighbors ~ ln_exports), nunn)

# Build a distance-based weight matrix (great-circle, 500 km cutoff)
# using your preferred geodistance function, then pass to a
# custom sandwich estimator. Dedicated spatial SE packages for
# Julia are in active development.

I have cross-sectional data with geographic coordinates and
suspect spatial autocorrelation in my residuals.

Model:            [dependent variable] ~ [predictors]
Unit of obs:      [countries / counties / grid cells / other]
Coordinates:      [latitude / longitude variable names]
Plausible range:  [how far, in km, do you expect correlation to reach?]
Software:         [R / Python / Stata / Julia]

1. Is Conley SEs the right approach, or would regional clustering
   suffice given my unit of observation?
2. How sensitive should I expect my results to be to the choice
   of bandwidth? Suggest a range of values to try.
3. Provide the code to implement Conley SEs with bandwidth h = [X] km
   and show how to check sensitivity across h = [X, Y, Z].

HAC Standard Errors

In time series and panel data, observations close in time are rarely independent. Today’s unemployment rate is shaped by yesterday’s; a firm’s quarterly earnings exhibit persistence; policy shocks carry over across periods. This serial correlation in errors biases OLS standard errors downward, making results appear more precise than they are.

HAC (heteroskedasticity and autocorrelation consistent) standard errors extend the sandwich estimator by weighting the cross-products of residuals with a kernel that decays as the lag distance grows:

\[\hat{V}_{HAC}(\hat{\boldsymbol{\beta}}) = (\mathbf{X}^\top\mathbf{X})^{-1} \hat{S} (\mathbf{X}^\top\mathbf{X})^{-1}\]

where the long-run variance is estimated as:

\[\hat{S} = \hat{\Gamma}_0 + \sum_{j=1}^{L} \kappa\!\left(\frac{j}{L+1}\right)\left(\hat{\Gamma}_j + \hat{\Gamma}_j^\top\right), \qquad \hat{\Gamma}_j = \frac{1}{T}\sum_{t=j+1}^{T} \mathbf{x}_t \mathbf{x}_{t-j}^\top \hat{e}_t \hat{e}_{t-j}\]

\(\kappa(\cdot)\) is a kernel function and \(L\) is the bandwidth (number of lags). The Newey and West (1987) estimator uses the Bartlett kernel \(\kappa(z) = 1 - |z|\) for \(|z| \leq 1\) (and \(0\) otherwise), giving lag-\(j\) autocovariance the weight \(1 - j/(L+1)\), which declines linearly from 1 at \(j=0\) to \(1/(L+1)\) at the last included lag \(j=L\) and is exactly 0 for \(j > L\). This tapering guarantees that \(\hat{S}\) is positive semi-definite. The \(\hat{\Gamma}_0\) term is the heteroskedasticity-robust (HC) component, so HAC nests HC.

The central applied choice is the lag length \(L\), which controls how many periods of autocorrelation the estimator accounts for. A common rule of thumb is \(L = \lfloor 4(T/100)^{2/9} \rfloor\) (Newey and West 1994), which gives roughly four lags at \(T = 100\); note that the MSE-optimal Bartlett bandwidth grows at rate \(T^{1/3}\), so for long series a faster-growing \(L\) may be warranted, and the two rules can give materially different \(L\). R’s NeweyWest() supports automatic, data-driven bandwidth selection via Andrews (1991) as an alternative to the rule of thumb, whereas Stata’s base newey command requires an explicit lag() and has no built-in data-driven bandwidth selection.² When results are sensitive to \(L\), report them across a range of values.

In short panels, clustering by time period is often preferable to HAC, as it makes fewer assumptions about the temporal correlation structure and sidesteps the bandwidth choice. HAC is the natural choice for long time series. When panel data exhibit both within-unit serial correlation and cross-sectional dependence, two-way clustering (by unit and by time) handles both simultaneously; Driscoll and Kraay (1998) standard errors are a more principled alternative that extends Newey-West to allow spatio-temporal covariance decay; plm::vcovSCC() in R implements them. Driscoll-Kraay relies on a long time dimension, however; with short \(T\) it is unreliable, and two-way clustering or time fixed effects are safer. Before applying any correction, a Breusch-Godfrey test (Breusch 1978; Godfrey 1978) can confirm that serial correlation is actually present; a Durbin-Watson statistic near 2 suggests no problem, while values well below 2 indicate positive serial correlation.

The Nunn and Wantchekon data is a cross-section: it has no time dimension for Newey-West to act on, so applying HAC to it would be meaningless. To demonstrate HAC where it actually applies, we instead simulate a short time series with serially correlated errors and run the corrections on it.

set.seed(123)
T_obs <- 200      # time periods
phi   <- 0.7      # AR(1) coefficient of the errors (serial correlation)
beta  <- 0.5      # true slope
x <- rnorm(T_obs)
e <- as.numeric(arima.sim(model = list(ar = phi), n = T_obs))  # serially correlated errors
y <- 1 + beta * x + e
hac_demo <- data.frame(t = seq_len(T_obs), y = y, x = x)
write.csv(hac_demo, "data/hac_demo_timeseries.csv", row.names = FALSE)

library(sandwich)
library(lmtest)

ts_data <- read.csv("data/hac_demo_timeseries.csv")
model   <- lm(y ~ x, data = ts_data)

# OLS standard errors (baseline)
coeftest(model)

t test of coefficients:

            Estimate Std. Error t value  Pr(>|t|)    
(Intercept) 1.014312   0.091465 11.0896 < 2.2e-16 ***
x           0.471439   0.097217  4.8494 2.501e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Newey-West with automatic lag selection
# NeweyWest() prewhitens by default (prewhite = TRUE); the explicit-lag call
# below turns it off, so the two lines differ in prewhitening as well as in lag
# selection. Set prewhite = FALSE here too to isolate the effect of the lag.
coeftest(model, vcov = NeweyWest(model))

t test of coefficients:

            Estimate Std. Error t value  Pr(>|t|)    
(Intercept) 1.014312   0.197251  5.1422 6.500e-07 ***
x           0.471439   0.091833  5.1336 6.768e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Newey-West with explicit lag
coeftest(model, vcov = NeweyWest(model, lag = 4, prewhite = FALSE))

t test of coefficients:

            Estimate Std. Error t value  Pr(>|t|)    
(Intercept) 1.014312   0.148021  6.8525 8.975e-11 ***
x           0.471439   0.090756  5.1946 5.079e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

import pandas as pd
import statsmodels.api as sm

ts = pd.read_csv("data/hac_demo_timeseries.csv")
# maxlags here is the rule-of-thumb L; statsmodels has no Andrews auto-bandwidth,
# so for data-driven L compute floor(4 * (T/100)**(2/9)) explicitly or use R's
# NeweyWest(). use_correction=True applies the small-sample dof adjustment.
reg = sm.OLS.from_formula("y ~ x", data=ts).fit(
    cov_type="HAC",
    cov_kwds={"maxlags": 4, "use_correction": True}
)
print(reg.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.106
Model:                            OLS   Adj. R-squared:                  0.102
Method:                 Least Squares   F-statistic:                     26.71
Date:                Tue, 18 Aug 2026   Prob (F-statistic):           5.74e-07
Time:                        14:55:51   Log-Likelihood:                -334.25
No. Observations:                 200   AIC:                             672.5
Df Residuals:                     198   BIC:                             679.1
Df Model:                           1                                         
Covariance Type:                  HAC                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      1.0143      0.149      6.818      0.000       0.723       1.306
x              0.4714      0.091      5.169      0.000       0.293       0.650
==============================================================================
Omnibus:                        3.354   Durbin-Watson:                   0.706
Prob(Omnibus):                  0.187   Jarque-Bera (JB):                3.315
Skew:                           0.271   Prob(JB):                        0.191
Kurtosis:                       2.677   Cond. No.                         1.06
==============================================================================

Notes:
[1] Standard Errors are heteroscedasticity and autocorrelation robust (HAC) using 4 lags and with small sample correction

import delimited "data/hac_demo_timeseries.csv", clear
tsset t

newey y x, lag(4)

using GLM, DataFrames, CSV, CovarianceMatrices

ts    = CSV.read("data/hac_demo_timeseries.csv", DataFrame)
model = lm(@formula(y ~ x), ts)

# Newey-West (Bartlett kernel, 4 lags)
stderror(model, HAC(BartlettKernel(), 4))

I have [time series / panel] data and suspect serial correlation
in my regression errors.

Model:        [dependent variable] ~ [predictors]
Data:         T = [time periods], N = [units if panel]
Structure:    [pure time series / balanced panel / unbalanced panel]
Software:     [R / Python / Stata / Julia]

1. Should I use HAC (Newey-West) or cluster by time period?
   What are the trade-offs given my T = [T]?
2. How do I choose the lag length L? Provide the rule-of-thumb
   calculation for my T and show the data-driven alternative.
3. Provide the code to implement HAC SEs and, if appropriate,
   two-way clustering (by unit and by time).

Bootstrap Standard Errors

When parametric corrections are unavailable, unreliable, or difficult to justify on distributional grounds, bootstrap standard errors are the fallback. Rather than deriving a formula for the sampling distribution analytically, the bootstrap approximates it empirically by resampling the data many times and observing how the estimate varies. Two variants cover the vast majority of applied situations.

for b in 1..B:
    resample the data        # pairs: draw rows w/ replacement
                             # wild:  perturb residuals by w_i
    refit the model
    store beta_b

SE = std(beta_1, ..., beta_B)

# simple percentile 95% CI
# let beta_(1) <= ... <= beta_(B) be the sorted estimates (1-based order statistics)
sort beta_1, ..., beta_B
ci_lower = beta_( floor(0.025 * (B+1)) )    # B = 999 -> 25th order statistic
ci_upper = beta_( ceil( 0.975 * (B+1)) )    # B = 999 -> 975th order statistic

library(boot)
library(fwildclusterboot)

nunn <- haven::read_dta("data/Nunn_Wantchekon_AER_2011.dta")

# Keep n constant across resamples by dropping rows with missing
# values in the model variables before resampling
nunn <- nunn[complete.cases(nunn[, c("trust_neighbors", "ln_exports")]), ]

# Pairs bootstrap: statistic function must return estimate AND variance
# for studentized (percentile-t) CIs, which have better asymptotic properties
boot_fn <- function(data, idx) {
  fit <- lm(trust_neighbors ~ ln_exports, data = data[idx, ])
  c(
    coef(fit)["ln_exports"],
    vcov(fit)["ln_exports", "ln_exports"]   # variance needed for studentized CI
  )
}
set.seed(42)
boot_out <- boot(data = nunn, statistic = boot_fn, R = 1999)

# Studentized (percentile-t) CIs; preferred, with O(n^{-1}) coverage error
boot.ci(boot_out, type = "stud")

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 1999 bootstrap replicates

CALL : 
boot.ci(boot.out = boot_out, type = "stud")

Intervals : 
Level    Studentized     
95%   (-0.0033,  0.0093 )  
Calculations and Intervals on Original Scale

# Simple percentile CIs shown for comparison
boot.ci(boot_out, type = "perc")

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 1999 bootstrap replicates

CALL : 
boot.ci(boot.out = boot_out, type = "perc")

Intervals : 
Level     Percentile     
95%   (-0.0034,  0.0093 )  
Calculations and Intervals on Original Scale

# Wild cluster bootstrap (few clusters)
# fwildclusterboot >= 0.13: seed is set globally, not passed to boottest()
regression <- lm(trust_neighbors ~ ln_exports, data = nunn)
set.seed(42)
boottest(regression, param = "ln_exports", B = 1999, clustid = "isocode")

import numpy as np
import pandas as pd
import statsmodels.api as sm

nunn = pd.read_stata("data/Nunn_Wantchekon_AER_2011.dta")
nunn = nunn.dropna(subset=["trust_neighbors", "ln_exports"])
y = nunn["trust_neighbors"].values
X = sm.add_constant(nunn[["ln_exports"]].values)  # column 0 = intercept, 1 = ln_exports

rng = np.random.default_rng(42)
B = 1999
boot_coefs = np.zeros(B)

for b in range(B):
    idx = rng.integers(0, len(y), size=len(y))
    boot_coefs[b] = sm.OLS(y[idx], X[idx]).fit().params[1]  # index 1 = ln_exports

print(f"Bootstrap SE:  {boot_coefs.std():.4f}")

Bootstrap SE:  0.0030

# Simple percentile CI: easier to implement but O(n^{-1/2}) coverage error.
# For studentized (percentile-t) CIs with O(n^{-1}) accuracy, refit the model
# on each resample, record its SE, and standardize before taking quantiles.
print(f"95% CI: [{np.percentile(boot_coefs, 2.5):.4f}, "
      f"{np.percentile(boot_coefs, 97.5):.4f}]")

95% CI: [-0.0033, 0.0088]

use "data/Nunn_Wantchekon_AER_2011.dta", clear

* Pairs bootstrap
bootstrap _b[ln_exports], reps(1999) seed(42): ///
  regress trust_neighbors ln_exports

* Wild cluster bootstrap (Roodman et al.)
* ssc install boottest
boottest ln_exports, cluster(isocode) reps(1999) seed(42)

using GLM, DataFrames, ReadStataTables, Bootstrap

nunn = DataFrame(readstat("data/Nunn_Wantchekon_AER_2011.dta"))

fit_fn(data) = coef(lm(@formula(trust_neighbors ~ ln_exports), data))[2]

bs = bootstrap(fit_fn, nunn, BasicSampling(1999))
confint(bs, BasicConfInt(0.95))

I need bootstrap standard errors for my regression model.

Model:       [dependent variable] ~ [predictors]
n =          [sample size]
Estimator:   [OLS / quantile / IV / other]
Clusters:    [none / G = ? clusters]
Software:    [R / Python / Stata / Julia]

1. Should I use pairs bootstrap or wild bootstrap given my
   estimator and data structure?
2. If I have G = [N] clusters, should I use the wild cluster
   bootstrap instead of cluster-robust SEs?
3. How many replications B do I need for stable SE estimates
   vs. stable confidence intervals?
4. Provide the code to implement it and report both bootstrap
   SEs and percentile-t confidence intervals.

How much does my choice matter?

With so many possible design choices surrounding standard errors, it can sometimes be difficult to know how important this decision actually is. How much do our results change as a result of these choices? To demonstrate how much your SE choice affects your results, we set up a simulation experiment.

The rules of our experiment are relatively straightforward. We created 2,000 datasets in which individuals were assigned to \(G = \{6, 12, 30, 50\}\) moderately correlated clusters (\(\rho = 0.3\)) with 30 observations each. We then assigned a treatment effect to these clusters and set the average treatment effect to zero. Finally, we fit a regression to each dataset to test whether our model correctly predicts a treatment effect of zero. We use a \(5\%\) significance level as our threshold, meaning that a trustworthy result should only raise a false positive about \(5\%\) of the time.

The results from our simulation experiment are shown above in Table 1. Had we used OLS or heteroskedasticity-robust SEs, we would have rejected the null in roughly \(50\%\) of cases, far above our target of 0.05. This result is not entirely surprising, given that our treatment is determined by cluster, rather than by individuals; nevertheless, it illustrates the dangers of ignoring group dependencies. The plain cluster-robust estimator does better than HC SEs, except when there are fewer clusters. CR2 and the wild cluster bootstrap perform best in our simulation. Both are essentially correct by \(G = 12\) and remain relatively close to the target even when the number of clusters is small (\(G = 6\)).

library(clubSandwich); library(fwildclusterboot); library(sandwich); library(lmtest)

# One replication of a cluster-randomized design under a TRUE null (beta = 0),
# so every rejection at the 5% level is a false positive.
sim_once <- function(G, n_g = 30, rho = 0.3, B = 999) {
  g  <- rep(seq_len(G), each = n_g)
  D  <- rep(rbinom(G, 1, 0.5), each = n_g)          # cluster-level treatment
  u  <- rep(rnorm(G, 0, sqrt(rho)), each = n_g)     # cluster random intercept
  v  <- rnorm(G * n_g, 0, sqrt(1 - rho))            # idiosyncratic error
  y  <- u + v                                       # beta = 0: the null is true
  dat <- data.frame(y, D, g = factor(g)); fit <- lm(y ~ D, data = dat)
  c(
    OLS  = coef(summary(fit))["D", 4],
    HC1  = coeftest(fit, vcov = vcovHC(fit, "HC1"))["D", 4],
    CRVE = coeftest(fit, vcov = vcovCL(fit, cluster = ~g, type = "HC1"),
                    df = G - 1)["D", 4],                       # t(G-1)
    CR2  = coef_test(fit, vcov = "CR2", cluster = dat$g)$p_Satt[2],
    Wild = boottest(fit, param = "D", clustid = "g",
                    B = B, type = "rademacher")$p_val          # null imposed
  ) < 0.05
}

# Each cell: rejection rate over M draws; Monte Carlo SE = sqrt(p(1 - p) / M).
# Full sweep, seeds, and CSV output live in scripts/se_monte_carlo.R.

What Corrected SEs Do and Do Not Guarantee

Robust standard errors come with a specific promise. They are consistent in the sense that as the sample grows, they converge to the true sampling variability, but they are not unbiased. Because our standard errors are built from fitted residuals, and these residuals are smaller than the true errors they stand in for, standard errors in finite samples tend to run a little small. This small sample bias is the reason why this guide carries so many corrections. HC3, CR2, the \(t_{G-1}\) reference distribution, the wild cluster bootstrap, Webb weights, etc., all function as a fix for this downward bias. The smaller the sample, the more this matters.

If your robust standard errors come out far larger than the ordinary OLS ones, treat the gap itself as information. King and Roberts (2015) show that under a correctly specified model the two should converge as the sample grows, so a large divergence is a sign the model is misspecified, not merely that the errors are heteroskedastic. A robust standard error patches the symptom. Before reporting it, check whether the outcome needs a transformation, whether the functional form fits, or whether a different model (a count or binary model, say) suits the data better. When robust and OLS standard errors agree closely, reporting the robust version is reassuring. When they diverge sharply, look at the model first.

How can I test if spatial correlation is affecting my results?

The standard formal test is Moran’s I, which asks whether residuals at nearby locations are more similar to one another than would be expected by chance. In R, spdep::moran.test() computes this after you define a spatial weights matrix. A significant result suggests spatially autocorrelated errors. Plotting the residuals on a map is often just as informative; geographic clustering should be visually apparent if it is strong enough to matter.

Theoretical knowledge of your data frequently tells you as much as the test. If you are analyzing county-level health outcomes, district-level agricultural yields, or neighborhood housing prices, spatial correlation is almost certainly present. The test is most useful when you are genuinely uncertain, not as a post-hoc justification for a correction you were already planning to apply.

What’s the difference between Conley standard errors and spatial clustering?

Both account for spatial correlation, but through different mechanisms. Spatial (regional) clustering groups observations into discrete geographic units (countries, provinces, watersheds) and allows arbitrary within-group correlation with no assumptions about its structure. Conley SEs allow correlation to decay continuously with geographic distance using a kernel function, treating nearby observations as more correlated than distant ones rather than assuming a hard within-group/between-group divide.

In practice, if your spatial correlation is well captured by the administrative boundaries you already have, clustering is simpler, requires no coordinate data, and is easier to explain. Conley SEs earn their complexity when the relevant correlation is genuinely continuous, as when a shock propagates geographically rather than stopping at a border. The two approaches should give similar results when the clusters are small and dense; they can diverge when clusters are large or when observations near a cluster boundary are more similar to observations in an adjacent cluster than to those on the far side of their own.

How do I choose the bandwidth (lag length) for Newey-West standard errors?

There is no universally correct answer, but a common starting point is \(L = \lfloor 4(T/100)^{2/9} \rfloor\), which gives roughly four lags for \(T = 100\) observations and grows slowly with the sample size. R’s NeweyWest() offers data-driven selection via Andrews (1991) as an alternative; Stata’s base newey requires an explicit lag() (see the HAC section).

The most defensible practice is to report results across several bandwidths (say \(L = 2, 4, 6, 8\)) and assess sensitivity. If your key coefficient remains significant and stable in magnitude across the range, the bandwidth choice is not driving your conclusions. If results flip with a different \(L\), that reflects genuine uncertainty about the temporal dependence structure and is worth acknowledging directly.

In panel data, should I cluster by unit, by time, or both?

The answer follows the same logic as all SE decisions: cluster where the dependence was introduced. If treatment varies at the unit level and units are independent of each other, cluster by unit. If a common shock (a policy reform, a financial crisis, a pandemic) affects all units in a given period, you also have cross-sectional dependence and need to account for it. Two-way clustering by unit and by time handles both sources simultaneously, but it is not an automatic default: cluster where the design introduced dependence, and add the time dimension only when a genuine common shock requires it.

The practical cost of two-way clustering is reduced precision: the SEs will be larger, sometimes substantially. With very few time periods (\(T < 10\)), clustering by time produces too few clusters to be reliable; consider time fixed effects or period dummies as an alternative to soaking up the common-shock variation rather than correcting for it in the SE. fixest::feols() in R handles all of these cases cleanly and computes the correct small-sample adjustments.

I have an unbalanced panel. How does this affect standard error corrections?

Cluster-robust SEs handle unbalanced panels naturally. The cluster-sandwich estimator sums within-cluster score contributions regardless of how many observations each cluster contributes, so an unbalanced panel is not a problem; ensure your cluster variable correctly identifies which observations belong together.

HAC estimators are more sensitive. The autocovariance terms \(\hat{\Gamma}_j\) are harder to estimate reliably when the time dimension varies across units, particularly if missingness is non-random. If certain units are missing periods for structural reasons (attrition, selective reporting, entry and exit), the missing data problem may be more consequential than the SE correction itself, and it warrants its own treatment before addressing temporal autocorrelation. When the panel is only mildly unbalanced due to incidental missingness, fixest handles the HAC computation correctly.

My robust standard errors are much larger than my OLS standard errors. Should I just report the robust ones?

Not automatically. The divergence itself is informative. King and Roberts (2015) argue that a large gap between robust and OLS standard errors is a diagnostic signal of model misspecification: if the model were correctly specified, the two estimators would converge as the sample grows. A large discrepancy means that heteroskedasticity is severe enough to suggest the mean function may be wrong, not merely that the error variance is unequal.

The right response depends on the size of the divergence. If the two are close, robust SEs provide a reliable correction and you can report them with confidence. If they differ substantially, particularly if the change flips a result from significant to not, the better path is to re-examine the model before choosing which SEs to report. Common fixes include transforming the outcome or predictors (log, square root, or an inverse hyperbolic sine for skewed variables), revisiting the functional form, or switching to a generalized linear model suited to the data-generating process. Robust SEs are a valid correction for a specific problem, not a general substitute for a well-specified model.

How should I report standard errors, and how do I check that my choice isn’t driving my conclusions?

Reporting. Every regression table should include a note that states explicitly which standard errors are used and why (Imbens 2021). Not just a label like “robust SEs” in a footnote, but a sentence that connects the choice to the research design. A good note looks something like: “Standard errors clustered by state, the level at which the policy variable varies.” Readers should be able to evaluate your SE choice without hunting through the text. If the table uses heteroskedasticity-robust SEs, note which variant (HC2, HC3) and the software defaults you relied on. If you cluster, state the cluster variable and report the number of clusters \(G\) somewhere in the table or caption; \(G\) is critical context for judging whether cluster-robust inference is reliable.

Robustness checking. Whenever your SE choice involves a judgment call, run the regression under alternative specifications and report the results. The most common comparisons are worth making explicit:

• If you cluster, also show results with heteroskedasticity-robust SEs and, if the number of clusters is small, with the wild cluster bootstrap. If your conclusions change between these, that is important information.
• If you use a particular bandwidth for HAC or Conley SEs, show sensitivity across a range of values (e.g., \(h = 250, 500, 1000\) km for Conley; \(L = 2, 4, 8\) for Newey-West).
• If you are uncertain whether to cluster at a higher or lower level (e.g., county vs. state), report both. The more conservative choice should generally be the headline.

The goal is not to find the specification that produces significance; it is to show that your key inferences are stable across defensible alternatives. A result that holds under OLS, HC3, and clustered SEs is more credible than one that requires a specific correction to remain significant. If clustering absorbs the significance, that is a finding about the data structure, not a reason to avoid clustering.

My cluster count is in the 20 to 50 range and my design is unbalanced. Is standard CRVE still reliable?

Possibly not. Cluster-robust SEs lean on the idea that your clusters each contribute a roughly comparable amount of information. That breaks down whenever a few clusters do most of the work. Some common skewed designs:

• A difference-in-differences study where only 3 of 40 states ever adopt the policy, so the effect is really identified off a handful of treated clusters.
• One cluster that dwarfs the rest, as in a national survey where a single large state or city contributes a big share of the sample.
• A binary treatment with many control clusters but only a few treated ones.

In each case the effective number of clusters is far smaller than the raw count \(G\), so CRVE can under-cover even when \(G\) is 20, 30, or 50. The Bell-McCaffrey diagnostic puts a number on this. Read \(K_{BM}\) as “how many clusters are actually doing the work” for your coefficient;⁴ you never compute it by hand, since clubSandwich::coef_test() reports it for you (Bell and McCaffrey 2002; Imbens and Kolesar 2016).

The fix is to use CR2 corrections (Imbens and Kolesar 2016), which adjust the meat of the sandwich estimator to account for leverage heterogeneity across clusters. clubSandwich in R implements CR2 and reports the effective degrees of freedom automatically via coef_test().

library(clubSandwich)

fit <- lm(trust_neighbors ~ ln_exports, data = nunn)

# CR2 (Bell-McCaffrey) with effective degrees of freedom, in one call
coef_test(fit, vcov = "CR2", cluster = nunn$isocode)

Compare the effective degrees of freedom reported by coef_test() to \(G - 1\). A large gap (the \(K_{BM} < G/3\) flag from above) is a signal that the wild cluster bootstrap would be more reliable than CRVE alone.

I observe most or all of the relevant clusters, not just a sample. Does that change how I compute standard errors?

Yes. Standard CRVE is derived under the assumption that your clusters are a random sample from a large population of clusters, so the uncertainty comes from which clusters happened to be drawn. When you observe a large fraction of all clusters, such as all 50 US states or all 16 countries in Nunn and Wantchekon, that source of sampling uncertainty largely disappears. The relevant uncertainty is now over the assignment mechanism (which clusters or units were treated), not over which clusters were sampled. In this regime, where treatment assignment rather than cluster sampling is the operative source of randomness, conventional CRVE tends to be conservative and heteroskedasticity-robust SEs anti-conservative, so neither is generally well calibrated. The exact direction and size of the gap depend on the estimand and on how much treatment effects vary across clusters (Abadie et al. 2023); the design-based framing below is what makes “too wide” or “too narrow” well defined.

Abadie et al. (2023) address this with two methods that target the actual design variance. CCV (Causal Cluster Variance) computes the variance directly from the potential-outcomes framework; it can be unstable in small samples. TSCB (Two-Stage Cluster Bootstrap) is generally preferred: Stage 1 resamples the cluster-level treatment assignments, holding the clusters themselves fixed; Stage 2 resamples units within clusters. The two-stage structure preserves the within-cluster correlation while only randomizing what is actually uncertain: the treatment assignment.

The exact CCV is computed from the potential-outcomes decomposition, but a convenient approximation you can build today interpolates between two matrices you already know how to get, the cluster-robust and the heteroskedasticity-robust, by the sampling fraction \(\rho\) (the share of each cluster’s target population you observe). As \(\rho \to 1\), design uncertainty dominates and the variance tends to the heteroskedasticity-robust matrix (CRVE is conservative there); as \(\rho \to 0\) it returns the cluster-robust matrix:

# Convex approximation to CCV (Abadie et al. 2023); the exact CCV comes from the
# potential-outcomes decomposition and may differ. rho = sampling fraction in (0, 1].
library(sandwich)
fit       <- lm(y ~ x, data = df)
V_cluster <- vcovCL(fit, cluster = df$g)   # Liang-Zeger
V_robust  <- vcovHC(fit, type = "HC2")     # heteroskedasticity-robust
V_ccv     <- V_robust + (1 - rho) * (V_cluster - V_robust)
se_ccv    <- sqrt(diag(V_ccv))["x"]        # use t(G-1) critical values

import numpy as np, statsmodels.formula.api as smf
fit_cl = smf.ols('y ~ x', data=df).fit(cov_type='cluster', cov_kwds={'groups': df['g']})
fit_hc = smf.ols('y ~ x', data=df).fit(cov_type='HC2')
rho    = 0.1
V_ccv  = fit_hc.cov_params() + (1 - rho) * (fit_cl.cov_params() - fit_hc.cov_params())
se_ccv = np.sqrt(np.diag(V_ccv))[list(fit_hc.params.index).index('x')]

TSCB does not have a one-line closed form, but the algorithm is fully specified, so you can implement it directly rather than cloning a private repository:

# Two-Stage Cluster Bootstrap (Abadie et al. 2023). Pseudocode, not runnable.
# ASSUMPTION: per-cluster sampling fraction rho is known/defensible.
#   rho ~ 0 -> ordinary cluster bootstrap; rho ~ 1 -> design uncertainty dominates.
# REQUIRES within-cluster treatment variation (treated AND control units per cluster).
for b in 1..B:
    # Stage 1: resample clusters (the design-uncertainty step).
    #   Stage 1 must mirror the ACTUAL assignment design (fixed # treated
    #   clusters, blocking/strata); a uniform permutation is valid only
    #   under complete cluster randomization.
    draw G clusters
    for each drawn cluster g:
        treated_g, control_g = units in g with d==1, d==0
        if treated_g empty OR control_g empty: skip g   # guard: split undefined
        # Stage 2: resample units within g WITH REPLACEMENT, preserving the
        #          treated/control split (resample treated and control separately)
        append resampled units to bootstrap sample
    if usable clusters < 2: discard replicate            # guard
    refit the model on the bootstrap sample; store the coefficient of interest
SE = sd(stored coefficients)
CI = percentile interval (2.5%, 97.5%) of stored coefficients

The CCV approximation above needs no external code. For TSCB there is no single base-package flag in fwildclusterboot; consult the Abadie et al. (2023) replication materials for a reference implementation, though the algorithm above is enough to reproduce the method directly.

Use CCV/TSCB when (a) the treatment varies at the cluster level, (b) you observe a substantial fraction of all clusters, and (c) the clusters were not themselves sampled at random. In observational studies covering a near-census of clusters, TSCB often yields narrower intervals than CRVE, reflecting the reduced uncertainty when the sampling step is largely absent. This is not guaranteed: the comparison depends on the estimand and the design, and where treatment-effect heterogeneity is substantial the design-based interval can be as wide as, or wider than, CRVE.

When should I use randomization inference instead of a sandwich or bootstrap?

When the randomness you care about is in assignment rather than sampling, or when you have too few clusters for any large-\(G\) method to be trusted, randomization inference (RI) gives exact finite-sample tests. Instead of assuming a sampling distribution, you re-run the study’s actual assignment mechanism many times under the sharp null of no effect and see where your observed statistic falls. Two cautions: the default null is sharp (zero effect for every unit, not just on average), and the permutation must mirror the real design, blocking, stratifying, or permuting at the cluster level exactly as treatment was assigned. Confidence intervals come from inverting the test over a grid of hypothesized effects \(\tau_0\): for each \(\tau_0\), subtract its implied effect from the outcome, test the sharp null, and keep the values not rejected at 5%.

library(ri2); library(randomizr)
decl <- declare_ra(clusters = df$cluster_id, m = treated_clusters)  # the ACTUAL design
grid <- seq(-2, 2, by = 0.02)
keep <- sapply(grid, function(tau0) {
  df$y_adj <- df$y - tau0 * df$z
  summary(conduct_ri(y_adj ~ z, declaration = decl,
                     assignment = "z", sharp_hypothesis = 0, data = df))$two_tailed_p_value > 0.05
})
ci_ri <- range(grid[keep])   # the set of tau0 not rejected at 5%

In Stata, ritest does this directly; Python has no mature finite-population RI package, so code the permutation loop by hand or run this step in R or Stata.

How do standard errors work for IV / 2SLS?

The asymptotic sandwich logic carries over: heteroskedasticity-robust and cluster-robust SEs are valid for two-stage least squares, computed from the first-stage-projected regressors (cluster at the design level). The finite-sample leverage corrections do not transfer mechanically, though, HC2/HC3 and CR2 are defined for OLS, so use an IV-specific implementation (clubSandwich supports CR2 for ivreg) rather than assuming the OLS weights apply. The bigger catch is weak instruments. When the first stage is weak, the usual 2SLS standard errors and \(t\)-tests are unreliable no matter how you compute the variance, because the estimator is not even approximately normal. Report a first-stage \(F\), and when it is low switch to weak-instrument-robust inference, the Anderson–Rubin test or the conditional likelihood ratio (CLR) test, which stay valid regardless of instrument strength. In R, estimate with fixest::feols(y ~ controls | endog ~ instr) or AER::ivreg() and run AR/CLR via ivmodel; in Stata, ivregress 2sls, vce(cluster …) plus weakiv or ivreg2; in Python, linearmodels.IV2SLS, though weak-instrument-robust inference is limited there, so prefer R ivmodel or Stata weakiv for that step.

How do I get robust standard errors for logit, Poisson, and other non-linear models?

The same sandwich, with a different bread. For a maximum-likelihood estimator the robust variance is \(\mathbf{A}^{-1}\mathbf{B}\mathbf{A}^{-1}\), where \(\mathbf{A} = -\mathbb{E}[\partial^2\ell/\partial\theta\partial\theta^\top]\) is the (expected) negative Hessian of the log-likelihood, i.e. the information matrix, and \(\mathbf{B}\) is the variance of the score; the cluster-robust version sums scores within clusters, exactly as for OLS. In R, apply sandwich::vcovHC() or sandwich::vcovCL() to a fitted glm object; in Python, pass cov_type="HC0" or cov_type="cluster" to the statsmodels .fit(); in Stata, add vce(robust) or vce(cluster …) to logit/poisson/glm. One extra step matters in non-linear models: the quantity you actually report, a marginal effect, a predicted probability, an odds ratio, is a non-linear function of the coefficients, so propagate its uncertainty with the delta method rather than reading a coefficient SE directly, using marginaleffects::avg_slopes() in R, .get_margeff() in statsmodels, or margins in Stata. The next question generalizes this delta-method step to any nonlinear quantity. (For GMM estimators the structure is identical, with the optimal-weighting-matrix meat; see R’s gmm package or Stata’s gmm.)

How do I get a standard error for a function of my coefficients?

Often the number you want to report is not a coefficient itself but a function of several: a ratio of two coefficients, an incidence-rate or odds ratio, an elasticity, a marginal effect, or the long-run multiplier in a dynamic model. The delta method gives its standard error. If \(g(\boldsymbol{\theta})\) is the quantity of interest and \(\hat{\mathbf{V}}\) is the variance matrix of \(\hat{\boldsymbol{\theta}}\), then

\[\mathrm{Var}\big(g(\hat{\boldsymbol{\theta}})\big) \;\approx\; \nabla g(\hat{\boldsymbol{\theta}})^\top\, \hat{\mathbf{V}}\, \nabla g(\hat{\boldsymbol{\theta}}),\]

where \(\nabla g\) is the gradient of \(g\) evaluated at the estimates. You linearize the transformation around \(\hat{\boldsymbol{\theta}}\) and push the coefficient uncertainty through it.

Two points keep this reliable. First, feed it the robust \(\hat{\mathbf{V}}\), the HC, cluster, or HAC matrix you would otherwise report, so the transformed quantity inherits the same correction; the delta method does not add robustness, it only propagates whatever variance you hand it. Second, the linear approximation degrades when \(g\) is sharply nonlinear near the estimates or the sample is small, the textbook case being a ratio whose denominator sits close to zero. There the bootstrap, resample, recompute \(g\), and read the interval off the distribution, is the safer choice and needs no derivatives.

You rarely differentiate by hand. In R, marginaleffects (avg_slopes(), predictions(), and hypotheses() for arbitrary functions) and car::deltaMethod() take the fitted model and a robust vcov; in Python, statsmodels exposes .get_margeff() and t_test/NonlinearDeltaCov; in Stata, margins and nlcom compute delta-method SEs for nonlinear combinations directly.

How do I test for serial correlation in a panel?

The Wooldridge test for autocorrelation in panel data is the standard quick check: plm::pwartest() in R (or plm::pbgtest() for a panel Breusch–Godfrey), and xtserial in Stata. A rejection says within-unit errors are serially correlated, which calls for clustering by unit (or Driscoll–Kraay when common time shocks are also present). As with the cross-sectional diagnostics, treat the test as confirmation rather than permission: in most panels with repeated measures on the same units you should cluster by unit regardless of the result.

Fixed or random effects, does the choice change my standard errors?

Fixed-versus-random effects is a question about identification, whether the unit effects may be correlated with your regressors (a Hausman test is the usual arbiter), not about standard errors. Under either choice you should still cluster by the unit whose repeated measures create the dependence, so the SE machinery in this guide applies to both the same way.

I’m running a staggered difference-in-differences design. What about the standard errors?

Fix the estimator before you touch the standard errors. When units adopt treatment at different times, the conventional two-way fixed-effects (TWFE) estimator is biased for the average treatment effect under heterogeneous or dynamic effects, because it implicitly uses already-treated units as controls; no SE correction repairs a biased point estimate (Roth et al. 2023). Use a heterogeneity-robust estimator instead: Callaway and Sant’Anna (Callaway and Sant’Anna 2021), Sun and Abraham (Sun and Abraham 2021), or the Borusyak–Jaravel–Spiess imputation estimator (Borusyak, Jaravel, and Spiess 2024). Their inference is then standard, cluster at the unit level (the level treatment is assigned and serial correlation lives), and they report cluster-robust or multiplier-bootstrap SEs natively. Event-study (dynamic) SEs are the same machinery applied period by period. In R: did (Callaway–Sant’Anna), fixest::sunab() (Sun–Abraham), didimputation; in Stata: csdid, eventstudyinteract, did_imputation; in Python: differences (Callaway–Sant’Anna) and pyfixest (Sun–Abraham) are usable but younger, so verify the inference options before relying on them.

I have very few clusters and cannot get more. What do I actually report?

This is a hard case, and there is no estimator that manufactures information that isn’t there. With a handful of clusters (\(G\) below about 10) and a fixed design, work down this list: (1) the wild cluster bootstrap with Webb six-point weights and the null imposed, which is the best-calibrated option in this regime (MacKinnon, Nielsen, and Webb 2023; Cameron, Gelbach, and Miller 2008); (2) randomization inference, if the assignment mechanism is known and you are willing to test a sharp null; (3) if neither applies, report the interval and stop treating the \(p\)-value as meaningful, and note that with so few clusters the test has little power and wide intervals are the correct conclusion, not a defect. What you should not do is cluster at a finer level to manufacture clusters, or fall back to OLS/HC SEs that ignore the design; both understate the true uncertainty. Pooling with related units to raise \(G\), if defensible, is the only way to genuinely narrow the interval.

How do standard errors work after matching or propensity-score matching?

Naively, people report the regression or weighted-mean SE from the matched sample as if the matching never happened, which is wrong: matching introduces its own sampling variability, and the matched units are not independent draws. Abadie and Imbens (Abadie and Imbens 2006) derived the correct large-sample variance for matching estimators, and it is what you should report. In R, Matching::Match() computes Abadie–Imbens SEs natively, and MatchIt pairs well with marginaleffects for the post-matching estimate; in Stata, teffects nnmatch provides the Abadie–Imbens SEs natively (psmatch2 is widely used but does not report the AI variance by default). Python tooling here is thin and largely unmaintained, so for matching standard errors specifically, prefer R or Stata. (Note that nearest-neighbor matching SEs are not bootstrap-consistent, so do not reach for the pairs bootstrap as a shortcut.)

I’m testing many outcomes or subgroups. How do I adjust the standard errors?

Standard errors themselves do not change, but the inferential threshold must, or the family-wise error rate inflates with the number of tests. Bonferroni and Benjamini–Hochberg are the familiar blunt instruments; the sharper tool for applied work is the Romano–Wolf stepdown procedure (Romano and Wolf 2005), which controls the family-wise error rate while exploiting the dependence across your test statistics (via the bootstrap), so it is far less conservative than Bonferroni when outcomes are correlated. In R use wildrwolf; in Stata, rwolf or wyoung. Python offers only the independence-based corrections (statsmodels.stats.multitest), not Romano–Wolf, so run that step in R or Stata. Decide the family of tests before looking at the results, and report both adjusted and unadjusted \(p\)-values.

My standard errors look fine, but I’m worried about omitted variables. Can SEs help?

No, standard errors quantify sampling variability around the estimate your model implies; they say nothing about whether that estimate is biased by confounding, selection, or measurement error. A tight, correctly clustered confidence interval around a confounded estimate is precisely measured nonsense. The complement to a good standard error is a sensitivity analysis that asks how strong unobserved confounding would have to be to overturn your conclusion: Oster’s \(\delta\) for coefficient stability under proportional selection (Oster 2019) (robomit in R, psacalc in Stata), or the more general Cinelli–Hazlett framework (Cinelli and Hazlett 2020) (sensemakr, maintained in both R and Python). Report one of these alongside your SEs whenever you make a causal claim; it addresses the uncertainty your standard error deliberately ignores.

What are the most common standard-error mistakes?

Most errors are decisions, not typos. The recurring ones:

1. Clustering at the wrong level. Clustering finer than the design (individuals when treatment is at the school) understates uncertainty; not clustering at all when the design induces dependence does the same. Cluster where sampling or assignment introduced the correlation.
2. Treating a large robust-vs-OLS gap as a reason to “use robust.” It is usually a sign the mean function is misspecified; fix the model first (King and Roberts 2015).
3. Clustering at the firm-by-year intersection instead of using the additive two-way estimator. The intersection is a bias-correction term, not the estimator.
4. Trusting normal critical values with few clusters. Below ~30 clusters, use \(t_{G-1}\) and CR2 or the wild cluster bootstrap; the raw CRVE over-rejects.
5. Reporting naive SEs after matching, ignoring the variability the matching step introduces (use Abadie–Imbens).
6. Clustered TWFE in a staggered DiD. The estimator is biased before the SE is even computed; switch to a heterogeneity-robust estimator first.
7. Reading a tight standard error as evidence against confounding. Precision is not validity; pair causal claims with a sensitivity analysis.

Standard-error selection table

This condenses the interactive decision tree at the top of the guide into a static reference. Find the row that matches how your data were sampled and how treatment was assigned; each section above is an extended explanation of the corresponding cell. Where a language has no trustworthy off-the-shelf tool, that is stated rather than papered over. (The table scrolls horizontally on narrow screens.)

Assumptions and consistency by estimator

What separates these estimators is not the sandwich structure (shared by all; see the unification box in the HC section) but which dependence each one assumes away, the asymptotics it needs to be consistent, and how it behaves before those asymptotics kick in.

Every row above varies the sandwich meat (which correlations it admits). Non-linear and quasi-likelihood models, logit, Poisson, GMM, keep the same sandwich but change the bread: the inverse information matrix replaces \((\mathbf{X}^\top\mathbf{X})^{-1}\), and the robust SE there guards against a misspecified variance (e.g., Poisson overdispersion) rather than against a dependence structure. See the FAQ “How do I get robust standard errors for logit, Poisson, and other non-linear models?” for the formula and code, and the question after it for delta-method SEs of any nonlinear quantity.

Reproducibility

Every result, figure, and code block in this guide is reproducible from the materials in the repository. A few specifics:

• Software. Rendered with R 4.5.2 and Python 3.12.13. Key package versions: sandwich 3.1.1, lmtest 0.9-40, fixest 0.13.2, clubSandwich 0.6.2, fwildclusterboot 0.14.3, boot 1.3-32, plm 2.6-7, survey 4.5; statsmodels 0.14.6, numpy 2.4.6, pandas 3.0.3. The R environment is pinned in the committed renv.lock; the Python environment is the project virtual environment (./.venv), from which a requirements.txt can be exported with pip freeze for exact reproduction.
• Seeds. Every stochastic step (the bootstraps, the simulated HAC series, the animations) sets a seed explicitly; we use set.seed(42) for the bootstrap examples. Bootstrap and simulation results are reproducible bit-for-bit given the versions above.
• A version-sensitive gotcha. fwildclusterboot \(\ge\) 0.13 sets the bootstrap seed globally (set.seed() before the call) rather than through a seed = argument to boottest(). The pinned 0.14.3 follows this convention; older code that passes seed = will error.
• Data provenance. The running example uses the Nunn & Wantchekon (2011) replication file (data/Nunn_Wantchekon_AER_2011.dta); sources are linked in the Download the Data callout. The estimation sample is 16 countries, not the raw 17, because trust_neighbors is missing for all of Zimbabwe (see footnote 1).
• Figures. The animations are built by the scripts in scripts/ (e.g., se_comparison.R) and saved to images/; rerun a script to regenerate its figure.

Papers and Books

Software Documentation

Tutorials and Course Materials