From 924a0ce0130e66fea8c19f9cafcfb8bebe04d4b1 Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Mon, 12 Oct 2026 15:04:42 +0100 Subject: [PATCH 1/4] BUG: Fix the degrees of freedom of the RandomEffects variance components The Swamy-Arora estimators of the variance of the idiosyncratic error and of the variance of the effects divided the residual sums of squares of the within and the between regressions by n - N - K + 1 and N - K, where K is the number of columns of exog. The regressions do not estimate a coefficient for every column. A regressor that is constant within entities (a constant, years of education) or that is collinear with the others once the means of the entities are removed (experience, when the model has year dummies) is not in the within regression, and the dummy variables of the years are collinear with the constant in the between regression of a balanced panel. The "+ 1" is also only right for a model with a constant, as the documentation of the estimator says. The variances were slightly wrong for these models, and so were theta, the estimates and the standard errors. The degrees of freedom are now the numbers of coefficients that the regressions estimate, the ranks of the within-demeaned regressors and of the means of the regressors. The trace in the estimator of Baltagi and Chang (1994) uses a pseudo inverse, since the means of the regressors can be collinear. A model with a constant and regressors that all vary within entities and are not collinear in the means, such as the models that are compared with Stata in the tests, does not change. For the wage panel and the model of the panel_examples notebook (year dummies; black, hisp and educ do not vary over time and exper is collinear with the year dummies) the variances are now 0.1231940 and 0.1053672 instead of 0.1233236 and 0.1069460, theta is 0.64291 instead of 0.64506, and the coefficient of union is 0.10613 instead of 0.10587. These agree with plm (ercomp and coef) to 7 digits. The Hausman statistic of the fixed effects and random effects models of the wage panel with year dummies is now 75.3108, which is the statistic of plm::phtest, instead of 74.5621. Stata's xtreg, re counts in the same way: the sigma_e of a model with a time-invariant regressor is the sigma_e of xtreg, fe, which omits it. In A. C. Cameron's log for the Cornwell and Rupert wage panel (https://cameron.econ.ucdavis.edu/bgpe2013/ct_panel.txt), xtreg lwage exp exp2 wks ed, re has sigma_u .31951859, sigma_e .15220316 and rho .81505521. RandomEffects had .31951715, .15222450 and .81501158, and now has .31951860, .15220316 and .81505522. The tests compare the idiosyncratic variance with the residual variance of the regression on the dummy variables of the entities, the estimators written out in NumPy, and plm (results/re-variance-components-reference.R). --- linearmodels/panel/model.py | 42 ++- .../re-variance-components-reference.R | 91 ++++++ .../re-variance-components-reference.csv | 73 +++++ .../panel/test_random_effects_variance.py | 281 ++++++++++++++++++ 4 files changed, 481 insertions(+), 6 deletions(-) create mode 100644 linearmodels/tests/panel/results/re-variance-components-reference.R create mode 100644 linearmodels/tests/panel/results/re-variance-components-reference.csv create mode 100644 linearmodels/tests/panel/test_random_effects_variance.py diff --git a/linearmodels/panel/model.py b/linearmodels/panel/model.py index b2c36d263c..d2b2b277d9 100644 --- a/linearmodels/panel/model.py +++ b/linearmodels/panel/model.py @@ -2982,6 +2982,25 @@ def fit( estimator similar to the Newey-West covariance estimator. * ``bandwidth`` - Bandwidth to use when computing the kernel. If not provided, a naive default is used. + + The variance of the idiosyncratic error is the residual variance of the + within regression, and the variance of the effects is estimated from the + residuals of the between regression. The degrees of freedom of the + variances are the numbers of coefficients that the regressions estimate. A + regressor that is constant within entities (the constant, or years of + education in a wage panel) or that is collinear with the others once the + means of the entities are removed (experience, if the model has year + dummies) is not in the within regression, and a regressor that is collinear + with the others in the means of the entities (the year dummies and the + constant in a balanced panel) is not in the between regression. + + .. versionchanged:: 7.1 + + The degrees of freedom of the variances count the coefficients that the + regressions estimate. Earlier versions counted every column of ``exog`` + (less one for the constant in the within regression), so that the + variances, theta and the estimates were slightly different for models with + time-invariant or collinear regressors, and for models without a constant. """ w = self.weights.values2d root_w = cast("linearmodels.typing.data.Float64Array", np.sqrt(w)) @@ -2991,6 +3010,10 @@ def fit( assert isinstance(demeaned_exog, PanelData) y = demeaned_dep.values2d x = demeaned_exog.values2d + # The regressors that do not vary within entities (a constant, a time-invariant + # regressor) or that are collinear once the means are removed have no + # coefficient in the within regression, and do not use a degree of freedom + rank_within = _lstsq(x, y, rcond=None)[2] if self.has_constant: w_sum = w.sum() y_gm = (w * self.dependent.values2d).sum(0) / w_sum @@ -3002,13 +3025,14 @@ def fit( wybar = self.dependent.mean("entity", weights=self.weights) wxbar = self.exog.mean("entity", weights=self.weights) - params = _lstsq(np.asarray(wxbar), np.asarray(wybar), rcond=None)[0] + params, _, rank_between, _ = _lstsq( + np.asarray(wxbar), np.asarray(wybar), rcond=None + ) wu = np.asarray(wybar) - np.asarray(wxbar) @ params nobs = weps.shape[0] neffects = wu.shape[0] - nvar = x.shape[1] - sigma2_e = float(np.squeeze(weps.T @ weps)) / (nobs - nvar - neffects + 1) + sigma2_e = float(np.squeeze(weps.T @ weps)) / (nobs - neffects - rank_within) ssr = float(np.squeeze(wu.T @ wu)) t = np.asarray(self.dependent.count("entity")) unbalanced = np.ptp(t) != 0 @@ -3019,11 +3043,17 @@ def fit( denom = means.T @ means sums = wx_df.groupby(level=0).sum().values num = sums.T @ sums - tr = np.trace(np.linalg.inv(denom) @ num) - sigma2_u = max(0, (ssr - (neffects - nvar) * sigma2_e) / (nobs - tr)) + # The means of the regressors can be collinear, e.g., the dummy variables + # of two periods that are observed for every entity have the same mean + # in all entities. The pseudo inverse uses the regressors that are not + # collinear + tr = np.trace(np.linalg.pinv(denom, hermitian=True) @ num) + sigma2_u = max( + 0, (ssr - (neffects - rank_between) * sigma2_e) / (nobs - tr) + ) else: t_bar = neffects / ((1.0 / t).sum()) - sigma2_u = max(0, ssr / (neffects - nvar) - sigma2_e / t_bar) + sigma2_u = max(0, ssr / (neffects - rank_between) - sigma2_e / t_bar) rho = sigma2_u / (sigma2_u + sigma2_e) theta = 1.0 - np.sqrt(sigma2_e / (t * sigma2_u + sigma2_e)) diff --git a/linearmodels/tests/panel/results/re-variance-components-reference.R b/linearmodels/tests/panel/results/re-variance-components-reference.R new file mode 100644 index 0000000000..12b267e6d4 --- /dev/null +++ b/linearmodels/tests/panel/results/re-variance-components-reference.R @@ -0,0 +1,91 @@ +# Reference values for the variance components of RandomEffects, from the R +# package plm, which are read by ../test_random_effects_variance.py. +# +# Run from this directory: Rscript re-variance-components-reference.R +# Requires: plm (tested with 2.6.7 on R 4.6.1) +# Input: ../../../datasets/wage_panel/wage_panel.csv.bz2 +# Output: re-variance-components-reference.csv +# +# The model is the Swamy-Arora random effects model of plm, +# plm(..., model = "random"), whose estimators of the variance components use the +# degrees of freedom of the regressions that estimate them: +# +# sigma2_idios = RSS(within) / (n - N - K_w) +# sigma2_id = RSS(between) / (N - K_b) - sigma2_idios / Tbar +# +# where K_w is the number of coefficients of the within regression, which are the +# regressors that vary within the entities and are not collinear after the +# entities' means are removed, and K_b is the number of coefficients of the +# between regression, including the constant. Time-invariant regressors (black, +# hisp, educ) are not in the within regression, and the dummy variables of the +# years are collinear with the constant in the between regression of a balanced +# panel. The degrees of freedom do not count them. +# +# Stata's xtreg, re uses the same degrees of freedom: for the Cornwell and Rupert +# wage panel (A. C. Cameron's log of psidextract.dta), xtreg lwage exp exp2 wks ed, +# re reports sigma_e = .15220316, which is also the sigma_e of xtreg lwage exp exp2 +# wks, fe, where ed is omitted because of collinearity (n - N - 3 degrees of +# freedom). The same model has sigma_e = 0.1522032 in plm. +# +# Cases, all of the wage panel: +# balanced-notebook the model of the panel_examples notebook: year dummies, +# time-invariant regressors, and exper, which is +# collinear with the year dummies in the within regression +# balanced-time-invariant time-invariant regressors without year dummies +# balanced-time-varying regressors that all vary within the entities + +suppressMessages(library(plm)) + +cat(R.version.string, "and plm", as.character(packageVersion("plm")), "\n") + +# CSV files with "\n" line endings on every platform +write_csv <- function(x, path, ...) { + con <- file(path, open = "wb") + on.exit(close(con)) + write.csv(x, con, row.names = FALSE, eol = "\n", ...) +} + +wages <- read.csv(bzfile("../../../datasets/wage_panel/wage_panel.csv.bz2")) + +formulas <- list( + notebook = lwage ~ black + hisp + exper + expersq + married + educ + union + + factor(year), + `time-invariant` = lwage ~ black + hisp + exper + expersq + married + educ + union, + `time-varying` = lwage ~ expersq + married + union +) +cases <- list( + list(name = "balanced-notebook", data = wages, formula = "notebook"), + list(name = "balanced-time-invariant", data = wages, formula = "time-invariant"), + list(name = "balanced-time-varying", data = wages, formula = "time-varying") +) + +rows <- list() +for (case in cases) { + pd <- pdata.frame(case$data, index = c("nr", "year")) + fit <- plm(formulas[[case$formula]], data = pd, model = "random") + components <- ercomp(fit) + theta <- as.numeric(components$theta) + est <- coef(fit) + se <- sqrt(diag(vcov(fit))) + within <- plm(formulas[[case$formula]], data = pd, model = "within") + between <- plm(formulas[[case$formula]], data = pd, model = "between") + rows[[case$name]] <- rbind( + data.frame(case = case$name, quantity = "nobs", value = nobs(fit)), + data.frame(case = case$name, quantity = "df_within", value = df.residual(within)), + data.frame(case = case$name, quantity = "df_between", value = df.residual(between)), + data.frame(case = case$name, quantity = "sigma2_idios", + value = unname(components$sigma2["idios"])), + data.frame(case = case$name, quantity = "sigma2_id", + value = unname(components$sigma2["id"])), + data.frame(case = case$name, quantity = "theta_mean", value = mean(theta)), + data.frame(case = case$name, quantity = paste0("coef", seq_along(est) - 1), + value = unname(est)), + data.frame(case = case$name, quantity = paste0("se", seq_along(est) - 1), + value = unname(se)) + ) + cat(sprintf("%-27s nobs %5d df within %5d df between %4d idios %.7f id %.7f\n", + case$name, nobs(fit), df.residual(within), df.residual(between), + components$sigma2["idios"], components$sigma2["id"])) +} +write_csv(do.call(rbind, rows), "re-variance-components-reference.csv") +cat("wrote re-variance-components-reference.csv\n") diff --git a/linearmodels/tests/panel/results/re-variance-components-reference.csv b/linearmodels/tests/panel/results/re-variance-components-reference.csv new file mode 100644 index 0000000000..7045a9a892 --- /dev/null +++ b/linearmodels/tests/panel/results/re-variance-components-reference.csv @@ -0,0 +1,73 @@ +"case","quantity","value" +"balanced-notebook","nobs",4360 +"balanced-notebook","df_within",3805 +"balanced-notebook","df_between",537 +"balanced-notebook","sigma2_idios",0.123193989960515 +"balanced-notebook","sigma2_id",0.105367200462246 +"balanced-notebook","theta_mean",0.642910879665865 +"balanced-notebook","coef0",0.0235864096542448 +"balanced-notebook","coef1",-0.139376734111937 +"balanced-notebook","coef2",0.0217317306549108 +"balanced-notebook","coef3",0.105754523245461 +"balanced-notebook","coef4",-0.00472394273755648 +"balanced-notebook","coef5",0.0639860359475824 +"balanced-notebook","coef6",0.091876272649444 +"balanced-notebook","coef7",0.106134422744568 +"balanced-notebook","coef8",0.0404619940706127 +"balanced-notebook","coef9",0.0309211619343808 +"balanced-notebook","coef10",0.0202806373503664 +"balanced-notebook","coef11",0.043118684873118 +"balanced-notebook","coef12",0.0578154064233405 +"balanced-notebook","coef13",0.0919475382668158 +"balanced-notebook","coef14",0.134928886286404 +"balanced-notebook","se0",0.150668257591444 +"balanced-notebook","se1",0.0477228164328931 +"balanced-notebook","se2",0.042606290027916 +"balanced-notebook","se3",0.0153668157460012 +"balanced-notebook","se4",0.00068949694503348 +"balanced-notebook","se5",0.0167742437355273 +"balanced-notebook","se6",0.0106597040958363 +"balanced-notebook","se7",0.0178538555251242 +"balanced-notebook","se8",0.024694610713068 +"balanced-notebook","se9",0.0323416128411864 +"balanced-notebook","se10",0.0415819882406935 +"balanced-notebook","se11",0.0513163474941151 +"balanced-notebook","se12",0.0612323120313417 +"balanced-notebook","se13",0.0712292614222332 +"balanced-notebook","se14",0.0813135284535609 +"balanced-time-invariant","nobs",4360 +"balanced-time-invariant","df_within",3811 +"balanced-time-invariant","df_between",537 +"balanced-time-invariant","sigma2_idios",0.123380320307742 +"balanced-time-invariant","sigma2_id",0.105343909168843 +"balanced-time-invariant","theta_mean",0.642640933868348 +"balanced-time-invariant","coef0",-0.107464203974044 +"balanced-time-invariant","coef1",-0.144130691111868 +"balanced-time-invariant","coef2",0.0201510730060844 +"balanced-time-invariant","coef3",0.112119493522384 +"balanced-time-invariant","coef4",-0.00406885475618792 +"balanced-time-invariant","coef5",0.0627951179727595 +"balanced-time-invariant","coef6",0.101224614698513 +"balanced-time-invariant","coef7",0.107378852592373 +"balanced-time-invariant","se0",0.110705725594299 +"balanced-time-invariant","se1",0.0476148274269489 +"balanced-time-invariant","se2",0.042601124174701 +"balanced-time-invariant","se3",0.00826087205583039 +"balanced-time-invariant","se4",0.000591825600027988 +"balanced-time-invariant","se5",0.0167728540561277 +"balanced-time-invariant","se6",0.00891328987417564 +"balanced-time-invariant","se7",0.0178300147703687 +"balanced-time-varying","nobs",4360 +"balanced-time-varying","df_within",3812 +"balanced-time-varying","df_between",541 +"balanced-time-varying","sigma2_idios",0.129581486663583 +"balanced-time-varying","sigma2_id",0.122493126638203 +"balanced-time-varying","theta_mean",0.658254629973685 +"balanced-time-varying","coef0",1.40520815095454 +"balanced-time-varying","coef1",0.00319171996193454 +"balanced-time-varying","coef2",0.131555169153797 +"balanced-time-varying","coef3",0.103450759968374 +"balanced-time-varying","se0",0.0192679695028832 +"balanced-time-varying","se1",0.000181336837300841 +"balanced-time-varying","se2",0.0168817492705327 +"balanced-time-varying","se3",0.0184859713861548 diff --git a/linearmodels/tests/panel/test_random_effects_variance.py b/linearmodels/tests/panel/test_random_effects_variance.py new file mode 100644 index 0000000000..b1bec04fd4 --- /dev/null +++ b/linearmodels/tests/panel/test_random_effects_variance.py @@ -0,0 +1,281 @@ +""" +The variance components of RandomEffects + +The idiosyncratic variance is the residual variance of the within regression, and +the variance of the effects is estimated from the residuals of the between +regression. The degrees of freedom of the variances are those of the regressions: +the within regression has the coefficients of the regressors that vary within +entities and are not collinear once the means of the entities are removed, and the +between regression has the coefficients of the regressors that are not collinear +across entities. A time-invariant regressor, such as education in a wage panel, is +not in the within regression, and neither is a regressor that is collinear with the +dummy variables of the time periods within entities, such as experience. A +dummy variable of the time periods can be collinear with the constant in the between +regression. The degrees of freedom do not count them. + +The references are + +* the residual variance of the regression on the dummy variables of the entities, + which is the idiosyncratic variance of Swamy and Arora; +* the estimators written out in NumPy using the ranks of the regressions; +* plm, in results/re-variance-components-reference.csv, for the wage panel of Vella + and Verbeek. The R script that produced it is + results/re-variance-components-reference.R. Stata's xtreg, re counts the + coefficients in the same way: sigma_e of a model with a time-invariant regressor is + that of xtreg, fe, which omits the regressor. +""" + +import os + +import numpy as np +from numpy.testing import assert_allclose +import pandas as pd +import pytest + +from linearmodels.datasets import wage_panel +from linearmodels.panel import PanelOLS, RandomEffects + +RESULTS = os.path.join(os.path.split(os.path.abspath(__file__))[0], "results") +NPERIOD = 6 +DUMMIES = [f"d{p}" for p in range(1, NPERIOD)] +DESIGNS = { + "time-varying": ["const", "x1"], + "no-constant": ["x1"], + "time-invariant": ["const", "x1", "z"], + "time-invariant-no-constant": ["x1", "z"], + # exper is collinear with the dummies once the means of the entities are removed + "collinear-with-dummies": ["const", "x1", "z", "exper", *DUMMIES], +} + + +def make_data(balance, seed=0, nentity=40): + """ + A panel with a time-invariant regressor (z) and one that is collinear with the + dummies of the time periods within entities (exper), as the experience in a + wage panel. The unbalanced panels drop observations from the last three periods + ("collinear-means") or from all periods + """ + rg = np.random.default_rng(seed) + index = pd.MultiIndex.from_product( + [np.arange(nentity), np.arange(NPERIOD)], names=["entity", "time"] + ) + entity = index.codes[0] + period = index.codes[1] + nobs = nentity * NPERIOD + effect = rg.standard_normal(nentity)[entity] + x1 = 0.5 * effect + rg.standard_normal(nobs) + z = rg.standard_normal(nentity)[entity] + exper = period + rg.integers(0, 5, nentity)[entity] + shock = (1 + 0.3 * np.abs(x1)) * rg.standard_normal(nobs) + data = pd.DataFrame( + { + "y": 1 + 0.8 * x1 + 0.5 * z + 0.05 * exper + 1.2 * effect + shock, + "const": 1.0, + "x1": x1, + "z": z, + "exper": exper, + "w": rg.chisquare(4, nobs) / 4 + 0.25, + }, + index=index, + ) + for p in range(1, NPERIOD): + data[f"d{p}"] = (period == p).astype(float) + if balance != "balanced": + keep = rg.random((nentity, NPERIOD)) > 0.3 + if balance == "unbalanced-collinear-means": + # The dummies of the first periods have the same mean in all entities + keep[:, :3] = True + for i in range(nentity): + keep[i, rg.choice(NPERIOD, 3, replace=False)] = True + data = data.loc[keep.ravel()] + return data + + +@pytest.fixture( + params=["balanced", "unbalanced", "unbalanced-collinear-means"], scope="module" +) +def data(request): + return make_data(request.param) + + +def lsdv_variance(data, columns, weights): + """ + The residual variance of the regression on the dummy variables of the entities + + The regressors of the model that are not identified by the dummy variables + (the time-invariant ones, which they span) do not reduce the degrees of freedom. + """ + dummies = pd.get_dummies(data.index.get_level_values(0)).astype(float) + regressors = [c for c in columns if c != "const"] + exog = np.column_stack([data[regressors].to_numpy(), dummies.to_numpy()]) + root_w = np.ones(data.shape[0]) + if weights is not None: + root_w = np.sqrt((weights / weights.mean()).to_numpy()) + wexog, wy = root_w[:, None] * exog, root_w * data["y"].to_numpy() + beta, _, rank, _ = np.linalg.lstsq(wexog, wy, rcond=None) + return float(np.sum((wy - wexog @ beta) ** 2)) / (data.shape[0] - rank) + + +def swar_reference(y, x, weights=None): + """ + The estimators of Swamy and Arora written out, with the ranks of the regressions + + Returns the idiosyncratic variance, the variance of the effects and the numbers + of coefficients of the within and the between regressions + """ + entity = y.index.get_level_values(0) + weights = pd.Series(1.0, index=y.index) if weights is None else weights + weights = weights / weights.mean() + total = weights.groupby(entity).transform("sum") + + def weighted_means(frame): + weighted = frame.mul(weights, axis=0) + return weighted.groupby(entity).transform("sum").div(total, axis=0) + + root_w = np.sqrt(weights.to_numpy()) + xd = root_w[:, None] * (x - weighted_means(x)).to_numpy() + yd = root_w * (y - weighted_means(y.to_frame()).iloc[:, 0]).to_numpy() + beta, _, rank_within, _ = np.linalg.lstsq(xd, yd, rcond=None) + nobs, nentity = x.shape[0], entity.nunique() + sigma2_eps = float(np.sum((yd - xd @ beta) ** 2)) / (nobs - nentity - rank_within) + # The regression of the means of the entities + xb = weighted_means(x).groupby(entity).first().to_numpy() + yb = weighted_means(y.to_frame()).iloc[:, 0].groupby(entity).first().to_numpy() + count = y.groupby(entity).size().to_numpy().astype(float) + rank_between = int(np.linalg.matrix_rank(xb)) + resid = yb - xb @ np.linalg.lstsq(xb, yb, rcond=None)[0] + tbar = nentity / np.sum(1 / count) + sigma2_effects = float(np.sum(resid**2)) / (nentity - rank_between) + sigma2_effects = max(0.0, sigma2_effects - sigma2_eps / tbar) + return sigma2_eps, sigma2_effects, rank_within, rank_between + + +@pytest.mark.parametrize("weighted", [False, True], ids=["unweighted", "weighted"]) +@pytest.mark.parametrize("design", list(DESIGNS)) +def test_idiosyncratic_variance_against_dummy_regression(data, design, weighted): + columns = DESIGNS[design] + weights = data["w"] if weighted else None + res = RandomEffects(data["y"], data[columns], weights=weights).fit() + expected = lsdv_variance(data, columns, weights) + assert_allclose(res.variance_decomposition["Residual"], expected) + + +@pytest.mark.parametrize("design", ["time-invariant", "time-invariant-no-constant"]) +def test_idiosyncratic_variance_is_that_of_the_fixed_effects(data, design): + # Stata: sigma_e of xtreg, re with a time-invariant regressor is the sigma_e + # of xtreg, fe, which omits the regressor + columns = DESIGNS[design] + re_res = RandomEffects(data["y"], data[columns]).fit() + fe_columns = [c for c in columns if c not in ("const", "z")] + fe_res = PanelOLS(data["y"], data[fe_columns], entity_effects=True).fit() + assert_allclose(re_res.variance_decomposition["Residual"], fe_res.s2) + + +@pytest.mark.parametrize("design", list(DESIGNS)) +def test_variance_components_against_estimators_written_out(data, design): + columns = DESIGNS[design] + res = RandomEffects(data["y"], data[columns]).fit() + sigma2_eps, sigma2_effects, *_ = swar_reference(data["y"], data[columns]) + decomposition = res.variance_decomposition + assert_allclose(decomposition["Residual"], sigma2_eps) + assert_allclose(decomposition["Effects"], sigma2_effects) + rho = sigma2_effects / (sigma2_effects + sigma2_eps) + assert_allclose(decomposition["Percent due to Effects"], rho) + theta = res.theta["theta"] + count = data.groupby(level=0).size().reindex(theta.index) + expected_theta = 1 - np.sqrt(sigma2_eps / (count * sigma2_effects + sigma2_eps)) + assert_allclose(theta, expected_theta) + + +@pytest.mark.parametrize("design", list(DESIGNS)) +def test_weighted_variance_components_against_estimators_written_out(data, design): + columns = DESIGNS[design] + res = RandomEffects(data["y"], data[columns], weights=data["w"]).fit() + sigma2_eps, sigma2_effects, *_ = swar_reference( + data["y"], data[columns], weights=data["w"] + ) + assert_allclose(res.variance_decomposition["Residual"], sigma2_eps) + assert_allclose(res.variance_decomposition["Effects"], sigma2_effects) + + +def test_the_designs_have_regressors_that_the_regressions_do_not_estimate(): + # The tests are only informative if the numbers of coefficients of the + # regressions differ from the number of regressors. The within regression + # has x1 and the dummies: z and the constant do not vary within entities, and + # exper is collinear with the dummies. In the means of a balanced panel the + # dummies are the constant 1/T. If periods 0, 1 and 2 are observed for every + # entity, the means of their dummies are equal (1/T_i), and so the constant is + # 3 times the mean of the dummy of period 1 plus the means of the other dummies + columns = DESIGNS["collinear-with-dummies"] + for balance, between in [ + ("balanced", 4), + ("unbalanced", len(columns)), + ("unbalanced-collinear-means", len(columns) - 2), + ]: + sample = make_data(balance) + ranks = swar_reference(sample["y"], sample[columns])[2:] + assert ranks == (1 + len(DUMMIES), between), balance + + +def wage_data(): + """The wage panel of the notebook""" + data = wage_panel.load() + year = pd.Categorical(data["year"]) + data = data.set_index(["nr", "year"]) + data["year"] = year + data["const"] = 1.0 + return data + + +TIME_INVARIANT = ["black", "hisp", "exper", "expersq", "married", "educ", "union"] +PLM_COLUMNS = { + "notebook": ["const", *TIME_INVARIANT, "year"], + "time-invariant": ["const", *TIME_INVARIANT], + "time-varying": ["const", "expersq", "married", "union"], +} +PLM_CASES = [ + ("balanced", "notebook"), + ("balanced", "time-invariant"), + ("balanced", "time-varying"), +] + + +@pytest.fixture(scope="module") +def plm_reference(): + path = os.path.join(RESULTS, "re-variance-components-reference.csv") + reference = pd.read_csv(path) + return { + case: group.set_index("quantity")["value"] + for case, group in reference.groupby("case") + } + + +@pytest.mark.parametrize(("balance", "model"), PLM_CASES) +def test_random_effects_against_plm(plm_reference, balance, model): + expected = plm_reference[f"{balance}-{model}"] + data = wage_data() + mod = RandomEffects(data["lwage"], data[PLM_COLUMNS[model]]) + res = mod.fit() + assert res.nobs == expected["nobs"] + decomposition = res.variance_decomposition + assert_allclose(decomposition["Residual"], expected["sigma2_idios"], rtol=1e-8) + assert_allclose(decomposition["Effects"], expected["sigma2_id"], rtol=1e-8) + # plm reports the theta of each observation + entities = mod.dependent.index.get_level_values(0) + theta = res.theta["theta"].reindex(entities) + assert_allclose(theta.mean(), expected["theta_mean"], rtol=1e-8) + nparam = res.params.shape[0] + coefficients = expected[[f"coef{i}" for i in range(nparam)]].to_numpy() + std_errors = expected[[f"se{i}" for i in range(nparam)]].to_numpy() + assert_allclose(res.params.to_numpy(), coefficients, rtol=1e-7, atol=1e-10) + assert_allclose(res.std_errors.to_numpy(), std_errors, rtol=1e-7) + # The estimators written out use the numbers of coefficients that plm does + sigma2_eps, sigma2_effects, rank_within, rank_between = swar_reference( + mod.dependent.dataframe["lwage"], + mod.exog.dataframe, + ) + nentity = entities.nunique() + assert expected["df_within"] == res.nobs - nentity - rank_within + assert expected["df_between"] == nentity - rank_between + assert_allclose(sigma2_eps, expected["sigma2_idios"], rtol=1e-8) + assert_allclose(sigma2_effects, expected["sigma2_id"], rtol=1e-8) From 6eda416e80c7fa4ce1295eafb83cbbebf0d865d4 Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Mon, 12 Oct 2026 15:10:56 +0100 Subject: [PATCH 2/4] BUG: Weight the between regression of the small_sample RandomEffects estimator With small_sample=True and an unbalanced panel, RandomEffects estimates the variance of the effects with the estimator of Baltagi and Chang (1994), whose sum of squares of the between regression is the sum of squares of the regression of the data replaced by the means of the entities, and weights each mean by the number of observations of the entity. The trace that is subtracted from the number of observations is that of the same regression. The residuals came from the unweighted regression of the means instead, so the sum of squares and the trace did not belong to the same estimator. A balanced panel is not affected, since all the weights are equal, and neither is the default estimator, which uses the harmonic mean of the numbers of observations as Stata does. The residuals are now those of the regression that weights the means by the numbers of observations. The variance of the effects equals that of plm (ercomp, whose default for an unbalanced panel is this estimator) to 8 digits. For the unbalanced wage panels of the tests it is 0.1083583 and 0.1076269 instead of 0.1083683 and 0.1076332. The difference grows with the differences in the numbers of observations: it was 5% for a panel in which entities enter and leave at different dates. --- linearmodels/panel/model.py | 15 ++++- .../re-variance-components-reference.R | 30 +++++++-- .../re-variance-components-reference.csv | 56 ++++++++++++++++ .../panel/test_random_effects_variance.py | 64 +++++++++++++++---- 4 files changed, 146 insertions(+), 19 deletions(-) diff --git a/linearmodels/panel/model.py b/linearmodels/panel/model.py index d2b2b277d9..206e90117f 100644 --- a/linearmodels/panel/model.py +++ b/linearmodels/panel/model.py @@ -3001,6 +3001,11 @@ def fit( (less one for the constant in the within regression), so that the variances, theta and the estimates were slightly different for models with time-invariant or collinear regressors, and for models without a constant. + With ``small_sample`` and an unbalanced panel, the residuals of the between + regression that the estimator of Baltagi and Chang (1994) uses are those of + the regression that weights the entities by their numbers of observations. + Earlier versions used the residuals of the unweighted regression of the + means, which gave a different variance of the effects. """ w = self.weights.values2d root_w = cast("linearmodels.typing.data.Float64Array", np.sqrt(w)) @@ -3037,7 +3042,15 @@ def fit( t = np.asarray(self.dependent.count("entity")) unbalanced = np.ptp(t) != 0 if small_sample and unbalanced: - ssr = float(np.squeeze((t * wu).T @ wu)) + # The between regression of the estimator of Baltagi and Chang (1994) is + # the regression of the data replaced by the means of the entities. It + # weights the means by the numbers of observations, as the trace below + # does, and the residuals of the regression of the means do not + root_t = np.sqrt(t) + xbar, ybar = np.asarray(wxbar), np.asarray(wybar) + params_t = _lstsq(root_t * xbar, root_t * ybar, rcond=None)[0] + wu_t = ybar - xbar @ params_t + ssr = float(np.squeeze((t * wu_t).T @ wu_t)) wx_df = cast("DataFrame", root_w * self.exog.dataframe) means = wx_df.groupby(level=0).transform("mean").values denom = means.T @ means diff --git a/linearmodels/tests/panel/results/re-variance-components-reference.R b/linearmodels/tests/panel/results/re-variance-components-reference.R index 12b267e6d4..644ff1dfb3 100644 --- a/linearmodels/tests/panel/results/re-variance-components-reference.R +++ b/linearmodels/tests/panel/results/re-variance-components-reference.R @@ -7,11 +7,11 @@ # Output: re-variance-components-reference.csv # # The model is the Swamy-Arora random effects model of plm, -# plm(..., model = "random"), whose estimators of the variance components use the +# plm(..., model = "random"), whose estimator of the variance components uses the # degrees of freedom of the regressions that estimate them: # # sigma2_idios = RSS(within) / (n - N - K_w) -# sigma2_id = RSS(between) / (N - K_b) - sigma2_idios / Tbar +# sigma2_id = RSS(between) / (N - K_b) - sigma2_idios / Tbar (balanced) # # where K_w is the number of coefficients of the within regression, which are the # regressors that vary within the entities and are not collinear after the @@ -19,7 +19,13 @@ # between regression, including the constant. Time-invariant regressors (black, # hisp, educ) are not in the within regression, and the dummy variables of the # years are collinear with the constant in the between regression of a balanced -# panel. The degrees of freedom do not count them. +# panel. The degrees of freedom do not count them. For an unbalanced panel plm +# uses the estimator of Baltagi and Chang (1994) (dfcor = 3), which is the +# small_sample option of RandomEffects.fit. That branch of ercomp counts the +# columns of the model matrix of the within regression for K_w, so that it counts +# a regressor that is aliased (exper with the year dummies), and this is not the +# df.residual of the within model of plm. The unbalanced cases have no aliased +# regressor. # # Stata's xtreg, re uses the same degrees of freedom: for the Cornwell and Rupert # wage panel (A. C. Cameron's log of psidextract.dta), xtreg lwage exp exp2 wks ed, @@ -27,12 +33,14 @@ # wks, fe, where ed is omitted because of collinearity (n - N - 3 degrees of # freedom). The same model has sigma_e = 0.1522032 in plm. # -# Cases, all of the wage panel: +# Cases (the unbalanced cases drop observations by a fixed rule): # balanced-notebook the model of the panel_examples notebook: year dummies, # time-invariant regressors, and exper, which is # collinear with the year dummies in the within regression # balanced-time-invariant time-invariant regressors without year dummies # balanced-time-varying regressors that all vary within the entities +# unbalanced-year-dummies year dummies and time-invariant regressors, without exper +# unbalanced-time-invariant as balanced-time-invariant suppressMessages(library(plm)) @@ -47,16 +55,28 @@ write_csv <- function(x, path, ...) { wages <- read.csv(bzfile("../../../datasets/wage_panel/wage_panel.csv.bz2")) +# The observations that the unbalanced cases drop: 588 of 4,360, so that each +# entity has between 6 and 8 observations and the means of the year dummies differ +# between the entities +drop <- ((wages$nr * 7 + wages$year * 13) %% 5 == 0) & ((wages$nr + wages$year) %% 3 != 0) +cat("unbalanced cases drop", sum(drop), "of", nrow(wages), "observations\n") + formulas <- list( notebook = lwage ~ black + hisp + exper + expersq + married + educ + union + factor(year), + `year-dummies` = lwage ~ black + hisp + expersq + married + educ + union + + factor(year), `time-invariant` = lwage ~ black + hisp + exper + expersq + married + educ + union, `time-varying` = lwage ~ expersq + married + union ) cases <- list( list(name = "balanced-notebook", data = wages, formula = "notebook"), list(name = "balanced-time-invariant", data = wages, formula = "time-invariant"), - list(name = "balanced-time-varying", data = wages, formula = "time-varying") + list(name = "balanced-time-varying", data = wages, formula = "time-varying"), + list(name = "unbalanced-year-dummies", data = wages[!drop, ], + formula = "year-dummies"), + list(name = "unbalanced-time-invariant", data = wages[!drop, ], + formula = "time-invariant") ) rows <- list() diff --git a/linearmodels/tests/panel/results/re-variance-components-reference.csv b/linearmodels/tests/panel/results/re-variance-components-reference.csv index 7045a9a892..92c414c316 100644 --- a/linearmodels/tests/panel/results/re-variance-components-reference.csv +++ b/linearmodels/tests/panel/results/re-variance-components-reference.csv @@ -71,3 +71,59 @@ "balanced-time-varying","se1",0.000181336837300841 "balanced-time-varying","se2",0.0168817492705327 "balanced-time-varying","se3",0.0184859713861548 +"unbalanced-year-dummies","nobs",3772 +"unbalanced-year-dummies","df_within",3217 +"unbalanced-year-dummies","df_between",531 +"unbalanced-year-dummies","sigma2_idios",0.128302344184146 +"unbalanced-year-dummies","sigma2_id",0.108358317305489 +"unbalanced-year-dummies","theta_mean",0.618184470344336 +"unbalanced-year-dummies","coef0",0.605548790479954 +"unbalanced-year-dummies","coef1",-0.125821461978439 +"unbalanced-year-dummies","coef2",0.0174929636180146 +"unbalanced-year-dummies","coef3",-0.00130071500004915 +"unbalanced-year-dummies","coef4",0.0837165981326888 +"unbalanced-year-dummies","coef5",0.0646943293668789 +"unbalanced-year-dummies","coef6",0.127259889974332 +"unbalanced-year-dummies","coef7",0.112573187254709 +"unbalanced-year-dummies","coef8",0.188019604193127 +"unbalanced-year-dummies","coef9",0.242700240753059 +"unbalanced-year-dummies","coef10",0.333987251802706 +"unbalanced-year-dummies","coef11",0.394081555168255 +"unbalanced-year-dummies","coef12",0.470921208312212 +"unbalanced-year-dummies","coef13",0.562894669724611 +"unbalanced-year-dummies","se0",0.127130321119421 +"unbalanced-year-dummies","se1",0.0491356466494286 +"unbalanced-year-dummies","se2",0.0439467291801474 +"unbalanced-year-dummies","se3",0.00053586489757734 +"unbalanced-year-dummies","se4",0.0182603878294297 +"unbalanced-year-dummies","se5",0.0102145931419279 +"unbalanced-year-dummies","se6",0.0195234302249132 +"unbalanced-year-dummies","se7",0.0242011174360354 +"unbalanced-year-dummies","se8",0.0252961420112048 +"unbalanced-year-dummies","se9",0.0282487685579614 +"unbalanced-year-dummies","se10",0.0324318401001503 +"unbalanced-year-dummies","se11",0.0383367319727384 +"unbalanced-year-dummies","se12",0.0457883832019504 +"unbalanced-year-dummies","se13",0.0547496534149072 +"unbalanced-time-invariant","nobs",3772 +"unbalanced-time-invariant","df_within",3223 +"unbalanced-time-invariant","df_between",537 +"unbalanced-time-invariant","sigma2_idios",0.128427843646611 +"unbalanced-time-invariant","sigma2_id",0.107626857054882 +"unbalanced-time-invariant","theta_mean",0.6169198975766 +"unbalanced-time-invariant","coef0",-0.125812807588583 +"unbalanced-time-invariant","coef1",-0.151894553284766 +"unbalanced-time-invariant","coef2",0.0157023105776524 +"unbalanced-time-invariant","coef3",0.115990971886776 +"unbalanced-time-invariant","coef4",-0.00431687671630469 +"unbalanced-time-invariant","coef5",0.0683264347224204 +"unbalanced-time-invariant","coef6",0.100948952871547 +"unbalanced-time-invariant","coef7",0.124414893965101 +"unbalanced-time-invariant","se0",0.113929832936069 +"unbalanced-time-invariant","se1",0.0487220333669172 +"unbalanced-time-invariant","se2",0.0435788643259302 +"unbalanced-time-invariant","se3",0.00910200119006807 +"unbalanced-time-invariant","se4",0.000653075903819231 +"unbalanced-time-invariant","se5",0.0182689960010359 +"unbalanced-time-invariant","se6",0.00913644762476636 +"unbalanced-time-invariant","se7",0.0193869180875927 diff --git a/linearmodels/tests/panel/test_random_effects_variance.py b/linearmodels/tests/panel/test_random_effects_variance.py index b1bec04fd4..b8c2707a0e 100644 --- a/linearmodels/tests/panel/test_random_effects_variance.py +++ b/linearmodels/tests/panel/test_random_effects_variance.py @@ -31,6 +31,7 @@ from numpy.testing import assert_allclose import pandas as pd import pytest +from scipy.linalg import qr from linearmodels.datasets import wage_panel from linearmodels.panel import PanelOLS, RandomEffects @@ -116,7 +117,7 @@ def lsdv_variance(data, columns, weights): return float(np.sum((wy - wexog @ beta) ** 2)) / (data.shape[0] - rank) -def swar_reference(y, x, weights=None): +def swar_reference(y, x, weights=None, small_sample=False): """ The estimators of Swamy and Arora written out, with the ranks of the regressions @@ -143,11 +144,26 @@ def weighted_means(frame): yb = weighted_means(y.to_frame()).iloc[:, 0].groupby(entity).first().to_numpy() count = y.groupby(entity).size().to_numpy().astype(float) rank_between = int(np.linalg.matrix_rank(xb)) - resid = yb - xb @ np.linalg.lstsq(xb, yb, rcond=None)[0] - tbar = nentity / np.sum(1 / count) - sigma2_effects = float(np.sum(resid**2)) / (nentity - rank_between) - sigma2_effects = max(0.0, sigma2_effects - sigma2_eps / tbar) - return sigma2_eps, sigma2_effects, rank_within, rank_between + if not small_sample: + resid = yb - xb @ np.linalg.lstsq(xb, yb, rcond=None)[0] + tbar = nentity / np.sum(1 / count) + sigma2_effects = float(np.sum(resid**2)) / (nentity - rank_between) + sigma2_effects = max(0.0, sigma2_effects - sigma2_eps / tbar) + return sigma2_eps, sigma2_effects, rank_within, rank_between + # Baltagi and Chang (1994) for unbalanced panels: the means are weighted by the + # number of observations. This is the regression of the data replaced by the + # means of the entities. The trace uses columns that are not collinear + root_t = np.sqrt(count) + coef = np.linalg.lstsq(root_t[:, None] * xb, root_t * yb, rcond=None)[0] + rss = float(np.sum(count * (yb - xb @ coef) ** 2)) + basis = qr(xb, pivoting=True, mode="economic")[2][:rank_between] + xr = xb[:, basis] + trace = np.trace( + np.linalg.inv((count[:, None] * xr).T @ xr) + @ ((count[:, None] * xr).T @ (count[:, None] * xr)) + ) + sigma2_effects = (rss - (nentity - rank_between) * sigma2_eps) / (nobs - trace) + return sigma2_eps, max(0.0, float(sigma2_effects)), rank_within, rank_between @pytest.mark.parametrize("weighted", [False, True], ids=["unweighted", "weighted"]) @@ -171,11 +187,14 @@ def test_idiosyncratic_variance_is_that_of_the_fixed_effects(data, design): assert_allclose(re_res.variance_decomposition["Residual"], fe_res.s2) +@pytest.mark.parametrize("small_sample", [False, True]) @pytest.mark.parametrize("design", list(DESIGNS)) -def test_variance_components_against_estimators_written_out(data, design): +def test_variance_components_against_estimators_written_out(data, design, small_sample): columns = DESIGNS[design] - res = RandomEffects(data["y"], data[columns]).fit() - sigma2_eps, sigma2_effects, *_ = swar_reference(data["y"], data[columns]) + res = RandomEffects(data["y"], data[columns]).fit(small_sample=small_sample) + sigma2_eps, sigma2_effects, *_ = swar_reference( + data["y"], data[columns], small_sample=small_sample + ) decomposition = res.variance_decomposition assert_allclose(decomposition["Residual"], sigma2_eps) assert_allclose(decomposition["Effects"], sigma2_effects) @@ -217,9 +236,13 @@ def test_the_designs_have_regressors_that_the_regressions_do_not_estimate(): assert ranks == (1 + len(DUMMIES), between), balance -def wage_data(): - """The wage panel of the notebook""" +def wage_data(unbalanced): + """The wage panel of the notebook, and its unbalanced version in the R script""" data = wage_panel.load() + if unbalanced: + year, nr = data["year"], data["nr"] + drop = ((nr * 7 + year * 13) % 5 == 0) & ((nr + year) % 3 != 0) + data = data.loc[~drop] year = pd.Categorical(data["year"]) data = data.set_index(["nr", "year"]) data["year"] = year @@ -230,6 +253,16 @@ def wage_data(): TIME_INVARIANT = ["black", "hisp", "exper", "expersq", "married", "educ", "union"] PLM_COLUMNS = { "notebook": ["const", *TIME_INVARIANT, "year"], + "year-dummies": [ + "const", + "black", + "hisp", + "expersq", + "married", + "educ", + "union", + "year", + ], "time-invariant": ["const", *TIME_INVARIANT], "time-varying": ["const", "expersq", "married", "union"], } @@ -237,6 +270,8 @@ def wage_data(): ("balanced", "notebook"), ("balanced", "time-invariant"), ("balanced", "time-varying"), + ("unbalanced", "year-dummies"), + ("unbalanced", "time-invariant"), ] @@ -253,9 +288,11 @@ def plm_reference(): @pytest.mark.parametrize(("balance", "model"), PLM_CASES) def test_random_effects_against_plm(plm_reference, balance, model): expected = plm_reference[f"{balance}-{model}"] - data = wage_data() + unbalanced = balance == "unbalanced" + data = wage_data(unbalanced) mod = RandomEffects(data["lwage"], data[PLM_COLUMNS[model]]) - res = mod.fit() + # plm uses the estimator of Baltagi and Chang (1994) for an unbalanced panel + res = mod.fit(small_sample=unbalanced) assert res.nobs == expected["nobs"] decomposition = res.variance_decomposition assert_allclose(decomposition["Residual"], expected["sigma2_idios"], rtol=1e-8) @@ -273,6 +310,7 @@ def test_random_effects_against_plm(plm_reference, balance, model): sigma2_eps, sigma2_effects, rank_within, rank_between = swar_reference( mod.dependent.dataframe["lwage"], mod.exog.dataframe, + small_sample=unbalanced, ) nentity = entities.nunique() assert expected["df_within"] == res.nobs - nentity - rank_within From 65a8096072f9806263dac54ecb30ac2864617028 Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Mon, 12 Oct 2026 20:00:01 +0100 Subject: [PATCH 3/4] BUG: Raise a clear error when a RandomEffects variance has no degrees of freedom The variance of the idiosyncratic error divides the residual sum of squares of the within regression by nobs - N - rank, and the variance of the effects divides that of the between regression by N - rank, where rank is the rank of the regressors once the means of the entities are removed, and of the means of the regressors. A panel with one observation of each entity, or with as many entities as the means have rank, has no degrees of freedom left. The model raised a ZeroDivisionError, SciPy's "array must not contain infs or NaNs", or returned an infinite variance of the effects (and a theta of one, the within estimator). A panel with a single entity returned a variance of the effects of zero, which is not an estimate. Both degrees of freedom are now checked before they are used. A ValueError says which variance cannot be estimated, the numbers of observations or entities and the ranks, and how many observations or entities are needed. The tests compare the errors with the degrees of freedom that NumPy finds from the ranks of the demeaned regressors and of the means for a grid of balanced and unbalanced panels with and without a constant, and the idiosyncratic variance of the models that can be estimated with the residual variance of the regression on the dummy variables of the entities. --- linearmodels/panel/model.py | 35 +++- .../tests/panel/test_random_effects_no_df.py | 150 ++++++++++++++++++ 2 files changed, 178 insertions(+), 7 deletions(-) create mode 100644 linearmodels/tests/panel/test_random_effects_no_df.py diff --git a/linearmodels/panel/model.py b/linearmodels/panel/model.py index ff01ff9255..9e9d3544af 100644 --- a/linearmodels/panel/model.py +++ b/linearmodels/panel/model.py @@ -3080,7 +3080,9 @@ def fit( means of the entities are removed (experience, if the model has year dummies) is not in the within regression, and a regressor that is collinear with the others in the means of the entities (the year dummies and the - constant in a balanced panel) is not in the between regression. + constant in a balanced panel) is not in the between regression. A + ``ValueError`` is raised if either regression has no residual degrees of + freedom, since the variances cannot be estimated. .. versionchanged:: 7.1 @@ -3093,7 +3095,10 @@ def fit( regression that the estimator of Baltagi and Chang (1994) uses are those of the regression that weights the entities by their numbers of observations. Earlier versions used the residuals of the unweighted regression of the - means, which gave a different variance of the effects. + means, which gave a different variance of the effects. Models that do not + have residual degrees of freedom in the regressions raise a + ``ValueError``. Earlier versions raised an error from NumPy or SciPy, + or returned variances that were not defined. """ w = self.weights.values2d root_w = cast("linearmodels.typing.data.Float64Array", np.sqrt(w)) @@ -3125,7 +3130,25 @@ def fit( nobs = weps.shape[0] neffects = wu.shape[0] - sigma2_e = float(np.squeeze(weps.T @ weps)) / (nobs - neffects - rank_within) + df_within = nobs - neffects - rank_within + df_between = neffects - rank_between + if df_within <= 0: + raise ValueError( + "The variance of the idiosyncratic error cannot be estimated since " + "the within regression has no residual degrees of freedom. The model " + f"has {nobs} observations of {neffects} entities, and the regressors " + f"that vary within entities have a rank of {rank_within}, so that " + f"at least {neffects + rank_within + 1} observations are needed." + ) + if df_between <= 0: + raise ValueError( + "The variance of the effects cannot be estimated since the between " + "regression has no residual degrees of freedom. The model has " + f"{neffects} entities, and the means of the regressors have a rank " + f"of {rank_between}, so that at least {rank_between + 1} entities " + "are needed." + ) + sigma2_e = float(np.squeeze(weps.T @ weps)) / df_within ssr = float(np.squeeze(wu.T @ wu)) t = np.asarray(self.dependent.count("entity")) unbalanced = np.ptp(t) != 0 @@ -3149,12 +3172,10 @@ def fit( # in all entities. The pseudo inverse uses the regressors that are not # collinear tr = np.trace(np.linalg.pinv(denom, hermitian=True) @ num) - sigma2_u = max( - 0, (ssr - (neffects - rank_between) * sigma2_e) / (nobs - tr) - ) + sigma2_u = max(0, (ssr - df_between * sigma2_e) / (nobs - tr)) else: t_bar = neffects / ((1.0 / t).sum()) - sigma2_u = max(0, ssr / (neffects - rank_between) - sigma2_e / t_bar) + sigma2_u = max(0, ssr / df_between - sigma2_e / t_bar) rho = sigma2_u / (sigma2_u + sigma2_e) theta = 1.0 - np.sqrt(sigma2_e / (t * sigma2_u + sigma2_e)) diff --git a/linearmodels/tests/panel/test_random_effects_no_df.py b/linearmodels/tests/panel/test_random_effects_no_df.py new file mode 100644 index 0000000000..3bcb10ad0b --- /dev/null +++ b/linearmodels/tests/panel/test_random_effects_no_df.py @@ -0,0 +1,150 @@ +""" +RandomEffects where the regressions of the variance components have no degrees of freedom + +The variance of the idiosyncratic error is the residual sum of squares of the within +regression divided by its degrees of freedom, nobs - N - rank, where rank is the rank +of the regressors once the means of the entities are removed. The variance of the +effects is estimated from the between regression, which has N - rank observations of +the means of the regressors left. The variances are not defined if there are no +degrees of freedom, e.g., if there is a single observation of each entity or fewer +entities than regressors. These models raised an error from the division, from NumPy or +from SciPy, or returned an infinite variance of the effects. +""" + +import itertools + +import numpy as np +from numpy.testing import assert_allclose +import pandas as pd +import pytest + +from linearmodels.panel import RandomEffects + + +def make_panel(counts, nvar, const, seed=0): + rg = np.random.default_rng(seed) + entities = np.repeat(np.arange(len(counts)), counts) + times = np.concatenate([np.arange(count) for count in counts]) + index = pd.MultiIndex.from_arrays([entities, times], names=["entity", "time"]) + x = pd.DataFrame( + rg.standard_normal((entities.shape[0], nvar)), + index=index, + columns=[f"x{i}" for i in range(nvar)], + ) + if const: + x.insert(0, "const", 1.0) + y = pd.Series(rg.standard_normal(entities.shape[0]), index=index, name="y") + return y, x + + +def degrees_of_freedom(x): + """The residual degrees of freedom of the within and the between regressions""" + nobs = x.shape[0] + means = x.groupby(level=0).mean() + demeaned = x - x.groupby(level=0).transform("mean") + within = nobs - means.shape[0] - np.linalg.matrix_rank(demeaned.to_numpy()) + between = means.shape[0] - np.linalg.matrix_rank(means.to_numpy()) + return within, between + + +def lsdv_variance(y, x): + """Residual variance of the regression on the dummy variables of the entities""" + codes = pd.factorize(x.index.get_level_values(0))[0] + dummies = np.eye(codes.max() + 1)[codes] + design = np.column_stack([x.to_numpy(), dummies]) + resid = y.to_numpy() - design @ np.linalg.lstsq(design, y.to_numpy(), rcond=None)[0] + return resid @ resid / (design.shape[0] - np.linalg.matrix_rank(design)) + + +WITHIN = "variance of the idiosyncratic error cannot be estimated" +BETWEEN = "variance of the effects cannot be estimated" +BALANCED = [[t] * n for n in (1, 2, 3, 4) for t in (1, 2, 3)] +UNBALANCED = [[1, 2, 1], [1, 1, 2], [2, 2, 1], [1, 3, 2, 1], [3, 1, 4, 2, 2]] + + +def full_rank_designs(): + """The designs where the regressors are not collinear, as the model requires""" + for counts, nvar, const in itertools.product( + BALANCED + UNBALANCED, [1, 2, 3], [False, True] + ): + x = make_panel(counts, nvar, const)[1] + if np.linalg.matrix_rank(x.to_numpy()) == x.shape[1]: + yield counts, nvar, const + + +DESIGNS = list(full_rank_designs()) +DESIGN_IDS = [f"{c}-{k}-{'const' if const else 'noconst'}" for c, k, const in DESIGNS] + + +@pytest.mark.parametrize("small_sample", [False, True]) +@pytest.mark.parametrize(("counts", "nvar", "const"), DESIGNS, ids=DESIGN_IDS) +def test_variances_need_degrees_of_freedom(counts, nvar, const, small_sample): + y, x = make_panel(counts, nvar, const) + within, between = degrees_of_freedom(x) + mod = RandomEffects(y, x) + if within <= 0: + with pytest.raises(ValueError, match=WITHIN): + mod.fit(small_sample=small_sample) + elif between <= 0: + with pytest.raises(ValueError, match=BETWEEN): + mod.fit(small_sample=small_sample) + else: + res = mod.fit(small_sample=small_sample) + variances = res.variance_decomposition + assert np.all(np.isfinite(variances)) + assert variances["Residual"] > 0 + assert variances["Effects"] >= 0 + assert np.all(np.isfinite(res.params)) + assert_allclose(variances["Residual"], lsdv_variance(y, x)) + + +def test_each_cell_of_the_grid_is_used(): + # Some designs of the grid fail for each reason, and some are estimated + reasons = [] + for counts, nvar, const in DESIGNS: + within, between = degrees_of_freedom(make_panel(counts, nvar, const)[1]) + reasons.append("within" if within <= 0 else "between" if between <= 0 else "ok") + assert {"within", "between", "ok"} == set(reasons) + assert all(reasons.count(reason) >= 10 for reason in set(reasons)) + + +@pytest.mark.parametrize("small_sample", [False, True]) +def test_one_observation_of_each_entity(small_sample): + y, x = make_panel([1] * 10, 1, True) + # The 10 means leave nothing of the 10 observations, and the regressor does not vary + with pytest.raises(ValueError, match=WITHIN) as exc_info: + RandomEffects(y, x).fit(small_sample=small_sample) + assert "at least 11 observations are needed" in str(exc_info.value) + assert "10 observations of 10 entities" in str(exc_info.value) + + +@pytest.mark.parametrize("small_sample", [False, True]) +def test_as_many_regressors_as_the_observations_of_the_entities_leave(small_sample): + # 3 entities with 2 observations: 3 regressors use the 3 degrees of freedom + y, x = make_panel([2, 2, 2], 3, False) + with pytest.raises(ValueError, match=WITHIN) as exc_info: + RandomEffects(y, x).fit(small_sample=small_sample) + assert "rank of 3" in str(exc_info.value) + assert "at least 7 observations are needed" in str(exc_info.value) + + +@pytest.mark.parametrize("small_sample", [False, True]) +@pytest.mark.parametrize("counts", [[10], [3, 3], [4, 4]]) +def test_not_more_entities_than_the_rank_of_the_means(counts, small_sample): + # The constant and a regressor have a rank of 2 in the means of 1 or 2 entities + y, x = make_panel(counts, 1, True) + with pytest.raises(ValueError, match=BETWEEN) as exc_info: + RandomEffects(y, x).fit(small_sample=small_sample) + rank = min(len(counts), 2) + assert f"{len(counts)} entities" in str(exc_info.value) + assert f"at least {rank + 1} entities are needed" in str(exc_info.value) + + +@pytest.mark.parametrize("small_sample", [False, True]) +def test_one_more_entity_than_the_rank_of_the_means(small_sample): + # The smallest panel that can be estimated: 3 entities, a constant and a regressor + y, x = make_panel([3, 3, 3], 1, True) + assert degrees_of_freedom(x) == (9 - 3 - 1, 1) + res = RandomEffects(y, x).fit(small_sample=small_sample) + assert np.all(np.isfinite(res.variance_decomposition)) + assert_allclose(res.variance_decomposition["Residual"], lsdv_variance(y, x)) From c77c62016151d8c1da4f6ce400995459cdc2da22 Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Mon, 12 Oct 2026 20:34:29 +0100 Subject: [PATCH 4/4] BUG: Find the rank of the within regression of RandomEffects in the units of the regressors The degrees of freedom of the idiosyncratic variance use the rank of the regressors once the means of the entities are removed. A regressor that does not vary within entities is not left at exactly zero when the means are removed, because the mean of the same number repeated T times is not always that number in floating point arithmetic. It is left with rounding error, about 1e-16 times its size, and that is not small next to the other regressors if its units are large. The rank counted such a regressor, for example income in dollars next to the other regressors, which made the degrees of freedom one too small for each of them, the same mistake that the degrees of freedom of earlier versions made for every time-invariant regressor, and the variance depended on the units: for the wage panel with the years of education multiplied by 3.14e6 the idiosyncratic variance was 0.133625 instead of 0.133584. The rank is now that of the demeaned regressors divided by the norms of the regressors before the means were removed, so that the rounding error is 1e-16 whatever the units. Regressors that are not time-invariant and the constant do not change. The tests multiply a time-invariant regressor by 3.14e3, 3.14e6 and 3.14e8 for the balanced and unbalanced panels, and compare the idiosyncratic variance with that of the regression on the dummy variables of the entities, which does not depend on the units. --- linearmodels/panel/model.py | 9 +++++++-- .../tests/panel/test_random_effects_variance.py | 14 ++++++++++++++ 2 files changed, 21 insertions(+), 2 deletions(-) diff --git a/linearmodels/panel/model.py b/linearmodels/panel/model.py index 9e9d3544af..a377898bc8 100644 --- a/linearmodels/panel/model.py +++ b/linearmodels/panel/model.py @@ -3110,8 +3110,13 @@ def fit( x = demeaned_exog.values2d # The regressors that do not vary within entities (a constant, a time-invariant # regressor) or that are collinear once the means are removed have no - # coefficient in the within regression, and do not use a degree of freedom - rank_within = _lstsq(x, y, rcond=None)[2] + # coefficient in the within regression, and do not use a degree of freedom. + # Removing the means leaves rounding error in a regressor that does not vary, + # which is large if its units are large, so the rank is that of the regressors + # in the units that they had before the means were removed + x_norm = np.sqrt(((root_w * self.exog.values2d) ** 2).sum(0)) + x_norm[x_norm == 0.0] = 1.0 + rank_within = _lstsq(x / x_norm, y, rcond=None)[2] if self.has_constant: w_sum = w.sum() y_gm = (w * self.dependent.values2d).sum(0) / w_sum diff --git a/linearmodels/tests/panel/test_random_effects_variance.py b/linearmodels/tests/panel/test_random_effects_variance.py index b8c2707a0e..95b4dc40d4 100644 --- a/linearmodels/tests/panel/test_random_effects_variance.py +++ b/linearmodels/tests/panel/test_random_effects_variance.py @@ -187,6 +187,20 @@ def test_idiosyncratic_variance_is_that_of_the_fixed_effects(data, design): assert_allclose(re_res.variance_decomposition["Residual"], fe_res.s2) +@pytest.mark.parametrize("scale", [1e3, 1e6, 1e8]) +def test_variance_does_not_depend_on_the_units_of_a_regressor_that_does_not_vary( + data, scale +): + # Removing the means of a regressor that does not vary within entities leaves + # rounding error, which is large next to the other regressors if its units are + # large, and it is not a coefficient of the within regression + columns = DESIGNS["time-invariant"] + scaled = data.assign(z=data["z"] * np.pi * scale) + res = RandomEffects(scaled["y"], scaled[columns]).fit() + expected = lsdv_variance(scaled, columns, None) + assert_allclose(res.variance_decomposition["Residual"], expected) + + @pytest.mark.parametrize("small_sample", [False, True]) @pytest.mark.parametrize("design", list(DESIGNS)) def test_variance_components_against_estimators_written_out(data, design, small_sample):