From 546e8b3119dd0cb7f1f52e66cc7a4534479dd1b6 Mon Sep 17 00:00:00 2001
From: Matthieu Gomez <gomez.matthieu@gmail.com>
Date: Fri, 30 Oct 2026 09:57:49 -0400
Subject: [PATCH 1/4] Fix absorbed degrees of freedom and first-stage F

- Do not count fixed-effect levels that are redundant across fixed effects
- Count fixed effects nested in a cluster only once, for the constant
- Use the first-stage degrees of freedom in F_kp when over-identified
- Require Vcov 0.8.4 and bump version to 2.2.0
---
 CHANGELOG.md              |  7 +++++
 Project.toml              |  4 +--
 src/FixedEffectModel.jl   |  3 ++-
 src/fit.jl                | 25 +++++++-----------
 src/partial_out.jl        |  2 +-
 src/utils/fixedeffects.jl | 55 +++++++++++++++++++++++++++++++++++++++
 test/fit.jl               | 53 +++++++++++++++++++++++++++----------
 7 files changed, 115 insertions(+), 34 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index fdcf8482..7dbe8397 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,12 @@
 # Changelog
 
+## 2.2.0
+
+- Fix the degrees of freedom absorbed by fixed effects: levels that are redundant across fixed effects no longer count (one per connected component), and fixed effects nested in a cluster count only once, for the constant. With two or more fixed effects, this changes `dof_fes`, `dof_residual`, standard errors, the F-statistic, and adjusted R².
+- Fix `F_kp`: it uses the first-stage degrees of freedom when the model is over-identified, applies a degrees-of-freedom correction under simple covariance (it was too high, by more with many fixed effects), and no longer depends on the order of the endogenous variables. `p_kp` changes accordingly.
+- Require Vcov 0.8.4, which contains the two changes to the rank test.
+- Improve the results display: show the dependent variable in the header, mark IV models in the title, format the F-test p-value like coefficient p-values, add a thousands separator to the number of observations, and warn on non-convergence instead of showing a Converged row (#284).
+
 ## 2.1.0
 
 - Reduce compilation latency: formula parsing, schema and model matrices go through a vector-based term pipeline, so methods compile per term type rather than per formula shape or column names (#283).
diff --git a/Project.toml b/Project.toml
index e4f039a1..72fe91d5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "FixedEffectModels"
 uuid = "9d5cd8c9-2029-5cab-9928-427838db53e3"
-version = "2.1.0"
+version = "2.2.0"
 
 [deps]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
@@ -28,7 +28,7 @@ StatsBase = "0.33, 0.34"
 StatsFuns = "0.9, 1, 2"
 StatsModels = "0.7"
 Tables = "1"
-Vcov = "0.8"
+Vcov = "0.8.4"
 julia = "1.10"
 
 
diff --git a/src/FixedEffectModel.jl b/src/FixedEffectModel.jl
index f56a427d..185fc0d4 100644
--- a/src/FixedEffectModel.jl
+++ b/src/FixedEffectModel.jl
@@ -99,7 +99,8 @@ StatsModels.formula(m::FixedEffectModel) = m.formula_schema
     dof_fes(m::FixedEffectModel)
 
 Return the number of degrees of freedom absorbed by the fixed effects (the total number of
-fixed-effect levels across all `fe(...)` dimensions).
+fixed-effect levels across all `fe(...)` dimensions, minus the levels that are redundant across
+them).
 """
 dof_fes(m::FixedEffectModel) = m.dof_fes
 
diff --git a/src/fit.jl b/src/fit.jl
index 2d097e1b..bd69ad82 100644
--- a/src/fit.jl
+++ b/src/fit.jl
@@ -323,21 +323,13 @@ function StatsAPI.fit(::Type{FixedEffectModel},
         end
     end
 
-    # Compute degrees of freedom absorbed by fixed effects.
-    # When an FE is nested within a cluster variable, it only absorbs 1 dof
-    # (its mean is not identified separately from the cluster effect),
-    # rather than the full number of groups.
-    ngroups_fes = [nunique(fe) for fe in subfes]
-    dof_fes = sum(ngroups_fes)
+    # Compute degrees of freedom absorbed by fixed effects: in total (for the adjusted R²),
+    # and for the standard errors, where an FE nested within a cluster variable does not count
+    # (its mean is not identified separately from the cluster effect).
+    dof_fes_total = dof_absorbed(subfes)
+    dof_fes = dof_fes_total
     if vcov isa Vcov.ClusterCovariance
-        dof_fes = 0
-        for i in 1:length(subfes)
-            if any(isnested(subfes[i], v.groups) for v in values(vcov_method.clusters))
-                dof_fes += 1
-            else
-                dof_fes += ngroups_fes[i]
-            end
-        end
+        dof_fes = dof_absorbed(subfes, [v.groups for v in values(vcov_method.clusters)])
     end
 
     # Compute standard error
@@ -363,8 +355,9 @@ function StatsAPI.fit(::Type{FixedEffectModel},
     if has_iv && first_stage
         Pip = Pi[(size(Pi, 1) - size(Z_res, 2) + 1):end, :]
         try 
+            # the residuals are those of the first stage, which has one regressor per instrument
             r_kp = Vcov.ranktest!(Xendo_res, Z_res, Pip,
-                              vcov_method, size(X, 2), dof_fes)
+                              vcov_method, size(Xexo, 2) + size(Z_res, 2), dof_fes)
             p_kp = chisqccdf(size(Z_res, 2) - size(Xendo_res, 2) + 1, r_kp)
             F_kp = r_kp / size(Z_res, 2)
         catch
@@ -403,5 +396,5 @@ function StatsAPI.fit(::Type{FixedEffectModel},
         esample = trues(nrows)
     end
 
-    return FixedEffectModel(coef, matrix_vcov, vcov, nclusters, esample, residuals2, augmentdf, fekeys, coef_names, response_name, formula_origin, formula_schema, contrasts, nobs, dof_, sum(ngroups_fes), dof_tstat_, rss, tss_total, F, p, iterations, converged, r2_within, F_kp, p_kp)
+    return FixedEffectModel(coef, matrix_vcov, vcov, nclusters, esample, residuals2, augmentdf, fekeys, coef_names, response_name, formula_origin, formula_schema, contrasts, nobs, dof_, dof_fes_total, dof_tstat_, rss, tss_total, F, p, iterations, converged, r2_within, F_kp, p_kp)
 end
diff --git a/src/partial_out.jl b/src/partial_out.jl
index 1ab67d6c..5615ea12 100644
--- a/src/partial_out.jl
+++ b/src/partial_out.jl
@@ -153,7 +153,7 @@ function partial_out(
         residuals .= residuals .+ m
     end
 
-    dof_fes = mapreduce(nunique, +, fes, init=0)
+    dof_fes = dof_absorbed(fes)
 
     # Return a dataframe
     out = DataFrame()
diff --git a/src/utils/fixedeffects.jl b/src/utils/fixedeffects.jl
index 50800963..2ba33947 100644
--- a/src/utils/fixedeffects.jl
+++ b/src/utils/fixedeffects.jl
@@ -153,3 +153,58 @@ function isnested(fe::FixedEffect, prefs)
     end
     return true
 end
+
+
+##############################################################################
+##
+## Degrees of freedom absorbed by fixed effects
+##
+##############################################################################
+
+# An FE nested within a cluster absorbs nothing of its own, and each intercept FE
+# after the first has as many redundant levels as the connected components it forms with an
+# earlier one (exact with two FEs, an upper bound on the degrees of freedom with more).
+function dof_absorbed(fes, clusters = ())
+    kept = [fe for fe in fes if !any(isnested(fe, g) for g in clusters)]
+    intercepts = [fe for fe in kept if fe.interaction isa UnitWeights]
+    out = sum(nunique, kept; init = 0)
+    for j in 2:length(intercepts)
+        out -= maximum(ncomponents(intercepts[i], intercepts[j]) for i in 1:(j - 1))
+    end
+    # the constant, absorbed by a nested FE, still counts
+    if isempty(intercepts) && any(fe.interaction isa UnitWeights for fe in fes)
+        out += 1
+    end
+    return out
+end
+
+# Number of connected components of the graph linking the levels of two FEs (union-find)
+function ncomponents(fe1::FixedEffect, fe2::FixedEffect)
+    parent = zeros(Int, fe1.n + fe2.n) # 0 for a level not seen yet
+    out = 0
+    @inbounds for (ref1, ref2) in zip(fe1.refs, fe2.refs)
+        i, j = Int(ref1), fe1.n + Int(ref2)
+        if parent[i] == 0
+            parent[i] = i
+            out += 1
+        end
+        if parent[j] == 0
+            parent[j] = j
+            out += 1
+        end
+        i, j = findroot!(parent, i), findroot!(parent, j)
+        if i != j
+            parent[i] = j
+            out -= 1
+        end
+    end
+    return out
+end
+
+function findroot!(parent, i)
+    @inbounds while parent[i] != i
+        parent[i] = parent[parent[i]]
+        i = parent[i]
+    end
+    return i
+end
diff --git a/test/fit.jl b/test/fit.jl
index e4ac7b75..dbb651e9 100644
--- a/test/fit.jl
+++ b/test/fit.jl
@@ -519,16 +519,16 @@ end
 	##############################################################################
 	m = @formula Sales ~ (Price ~ Pimin)
 	x = reg(df, m)
-	@test x.F_kp ≈ 52248.79247 atol = 1e-4
+	@test x.F_kp ≈ 52173.06958 atol = 1e-4
 	m = @formula Sales ~ NDI + (Price ~ Pimin)
 	x = reg(df, m)
-	@test x.F_kp ≈ 5159.812208193612 atol = 1e-4
+	@test x.F_kp ≈ 5148.59523 atol = 1e-4
 	m = @formula Sales ~ (Price ~ Pimin + CPI)
 	x = reg(df, m)
-	@test x.F_kp   ≈ 27011.44080 atol = 1e-4
+	@test x.F_kp   ≈ 26952.72028 atol = 1e-4
 	m = @formula Sales ~ (Price ~ Pimin) + fe(State)
 	x = reg(df, m)
-	@test x.F_kp ≈ 100927.75710 atol = 1e-4
+	@test x.F_kp ≈ 97490.36248 atol = 1e-4
 
 	# exactly same with ivreg2
 	m = @formula Sales ~ (Price ~ Pimin)
@@ -542,32 +542,26 @@ end
 	@test x.F_kp ≈ 2093.46609 atol = 1e-4
 	m = @formula Sales ~ (Price ~ Pimin + CPI)
 	x = reg(df, m, Vcov.robust())
-	@test x.F_kp  ≈ 16418.21196 atol = 1e-4
+	@test x.F_kp  ≈ 16406.29744 atol = 1e-4
 
 
 
-	# like in ivreg2 but += 5 difference. combination iv difference, degrees of freedom difference?
-	# iv + cluster - F_kp varies from ivreg2 about 7e-4 (SEs off by more than 2 df)
+	# exactly same with ivreg2 (SEs match with the option small)
 	m = @formula Sales ~ (Price ~ Pimin)
 	x = reg(df, m, Vcov.cluster(:State))
 	@test x.F_kp     ≈ 7249.88606 atol = 1e-4
-	# iv + cluster - F_kp varies from ivreg2 about 5e-6 (SEs off by more than 2 df)
 	m = @formula Sales ~ CPI + (Price ~ Pimin)
 	x = reg(df, m, Vcov.cluster(:State))
 	@test x.F_kp   ≈ 538.40393 atol = 1e-4
-	# Modified test values below after multiway clustering update
-	# iv + 2way clustering - F_kp matches ivreg2 (SEs match with df_add=-2)
 	m = @formula Sales ~ CPI + (Price ~ Pimin)
 	x = reg(df, m, Vcov.cluster(:State, :Year))
 	@test x.F_kp   ≈ 421.9651 atol = 1e-4
-	# multivariate iv + clustering - F_kp varies from ivreg2 about 3 (SEs off by more than 2 df)
 	m = @formula Sales ~ (Price ~ Pimin + CPI)
 	x = reg(df, m, Vcov.cluster(:State))
-	@test x.F_kp      ≈ 4080.66081 atol = 1e-4
-	# multivariate iv + multiway clustering - F_kp varies from ivreg2 about 2 (SEs off by more than 2 df)
+	@test x.F_kp      ≈ 4077.69952 atol = 1e-4
 	m = @formula Sales ~ (Price ~ Pimin + CPI)
 	x = reg(df, m, Vcov.cluster(:State, :Year))
-	@test x.F_kp  ≈ 2873.1405 atol = 1e-4
+	@test x.F_kp  ≈ 2871.05552 atol = 1e-4
 
 	############################################
 	##
@@ -900,3 +894,34 @@ end
 		@test isequal(df1[!, c], snapshot[!, c])
 	end
 end
+
+@testset "first-stage F and absorbed degrees of freedom" begin
+	df = DataFrame(CSV.File(joinpath(dirname(pathof(FixedEffectModels)), "../dataset/Cigar.csv")))
+	# F_kp with several endogenous variables does not depend on their order
+	for m in (@formula(Sales ~ (Price + NDI ~ Pimin + Pop)), @formula(Sales ~ (NDI + Price ~ Pimin + Pop)))
+		@test reg(df, m).F_kp ≈ 21.07294 atol = 1e-4
+		@test reg(df, m, Vcov.robust()).F_kp ≈ 16.90896 atol = 1e-4
+		@test reg(df, m, Vcov.cluster(:State)).F_kp ≈ 3.58888 atol = 1e-4
+	end
+	# over-identified, with fixed effects
+	m = @formula(Sales ~ (Price + NDI ~ Pimin + Pop + Pop16) + fe(State))
+	@test reg(df, m).F_kp ≈ 115.02860 atol = 1e-4
+	@test reg(df, m, Vcov.robust()).F_kp ≈ 12.82761 atol = 1e-4
+	# levels that are redundant across fixed effects do not count
+	m = @formula(Sales ~ Price + NDI + fe(State) + fe(Year))
+	x = reg(df, m)
+	@test FixedEffectModels.dof_fes(x) == 75
+	@test dof_residual(x) == 1303
+	@test stderror(x)[1] ≈ 0.0745977930 atol = 1e-7
+	@test stderror(reg(df, m, Vcov.cluster(:State)))[1] ≈ 0.1898193772 atol = 1e-7
+	@test stderror(reg(df, m, Vcov.cluster(:State, :Year)))[1] ≈ 0.2302891977 atol = 1e-7
+	# two connected components
+	df.yreg = df.Year .* 100 .+ (df.State .<= 23)
+	x = reg(df, @formula(Sales ~ Price + NDI + fe(State) + fe(yreg)))
+	@test FixedEffectModels.dof_fes(x) == 104
+	@test stderror(x)[1] ≈ 0.0757482286 atol = 1e-7
+	# a third fixed effect that is a function of the first adds no degree of freedom
+	df.region = mod.(df.State, 7)
+	x = reg(df, @formula(Sales ~ Price + NDI + fe(State) + fe(Year) + fe(region)))
+	@test FixedEffectModels.dof_fes(x) == 75
+end

From 96e58bf8fd1fdc596791bc814e161d593c68d899 Mon Sep 17 00:00:00 2001
From: Matthieu Gomez <gomez.matthieu@gmail.com>
Date: Fri, 30 Oct 2026 10:23:28 -0400
Subject: [PATCH 2/4] Use FixedEffects.components and require FixedEffects
 3.4.1

---
 CHANGELOG.md              |  2 +-
 Project.toml              |  2 +-
 src/utils/fixedeffects.jl | 33 +--------------------------------
 3 files changed, 3 insertions(+), 34 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7dbe8397..44868dfb 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,7 +4,7 @@
 
 - Fix the degrees of freedom absorbed by fixed effects: levels that are redundant across fixed effects no longer count (one per connected component), and fixed effects nested in a cluster count only once, for the constant. With two or more fixed effects, this changes `dof_fes`, `dof_residual`, standard errors, the F-statistic, and adjusted R².
 - Fix `F_kp`: it uses the first-stage degrees of freedom when the model is over-identified, applies a degrees-of-freedom correction under simple covariance (it was too high, by more with many fixed effects), and no longer depends on the order of the endogenous variables. `p_kp` changes accordingly.
-- Require Vcov 0.8.4, which contains the two changes to the rank test.
+- Require Vcov 0.8.4, which contains the two changes to the rank test, and FixedEffects 3.4.1.
 - Improve the results display: show the dependent variable in the header, mark IV models in the title, format the F-test p-value like coefficient p-values, add a thousands separator to the number of observations, and warn on non-convergence instead of showing a Converged row (#284).
 
 ## 2.1.0
diff --git a/Project.toml b/Project.toml
index 72fe91d5..ab5ad3e6 100644
--- a/Project.toml
+++ b/Project.toml
@@ -19,7 +19,7 @@ Vcov = "ec2bfdc2-55df-4fc9-b9ae-4958c2cf2486"
 
 [compat]
 DataFrames = "0.21, 0.22, 1"
-FixedEffects = "2.6, 3"
+FixedEffects = "3.4.1"
 PrecompileTools = "1"
 Reexport = "0.1, 0.2, 1"
 Statistics = "1"
diff --git a/src/utils/fixedeffects.jl b/src/utils/fixedeffects.jl
index 2ba33947..3f11a552 100644
--- a/src/utils/fixedeffects.jl
+++ b/src/utils/fixedeffects.jl
@@ -169,7 +169,7 @@ function dof_absorbed(fes, clusters = ())
     intercepts = [fe for fe in kept if fe.interaction isa UnitWeights]
     out = sum(nunique, kept; init = 0)
     for j in 2:length(intercepts)
-        out -= maximum(ncomponents(intercepts[i], intercepts[j]) for i in 1:(j - 1))
+        out -= maximum(FixedEffects.components([intercepts[i], intercepts[j]])[2] for i in 1:(j - 1))
     end
     # the constant, absorbed by a nested FE, still counts
     if isempty(intercepts) && any(fe.interaction isa UnitWeights for fe in fes)
@@ -177,34 +177,3 @@ function dof_absorbed(fes, clusters = ())
     end
     return out
 end
-
-# Number of connected components of the graph linking the levels of two FEs (union-find)
-function ncomponents(fe1::FixedEffect, fe2::FixedEffect)
-    parent = zeros(Int, fe1.n + fe2.n) # 0 for a level not seen yet
-    out = 0
-    @inbounds for (ref1, ref2) in zip(fe1.refs, fe2.refs)
-        i, j = Int(ref1), fe1.n + Int(ref2)
-        if parent[i] == 0
-            parent[i] = i
-            out += 1
-        end
-        if parent[j] == 0
-            parent[j] = j
-            out += 1
-        end
-        i, j = findroot!(parent, i), findroot!(parent, j)
-        if i != j
-            parent[i] = j
-            out -= 1
-        end
-    end
-    return out
-end
-
-function findroot!(parent, i)
-    @inbounds while parent[i] != i
-        parent[i] = parent[parent[i]]
-        i = parent[i]
-    end
-    return i
-end

From 09cda031e0253b7dc4824d65509552055946864b Mon Sep 17 00:00:00 2001
From: Matthieu Gomez <gomez.matthieu@gmail.com>
Date: Fri, 30 Oct 2026 14:09:47 -0400
Subject: [PATCH 3/4] Fix discrepancies with reghdfe found by stress tests

- Do not drop regressors or the constant because of the units of another regressor
- Do not drop regressors as collinear with fixed effects because of their mean
- Estimate models with a constant in deviations from the means
- Count slopes nested in a cluster in the degrees of freedom
- Test only linearly independent coefficients in the F-statistic
- Fix adjr2 with slope fixed effects and no constant
- Require Vcov 0.8.5
---
 CHANGELOG.md              |  8 ++++-
 Project.toml              |  2 +-
 src/FixedEffectModel.jl   |  2 +-
 src/fit.jl                | 56 +++++++++++++++++++++++------
 src/utils/basecol.jl      | 25 +++++--------
 src/utils/fixedeffects.jl |  8 ++---
 src/utils/formula.jl      |  5 +++
 src/utils/tss.jl          | 16 +++++++--
 test/fit.jl               | 75 +++++++++++++++++++++++++++++++++++++++
 9 files changed, 160 insertions(+), 37 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 44868dfb..8c8f1084 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,9 +2,15 @@
 
 ## 2.2.0
 
+- Fix regressors and the constant that were dropped as collinear because of the units of another regressor, e.g. the constant with a regressor around 1e12. Fix standard errors in the same cases.
+- Fix regressors that were dropped as collinear with fixed effects because of their mean, e.g. `year` in a panel of five years. The test now compares what remains after absorbing the fixed effects to the variation around the mean, up to the accuracy of the solver.
+- Fix the loss of precision with variables with a large mean: models with a constant are estimated in deviations from the means. This changes coefficients and standard errors only in these cases.
+- Fix the degrees of freedom of slopes nested in a cluster: they count, as in reghdfe, so standard errors, the F-statistic, and `F_kp` change in these models.
+- Fix the F-statistic when the covariance matrix is singular: it tests the linearly independent coefficients, rather than reporting values as large as 1e18.
+- Fix `adjr2` with slope fixed effects and no constant.
 - Fix the degrees of freedom absorbed by fixed effects: levels that are redundant across fixed effects no longer count (one per connected component), and fixed effects nested in a cluster count only once, for the constant. With two or more fixed effects, this changes `dof_fes`, `dof_residual`, standard errors, the F-statistic, and adjusted R².
 - Fix `F_kp`: it uses the first-stage degrees of freedom when the model is over-identified, applies a degrees-of-freedom correction under simple covariance (it was too high, by more with many fixed effects), and no longer depends on the order of the endogenous variables. `p_kp` changes accordingly.
-- Require Vcov 0.8.4, which contains the two changes to the rank test, and FixedEffects 3.4.1.
+- Require Vcov 0.8.5, which contains the two changes to the rank test and fixes standard errors of coefficients with small variances, and FixedEffects 3.4.1.
 - Improve the results display: show the dependent variable in the header, mark IV models in the title, format the F-test p-value like coefficient p-values, add a thousands separator to the number of observations, and warn on non-convergence instead of showing a Converged row (#284).
 
 ## 2.1.0
diff --git a/Project.toml b/Project.toml
index ab5ad3e6..9508abe0 100644
--- a/Project.toml
+++ b/Project.toml
@@ -28,7 +28,7 @@ StatsBase = "0.33, 0.34"
 StatsFuns = "0.9, 1, 2"
 StatsModels = "0.7"
 Tables = "1"
-Vcov = "0.8.4"
+Vcov = "0.8.5"
 julia = "1.10"
 
 
diff --git a/src/FixedEffectModel.jl b/src/FixedEffectModel.jl
index 185fc0d4..be248661 100644
--- a/src/FixedEffectModel.jl
+++ b/src/FixedEffectModel.jl
@@ -139,7 +139,7 @@ function StatsAPI.adjr2(model::FixedEffectModel, variant::Symbol=:devianceratio)
         n = nobs(model)
         dev  = deviance(model)
         dev0 = nulldeviance(model)
-        1 - (dev*(n - (has_int || has_fe(model)))) / (dev0 * max(n - k, 1))
+        1 - (dev*(n - (has_int || has_fe_intercept(model.formula)))) / (dev0 * max(n - k, 1))
     else
         throw(ArgumentError("variant must be one of :McFadden or :devianceratio"))
     end
diff --git a/src/fit.jl b/src/fit.jl
index bd69ad82..d35e3590 100644
--- a/src/fit.jl
+++ b/src/fit.jl
@@ -238,14 +238,34 @@ function StatsAPI.fit(::Type{FixedEffectModel},
     all(isfinite, Xendo) || throw(ArgumentError("Some observations for the endogenous variables are infinite"))
     all(isfinite, Z) || throw(ArgumentError("Some observations for the instrumental variables are infinite"))
 
-    iterations, converged = 0, true
-    if has_fes
-        # Save pre-demeaned data for FE recovery after estimation
-        if save_fes
-            oldy = copy(y)
-            oldX = hcat(Xexo, Xendo)
+    # Save pre-demeaned data for FE recovery after estimation
+    if has_fes && save_fes
+        oldy = copy(y)
+        oldX = hcat(Xexo, Xendo)
+    end
+
+    # When the model has a constant, work with deviations from the means. Cross-products keep
+    # their precision for variables with a large mean, and neither the convergence of the
+    # fixed-effect solver nor the collinearity test below depends on the mean.
+    has_centered_constant = has_intercept && !has_fes
+    ymean = 0.0
+    Xmeans = zeros(size(Xexo, 2) + size(Xendo, 2))
+    if has_fe_intercept || has_centered_constant
+        ymean = mean(y, weights)
+        y .-= ymean
+        for (j, x) in enumerate(vcat(eachcol(Xexo), eachcol(Xendo), eachcol(Z)))
+            has_centered_constant && j == 1 && continue
+            xmean = mean(x, weights)
+            sumsquares = sum(abs2, x)
+            x .-= xmean
+            j <= length(Xmeans) && (Xmeans[j] = xmean)
+            # only rounding errors remain in a constant variable
+            sum(abs2, x) < (1000 * eps())^2 * sumsquares && (x .= 0.0)
         end
+    end
 
+    iterations, converged = 0, true
+    if has_fes
         cols = vcat(eachcol(y), eachcol(Xexo), eachcol(Xendo), eachcol(Z))
         colnames = vcat(response_name, coefnames_exo, coefnames_endo, coefnames_iv)
         # compute 2-norm (sum of squares) for each variable 
@@ -257,9 +277,11 @@ function StatsAPI.fit(::Type{FixedEffectModel},
         # Partial out fixed effects from y, Xexo, Xendo, Z
         _, iterations, convergeds = solve_residuals!(cols, feM; maxiter = maxiter, tol = tol, progress_bar = progress_bar)
 
-        # set variables that are likely to be collinear with the fixed effects to zero
+        # set variables that are likely to be collinear with the fixed effects to zero:
+        # what remains of them is what the solver, which is accurate up to tol, cannot remove
+        collinear_tol = max(100 * tol^2, (1000 * eps(double_precision ? Float64 : Float32))^2)
         for i in 1:length(cols)
-            if sum(abs2, cols[i]) < tol * sumsquares_pre[i]
+            if sum(abs2, cols[i]) <= collinear_tol * sumsquares_pre[i]
                 if i == 1
                     @info "Dependent variable $(colnames[1]) is probably perfectly explained by fixed effects."
                 else
@@ -338,7 +360,7 @@ function StatsAPI.fit(::Type{FixedEffectModel},
     matrix_vcov = StatsAPI.vcov(vcov_data, vcov_method)
    
     # Compute Fstat
-    F = Fstat(coef, matrix_vcov, has_intercept)
+    F, rank_F = Fstat(coef, matrix_vcov, has_intercept)
    
     # dof_ is the number of estimated coefficients beyond the intercept (F-stat numerator).
     dof_ = size(X, 2) - has_intercept
@@ -348,7 +370,7 @@ function StatsAPI.fit(::Type{FixedEffectModel},
     # and minimum(nclusters)-1 = G-1 for clustering. It must NOT be decremented again for
     # the intercept (doing so reported N-K-1 instead of N-K, and G-2 instead of G-1).
     dof_tstat_ = max(1, Vcov.dof_residual(vcov_data, vcov_method))
-    p = fdistccdf(dof_, dof_tstat_, F)
+    p = fdistccdf(rank_F, dof_tstat_, F)
     
     # Compute Fstat of First Stage
     F_kp, p_kp = NaN, NaN
@@ -375,6 +397,20 @@ function StatsAPI.fit(::Type{FixedEffectModel},
     # add omitted variables and reorder IV-reclassified variables
     coef, matrix_vcov = reinsert_omitted!(coef, matrix_vcov, basis_coef, perm)
 
+    # the constant of the model in levels: y - x'b has mean ymean - Xmeans'b
+    if has_centered_constant
+        kept = findall(perm === nothing ? basis_coef : basis_coef[invperm(perm)])
+        t = .- Xmeans[kept]
+        t[1] = 1.0
+        coef[1] += ymean - Xmeans' * coef
+        matrix_vcov = Matrix(matrix_vcov)
+        row = matrix_vcov[kept, kept] * t
+        matrix_vcov[1, kept] = row
+        matrix_vcov[kept, 1] = row
+        matrix_vcov[1, 1] = t' * row
+        matrix_vcov = Symmetric(matrix_vcov)
+    end
+
     # Recover FE estimates by projecting (y - Xβ) onto the FE structure.
     # Uses pre-demeaned oldy/oldX (saved before weighting and demeaning).
     # coef has been expanded by reinsert_omitted!: omitted entries are zero
diff --git a/src/utils/basecol.jl b/src/utils/basecol.jl
index d3228a5a..11510251 100644
--- a/src/utils/basecol.jl
+++ b/src/utils/basecol.jl
@@ -6,11 +6,14 @@ Returns minus the symmetric inverse (negated so callers negate back).
 """
 function invsym!(X::Symmetric; has_intercept = false, setzeros = false, diagonal = 1:size(X, 2))
     # Per-column reference scale for the numerically-zero pivot test below.
-    tols = sweep_tolerances(X)
+    # A column is dropped when its pivot is small relative to its own sum of squares, which does
+    # not depend on the units of the other columns. A zero column has a zero tolerance.
+    sumsquares = diag(X)
+    tols = sumsquares .* sqrt(eps())
     buffer = zeros(size(X, 1))
     for j in diagonal
         d = X[j,j]
-        if setzeros && abs(d) < tols[j] * sqrt(eps())
+        if setzeros && abs(d) <= tols[j]
             X.data[1:j,j] .= 0
             X.data[j,(j+1):end] .= 0
         else
@@ -23,26 +26,14 @@ function invsym!(X::Symmetric; has_intercept = false, setzeros = false, diagonal
             X[j,j] = - 1 / d
         end
         if setzeros && has_intercept && j == 1
-            tols = sweep_tolerances(X)
+            # Relative to the sum of squares net of the constant, but no smaller than the rounding
+            # errors made in removing the constant.
+            tols = max.(diag(X) .* sqrt(eps()), sumsquares .* (1000 * eps()))
         end
     end
     return X
 end
 
-# Reference scale per column used by invsym! to flag a (post-sweep) pivot as numerically
-# zero: a column is dropped when |pivot| < scale_j * sqrt(eps). The scale is the column's
-# own diagonal, floored by a tiny multiple of the largest diagonal so the test is
-# scale-invariant. The previous absolute floor of 1 spuriously dropped genuinely
-# independent regressors whose total sum of squares was below ~sqrt(eps), while the
-# relative floor still drops all-zero columns (diag 0 -> floored to ref*eps > 0). On
-# well-scaled designs (all diagonals >= 1) this reproduces the old tolerances exactly.
-function sweep_tolerances(X::Symmetric)
-    d = diag(X)
-    ref = maximum(d; init = zero(eltype(d)))
-    floor_ = ref > 0 ? ref * eps(eltype(d)) : eps(eltype(d))
-    return max.(d, floor_)
-end
-
 """
     basis!(XX::Symmetric; has_intercept=false)
 
diff --git a/src/utils/fixedeffects.jl b/src/utils/fixedeffects.jl
index 3f11a552..dfcde3fa 100644
--- a/src/utils/fixedeffects.jl
+++ b/src/utils/fixedeffects.jl
@@ -161,18 +161,18 @@ end
 ##
 ##############################################################################
 
-# An FE nested within a cluster absorbs nothing of its own, and each intercept FE
+# An intercept FE nested within a cluster absorbs nothing of its own (a slope still does), and each intercept FE
 # after the first has as many redundant levels as the connected components it forms with an
 # earlier one (exact with two FEs, an upper bound on the degrees of freedom with more).
 function dof_absorbed(fes, clusters = ())
-    kept = [fe for fe in fes if !any(isnested(fe, g) for g in clusters)]
+    kept = [fe for fe in fes if !(fe.interaction isa UnitWeights && any(isnested(fe, g) for g in clusters))]
     intercepts = [fe for fe in kept if fe.interaction isa UnitWeights]
     out = sum(nunique, kept; init = 0)
     for j in 2:length(intercepts)
         out -= maximum(FixedEffects.components([intercepts[i], intercepts[j]])[2] for i in 1:(j - 1))
     end
-    # the constant, absorbed by a nested FE, still counts
-    if isempty(intercepts) && any(fe.interaction isa UnitWeights for fe in fes)
+    # the constant, absorbed by a nested FE, still counts if nothing else is absorbed
+    if isempty(kept) && !isempty(fes)
         out += 1
     end
     return out
diff --git a/src/utils/formula.jl b/src/utils/formula.jl
index 879cbce0..4403f7be 100644
--- a/src/utils/formula.jl
+++ b/src/utils/formula.jl
@@ -64,6 +64,11 @@ has_fe(@nospecialize(t::InteractionTerm)) = any(has_fe(x) for x in t.terms)
 has_fe(::AbstractTerm) = false
 has_fe(@nospecialize(t::FormulaTerm)) = any(has_fe(x) for x in eachterm(t.rhs))
 
+# A fixed effect that absorbs a constant, unlike a slope such as fe(id)&x
+has_fe_intercept(@nospecialize(t::AbstractTerm)) = has_fe(t)
+has_fe_intercept(@nospecialize(t::InteractionTerm)) = all(has_fe(x) for x in t.terms)
+has_fe_intercept(@nospecialize(t::FormulaTerm)) = any(has_fe_intercept(x) for x in eachterm(t.rhs))
+
 function parse_fe(@nospecialize(f::FormulaTerm))
     if has_fe(f)
         rhs_terms = eachterm(f.rhs)
diff --git a/src/utils/tss.jl b/src/utils/tss.jl
index f23013cc..221957ba 100644
--- a/src/utils/tss.jl
+++ b/src/utils/tss.jl
@@ -11,15 +11,25 @@ end
 function Fstat(coef::Vector{Float64}, matrix_vcov::AbstractMatrix{Float64}, has_intercept::Bool)
     coefF = copy(coef)
     # TODO: check I can't do better
-    length(coef) == has_intercept && return NaN
+    length(coef) == has_intercept && return NaN, 0
     if has_intercept
         coefF = coefF[2:end]
         matrix_vcov = matrix_vcov[2:end, 2:end]
     end
+    all(isfinite, matrix_vcov) || return NaN, length(coefF)
+    # Test only the linearly independent coefficients when the matrix is singular, e.g. after
+    # the Cameron, Gelbach & Miller adjustment or with fewer clusters than coefficients.
+    scale = sqrt.(diag(matrix_vcov))
+    scale[iszero.(scale)] .= 1
+    correlation = matrix_vcov ./ (scale * scale')
+    independent = basis!(Symmetric(Matrix(correlation)))
+    rank_vcov = count(independent)
+    rank_vcov == 0 && return NaN, 0
     try
-        return (coefF' * (cholesky(Symmetric(matrix_vcov)) \ coefF)) / length(coefF)
+        z = (coefF ./ scale)[independent]
+        return (z' * (cholesky(Symmetric(correlation[independent, independent])) \ z)) / rank_vcov, rank_vcov
     catch
         @info "The variance-covariance matrix is not invertible. F-statistic not computed "
-        return NaN
+        return NaN, rank_vcov
     end
 end
diff --git a/test/fit.jl b/test/fit.jl
index dbb651e9..bb34f803 100644
--- a/test/fit.jl
+++ b/test/fit.jl
@@ -925,3 +925,78 @@ end
 	x = reg(df, @formula(Sales ~ Price + NDI + fe(State) + fe(Year) + fe(region)))
 	@test FixedEffectModels.dof_fes(x) == 75
 end
+
+@testset "scale of the variables" begin
+	df = DataFrame(CSV.File(joinpath(dirname(pathof(FixedEffectModels)), "../dataset/Cigar.csv")))
+
+	# a trend with a large mean in a short panel is not collinear with the fixed effects
+	short = df[df.Year .<= 67, :]
+	short.year = short.Year .+ 1900.0
+	x = reg(short, @formula(Sales ~ NDI + year + fe(State)))
+	x0 = reg(short, @formula(Sales ~ NDI + Year + fe(State)))
+	@test coef(x) ≈ coef(x0) rtol = 1e-8
+	@test stderror(x) ≈ stderror(x0) rtol = 1e-8
+
+	# neither the constant nor a regressor is dropped because of the units of another
+	df.big = df.NDI .* 1e12
+	df.small = df.Pop .* 1e-6
+	for v in (Vcov.simple(), Vcov.robust(), Vcov.cluster(:State))
+		x = reg(df, @formula(Sales ~ big + small), v)
+		x0 = reg(df, @formula(Sales ~ NDI + Pop), v)
+		@test coef(x) ≈ coef(x0) ./ [1, 1e12, 1e-6] rtol = 1e-8
+		@test stderror(x) ≈ stderror(x0) ./ [1, 1e12, 1e-6] rtol = 1e-6
+	end
+
+	# a large mean does not affect precision, with or without weights
+	df.shifted = df.NDI .+ 1e10
+	for w in (nothing, :Pop)
+		for v in (Vcov.simple(), Vcov.robust())
+			x = reg(df, @formula(Sales ~ shifted + Price), v; weights = w)
+			x0 = reg(df, @formula(Sales ~ NDI + Price), v; weights = w)
+			@test coef(x)[2:3] ≈ coef(x0)[2:3] rtol = 1e-6
+			@test stderror(x)[2:3] ≈ stderror(x0)[2:3] rtol = 1e-6
+			@test coef(x)[1] ≈ coef(x0)[1] - 1e10 * coef(x0)[2] rtol = 1e-6
+			@test stderror(x)[1] ≈ sqrt(vcov(x0)[1, 1] - 2e10 * vcov(x0)[1, 2] + 1e20 * vcov(x0)[2, 2]) rtol = 1e-6
+		end
+	end
+
+	# a constant regressor is collinear with the constant, or with the fixed effects
+	df.constant .= 0.1
+	@test isnan(stderror(reg(df, @formula(Sales ~ constant + Price)))[2])
+	@test isnan(stderror(reg(df, @formula(Sales ~ constant + Price + fe(State))))[1])
+	# and an IV model has the same constant as a model in levels
+	x = reg(df, @formula(Sales ~ NDI + (Price ~ Pimin)))
+	df.NDIshift = df.NDI .+ 1e6
+	x2 = reg(df, @formula(Sales ~ NDIshift + (Price ~ Pimin)))
+	@test coef(x2)[1] ≈ coef(x)[1] - 1e6 * coef(x)[2] rtol = 1e-6
+	@test coef(x2)[2:3] ≈ coef(x)[2:3] rtol = 1e-6
+end
+
+@testset "slopes nested in a cluster" begin
+	df = DataFrame(CSV.File(joinpath(dirname(pathof(FixedEffectModels)), "../dataset/Cigar.csv")))
+	df.trend = df.Year .- 63.0
+	# matches reghdfe Sales NDI Pop, absorb(State##c.trend) vce(cluster State)
+	x = reg(df, @formula(Sales ~ NDI + Pop + fe(State) + fe(State)&trend), Vcov.cluster(:State); tol = 1e-11)
+	@test coef(x) ≈ [-0.0056133909, -0.0029599972] rtol = 1e-6
+	@test stderror(x) ≈ [0.0010845672, 0.0036060255] rtol = 1e-3
+end
+
+@testset "F-statistic with a singular covariance matrix" begin
+	df = DataFrame(CSV.File(joinpath(dirname(pathof(FixedEffectModels)), "../dataset/Cigar.csv")))
+	df.two = mod.(df.State, 2)
+	# two clusters identify a single linear combination of the coefficients
+	x = reg(df, @formula(Sales ~ NDI + Pop), Vcov.cluster(:two))
+	@test x.F ≈ (coef(x)[2] / stderror(x)[2])^2 rtol = 1e-6
+	# a single cluster identifies none
+	df.one .= 1
+	x = reg(df, @formula(Sales ~ NDI + Pop), Vcov.cluster(:one))
+	@test isnan(x.F)
+	@test all(isnan, stderror(x))
+end
+
+@testset "adjusted R2 with slopes and no constant" begin
+	df = DataFrame(CSV.File(joinpath(dirname(pathof(FixedEffectModels)), "../dataset/Cigar.csv")))
+	x = reg(df, @formula(Sales ~ 0 + NDI + fe(State)&Year))
+	n, k = nobs(x), 1 + 46
+	@test adjr2(x) ≈ 1 - (1 - r2(x)) * n / (n - k)
+end

From 2114e8a7eea06ce6e2de8028abdf0d949bd5cc83 Mon Sep 17 00:00:00 2001
From: Matthieu Gomez <gomez.matthieu@gmail.com>
Date: Sun, 1 Nov 2026 10:40:15 -0500
Subject: [PATCH 4/4] Demean only with fixed effects

---
 CHANGELOG.md |  1 -
 src/fit.jl   | 34 ++++++----------------------------
 test/fit.jl  | 13 -------------
 3 files changed, 6 insertions(+), 42 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8c8f1084..40a0432b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,7 +4,6 @@
 
 - Fix regressors and the constant that were dropped as collinear because of the units of another regressor, e.g. the constant with a regressor around 1e12. Fix standard errors in the same cases.
 - Fix regressors that were dropped as collinear with fixed effects because of their mean, e.g. `year` in a panel of five years. The test now compares what remains after absorbing the fixed effects to the variation around the mean, up to the accuracy of the solver.
-- Fix the loss of precision with variables with a large mean: models with a constant are estimated in deviations from the means. This changes coefficients and standard errors only in these cases.
 - Fix the degrees of freedom of slopes nested in a cluster: they count, as in reghdfe, so standard errors, the F-statistic, and `F_kp` change in these models.
 - Fix the F-statistic when the covariance matrix is singular: it tests the linearly independent coefficients, rather than reporting values as large as 1e18.
 - Fix `adjr2` with slope fixed effects and no constant.
diff --git a/src/fit.jl b/src/fit.jl
index d35e3590..815d1ec3 100644
--- a/src/fit.jl
+++ b/src/fit.jl
@@ -244,21 +244,13 @@ function StatsAPI.fit(::Type{FixedEffectModel},
         oldX = hcat(Xexo, Xendo)
     end
 
-    # When the model has a constant, work with deviations from the means. Cross-products keep
-    # their precision for variables with a large mean, and neither the convergence of the
-    # fixed-effect solver nor the collinearity test below depends on the mean.
-    has_centered_constant = has_intercept && !has_fes
-    ymean = 0.0
-    Xmeans = zeros(size(Xexo, 2) + size(Xendo, 2))
-    if has_fe_intercept || has_centered_constant
-        ymean = mean(y, weights)
-        y .-= ymean
-        for (j, x) in enumerate(vcat(eachcol(Xexo), eachcol(Xendo), eachcol(Z)))
-            has_centered_constant && j == 1 && continue
-            xmean = mean(x, weights)
+    # The fixed effects absorb the means, so remove them first: neither the convergence of the
+    # fixed-effect solver nor the collinearity test below then depends on the mean.
+    if has_fe_intercept
+        y .-= mean(y, weights)
+        for x in vcat(eachcol(Xexo), eachcol(Xendo), eachcol(Z))
             sumsquares = sum(abs2, x)
-            x .-= xmean
-            j <= length(Xmeans) && (Xmeans[j] = xmean)
+            x .-= mean(x, weights)
             # only rounding errors remain in a constant variable
             sum(abs2, x) < (1000 * eps())^2 * sumsquares && (x .= 0.0)
         end
@@ -397,20 +389,6 @@ function StatsAPI.fit(::Type{FixedEffectModel},
     # add omitted variables and reorder IV-reclassified variables
     coef, matrix_vcov = reinsert_omitted!(coef, matrix_vcov, basis_coef, perm)
 
-    # the constant of the model in levels: y - x'b has mean ymean - Xmeans'b
-    if has_centered_constant
-        kept = findall(perm === nothing ? basis_coef : basis_coef[invperm(perm)])
-        t = .- Xmeans[kept]
-        t[1] = 1.0
-        coef[1] += ymean - Xmeans' * coef
-        matrix_vcov = Matrix(matrix_vcov)
-        row = matrix_vcov[kept, kept] * t
-        matrix_vcov[1, kept] = row
-        matrix_vcov[kept, 1] = row
-        matrix_vcov[1, 1] = t' * row
-        matrix_vcov = Symmetric(matrix_vcov)
-    end
-
     # Recover FE estimates by projecting (y - Xβ) onto the FE structure.
     # Uses pre-demeaned oldy/oldX (saved before weighting and demeaning).
     # coef has been expanded by reinsert_omitted!: omitted entries are zero
diff --git a/test/fit.jl b/test/fit.jl
index bb34f803..e2d97fd6 100644
--- a/test/fit.jl
+++ b/test/fit.jl
@@ -947,19 +947,6 @@ end
 		@test stderror(x) ≈ stderror(x0) ./ [1, 1e12, 1e-6] rtol = 1e-6
 	end
 
-	# a large mean does not affect precision, with or without weights
-	df.shifted = df.NDI .+ 1e10
-	for w in (nothing, :Pop)
-		for v in (Vcov.simple(), Vcov.robust())
-			x = reg(df, @formula(Sales ~ shifted + Price), v; weights = w)
-			x0 = reg(df, @formula(Sales ~ NDI + Price), v; weights = w)
-			@test coef(x)[2:3] ≈ coef(x0)[2:3] rtol = 1e-6
-			@test stderror(x)[2:3] ≈ stderror(x0)[2:3] rtol = 1e-6
-			@test coef(x)[1] ≈ coef(x0)[1] - 1e10 * coef(x0)[2] rtol = 1e-6
-			@test stderror(x)[1] ≈ sqrt(vcov(x0)[1, 1] - 2e10 * vcov(x0)[1, 2] + 1e20 * vcov(x0)[2, 2]) rtol = 1e-6
-		end
-	end
-
 	# a constant regressor is collinear with the constant, or with the fixed effects
 	df.constant .= 0.1
 	@test isnan(stderror(reg(df, @formula(Sales ~ constant + Price)))[2])