BUG: Fix the degrees of freedom of the RandomEffects variance components

linearmodels/panel/model.py

@@ -3070,6 +3070,35 @@ def fit(
         ``debiased=True`` and ``group_debias=True`` give the clustered
         covariance of Stata's ``xtreg, re vce(cluster ...)``.
         See :ref:`panel-small-sample`.
+
+        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. A
+        ``ValueError`` is raised if either regression has no residual degrees of
+        freedom, since the variances cannot be estimated.
+
+        .. 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.
+           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. 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))
@@ -3079,6 +3108,15 @@ 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.
+        # 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
@@ -3090,28 +3128,59 @@ 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)
+        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
         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
             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 - df_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 / 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))

linearmodels/tests/panel/results/re-variance-components-reference.R

@@ -0,0 +1,111 @@
+# 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 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     (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
+# 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. 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,
+# 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 (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))
+
+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"))
+
+# 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 = "unbalanced-year-dummies", data = wages[!drop, ],
+       formula = "year-dummies"),
+  list(name = "unbalanced-time-invariant", data = wages[!drop, ],
+       formula = "time-invariant")
+)
+
+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")

linearmodels/tests/panel/results/re-variance-components-reference.csv

@@ -0,0 +1,129 @@
+"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
+"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