OLS regression in Python with robust and clustered standard errors

Ordinary least squares (OLS) is the default tool for relating an outcome to a set of explanatory variables. The coefficient estimates come from one formula, but the standard errors depend on what you assume about the errors: whether their variance is constant, and whether observations in the same firm, school or period share shocks. A wrong assumption leaves the coefficients unchanged and makes every test and confidence interval misleading. This page shows how to fit OLS and weighted least squares in Python with oe.ols, how to choose the covariance estimator, and how to test hypotheses after the fit.

What OLS estimates

OLS chooses the coefficients that minimize the sum of squared residuals. With weights, each squared residual is multiplied by its weight; this is weighted least squares (WLS):

β̂ = (X′WX)⁻¹ X′Wy,   W = diag(w₁, …, wₙ), and W = I for OLS

The estimates describe the best linear approximation to the conditional mean of y. They have a causal reading only if the regressors are uncorrelated with the error term. When that fails, because of omitted variables, measurement error or simultaneity, no choice of standard errors repairs the bias; see the instrumental variables page.

The covariance estimator is a separate decision, and it follows from how the data were generated:

• Classical standard errors assume independent errors with a common variance.
• Heteroskedasticity-robust (HC) standard errors drop the common-variance assumption (White 1980). HC2 and HC3 add finite-sample corrections that matter in small samples (MacKinnon and White 1985).
• Cluster-robust standard errors allow any correlation among observations in the same cluster, with independence across clusters. Multiway clustering extends this to two or more non-nested dimensions, such as firm and year (Cameron, Gelbach and Miller 2011).
• HAC standard errors allow heteroskedasticity and serial correlation up to a chosen lag in time-series data.

Clustering matters most for a regressor that varies only at the group level, such as a firm-wide policy studied with worker data. If the errors share a group component, classical standard errors treat each worker as independent information about the policy and can understate the uncertainty severely (Moulton 1990). The example below shows that case.

The guide to robust and clustered standard errors explains each estimator's formula, compares them on one dataset and covers the case of few clusters.

In OpenEconometrics

oe.ols fits OLS and WLS. Give the outcome and regressors as y="lwage", x=["educ", "exper"], or as a formula such as formula="lwage ~ educ + exper + I(exper**2) + C(region)". Formulas support +, interactions : and *, I(), C() for categorical terms, log(), exp(), sqrt(), and the time operators L() and D() with time=. Collinear terms are omitted and listed in fit.provenance["omitted_terms"]. Rows with missing values are dropped by default (missing="drop").

The model in the example corresponds to the Stata command regress lwage educ c.exper##c.exper union, vce(cluster firm). The main keywords are:

• covariance="nonrobust" (the default), "HC0", "HC1", "HC2" or "HC3". "robust" is an alias for "HC1", the estimator Stata's vce(robust) uses.
• covariance="cluster" with cluster="firm" or a list of up to four columns. Passing cluster= alone selects cluster covariance. "cluster_hc2" and "cluster_hc3" add a leverage correction for one-way clustering.
• covariance="hac" with time= (an integer period column), lags= (a number or "auto") and kernel="bartlett", "parzen", "quadratic_spectral" or "truncated".
• covariance="bootstrap" with reps= and seed=, or covariance="jackknife", by observation or by one-way cluster.
• weights="w" with weight_type="aweight", "fweight", "pweight" or "iweight". Frequency weights are integer counts of identical rows. Sampling weights (pweight) require a robust or cluster covariance and default to HC1.
• dfadjust=True applies Bell–McCaffrey degrees of freedom to single contrasts with HC2, HC3 or their cluster versions; hansen=True applies Hansen's HC3 adjustment.

After the fit, oe.test runs Wald tests of one or several linear restrictions, oe.testparm tests groups of terms by name or pattern, oe.lincom and oe.nlcom estimate linear and nonlinear combinations, and oe.margins computes average marginal effects or effects at the means. Prediction and diagnostics are methods of the result: fit.predict(), fit.vif(), fit.hettest(), fit.white_test(), fit.reset_test() and fit.breusch_godfrey(). fit.to_latex() exports the coefficient table.

oe.regress and oe.newey are thin wrappers that accept Stata's vocabulary and return the same result as the equivalent oe.ols call. The OLS documentation gives every convention, and the linear family guide describes the wrappers.

Supported options

Generated from oe.capabilities() of openecon 0.3.20a1. Stata names identify the corresponding workflow, not numerical parity.

Example

The script simulates 3,000 workers in 60 firms of 50 workers each. Union coverage is set at the firm level, and each firm adds a shared wage shock with standard deviation 0.15. The individual error's standard deviation rises from 0.2 to 0.5 with experience. The true coefficients are 0.08 for a year of schooling and 0.10 for union coverage; the experience profile is 0.03 · exper − 0.0005 · exper², so the true return to experience at 10 years is 0.02.

import numpy as np
import pandas as pd
import openecon as oe


def simulate(n_firms=60, per_firm=50, seed=11):
    """Log wages of workers in firms: a firm-level shock and heteroskedastic individual errors."""
    rng = np.random.default_rng(seed)
    firm = np.repeat(np.arange(n_firms), per_firm)
    union = np.repeat(rng.integers(0, 2, size=n_firms), per_firm)  # set at the firm level
    firm_shock = np.repeat(rng.normal(scale=0.15, size=n_firms), per_firm)
    n = n_firms * per_firm
    educ = rng.integers(10, 19, size=n)
    exper = rng.uniform(0, 30, size=n)
    noise = rng.normal(size=n) * (0.2 + 0.01 * exper)  # error spread grows with experience
    lwage = (1.2 + 0.08 * educ + 0.03 * exper - 0.0005 * exper**2 + 0.10 * union
             + firm_shock + noise)
    return pd.DataFrame({"lwage": lwage, "educ": educ, "exper": exper,
                         "union": union, "firm": firm})


if __name__ == "__main__":
    df = simulate()
    formula = "lwage ~ educ + exper + I(exper**2) + union"
    for cov in ["nonrobust", "HC3", "cluster"]:
        fit = oe.ols(data=df, formula=formula, covariance=cov,
                     cluster="firm" if cov == "cluster" else None)
        se = {c.term: c.std_error for c in fit.coefficients}
        print(f"{cov:>9}: se(educ) = {se['educ']:.4f}, se(union) = {se['union']:.4f}")
    print()
    print(fit.summary())
    joint = oe.test(fit, ["exper", "I(exper**2)"])
    print(f"Experience terms jointly zero: F({joint['df_num']}, {joint['df_denom']}) = "
          f"{joint['statistic']:.2f}, p = {joint['p_value']:.3g}")
    at10 = oe.lincom(fit, {"exper": 1, "I(exper**2)": 20})
    print(f"Return to experience at 10 years: {at10['estimate']:.4f} "
          f"(se {at10['std_error']:.4f})")
    ame = oe.margins(fit, variables=["exper"], method="ame")
    print(f"Average marginal effect of experience: {ame['estimate'][0]:.4f} "
          f"(se {ame['std_error'][0]:.4f})")
    print(f"True average marginal effect: {0.03 - 0.001 * df['exper'].mean():.4f}")

With openecon 0.3.20a1 from PyPI the script prints:

nonrobust: se(educ) = 0.0028, se(union) = 0.0145
      HC3: se(educ) = 0.0028, se(union) = 0.0143
  cluster: se(educ) = 0.0031, se(union) = 0.0359

OLS — lwage
Observations: 3000  |  Covariance: cluster  |  Confidence: 95%

Term             Estimate   Std. error         t     P>|stat|      CI lower      CI upper
-----------  ------------  -----------  --------  -----------  ------------  ------------
Intercept         1.26941    0.0484872   26.1803  3.57833e-34       1.17239       1.36643
educ            0.0782193   0.00305974    25.564  1.30456e-33     0.0720968     0.0843418
exper           0.0321155   0.00302642   10.6117  2.65327e-15     0.0260596     0.0381713
I(exper**2)  -0.000566533  9.95005e-05  -5.69377  4.14203e-07  -0.000765633  -0.000367433
union           0.0482095    0.0358931   1.34314     0.184371    -0.0236125      0.120032

r_squared: 0.277119  |  adjusted_r_squared: 0.276154  |  rmse: 0.388306  |  df_model: 4
df_resid: 2995  |  ss_resid: 451.59  |  ss_model: 173.119  |  ss_total: 624.709
log_likelihood: -1416.43  |  condition_number_scaled: 7.8557  |  aic: 2842.85  |  bic: 2872.89
Model F test: F(4, 59) = 372.302, p = 4.057e-41
Experience terms jointly zero: F(2, 59) = 221.59, p = 3.67e-28
Return to experience at 10 years: 0.0208 (se 0.0012)
Average marginal effect of experience: 0.0151 (se 0.0008)
True average marginal effect: 0.0150

Reading the output

• Clustering changes the union standard error the most. It rises from 0.0145 (classical) to 0.0359 (clustered by firm), a factor of about 2.5. Union status and the wage shock are both shared by the 50 workers of a firm, so the data hold about 60 independent pieces of information on union coverage, not 3,000. The classical interval, roughly 0.02 to 0.08, would exclude the true 0.10.
• The clustered estimate is honest but imprecise. The union coefficient is 0.0482 with a 95% interval from −0.024 to 0.120. The interval contains 0.10, but with 60 firms the data cannot distinguish the effect from zero (p = 0.18).
• Heteroskedasticity barely moves these two standard errors. HC3 gives 0.0028 for schooling and 0.0143 for union, close to the classical values. The error variance grows with experience, which is unrelated to schooling and union here; robust and classical errors diverge when the variance is related to the regressor in question.
• Schooling is estimated precisely. The estimate is 0.0782 with an interval from 0.0721 to 0.0843, which contains 0.08. Clustered inference uses G − 1 = 59 degrees of freedom, which is why the tests show F(·, 59).
• Postestimation uses the clustered covariance. The experience terms are jointly significant (F(2, 59) = 221.59). oe.lincom puts the return to experience at 10 years at 0.0208 against a true 0.02. oe.margins differentiates through the squared term and averages over the sample, giving 0.0151 against a true 0.0150.

Limitations

• Cluster-robust covariance for oe.ols allows at most 4 clustering dimensions.
• HAC covariance needs an explicit integer time column; oe.ols never assumes that rows are consecutive periods. lags="auto" (the Newey–West plug-in rule) covers the Bartlett, Parzen and quadratic-spectral kernels; the truncated kernel needs an explicit lags value.
• Bootstrap and jackknife covariances resample observations or one-way clusters only; multiway resampling is rejected.
• dfadjust=True and hansen=True adjust single coefficients and contrasts; joint Wald tests keep their conventional degrees of freedom.
• Heteroskedasticity tests are available for unweighted and frequency-weighted models, and serial-correlation tests for unweighted models only.
• Absorbed high-dimensional fixed effects, restricted least squares and complex survey designs are separate functions, such as oe.reghdfe and oe.cnsreg.
• No numerical parity with Stata is claimed: the inventory's stata_parity_validated flag is False.

References

• Cameron, A. C., Gelbach, J. B. and Miller, D. L. (2011). Robust inference with multiway clustering. Journal of Business & Economic Statistics, 29(2), 238–249. doi:10.1198/jbes.2010.07136
• MacKinnon, J. G. and White, H. (1985). Some heteroskedasticity-consistent covariance matrix estimators with improved finite sample properties. Journal of Econometrics, 29(3), 305–325. doi:10.1016/0304-4076(85)90158-7
• Moulton, B. R. (1990). An illustration of a pitfall in estimating the effects of aggregate variables on micro units. Review of Economics and Statistics, 72(2), 334–338. doi:10.2307/2109724
• White, H. (1980). A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity. Econometrica, 48(4), 817–838. doi:10.2307/1912934