Fix degrees-of-freedom reporting and collinearity edge cases

In order of importance:

- fit.jl: residual dof no longer subtracts the intercept twice. The reported
  dof_residual (and the t-stat p-values, confidence intervals, and F-stat
  p-value derived from it) was off by one: N-K-1 instead of N-K for
  simple/robust, and G-2 instead of G-1 for cluster-robust SEs. Standard
  errors themselves were already correct. Now matches Stata.

- basecol.jl: make rank detection in invsym! scale-invariant. The absolute
  tolerance floor of 1 dropped genuinely independent regressors whose total
  sum of squares fell below ~sqrt(eps); the new sweep_tolerances helper floors
  each column's tolerance by a relative max(diag)*eps and still drops all-zero
  columns. Bit-identical on well-scaled designs.

- fixedeffects.jl: nunique no longer counts continuous-interaction FE groups
  (fe(id)&x) whose interaction is identically zero, so dof_fes is not
  overstated for those groups.

- FixedEffectModel.jl: coeftable labels the confidence-interval columns with
  the requested level instead of a hard-coded 95%.

- basecol.jl: simplify reinsert_omitted! (reinsertion of collinear-dropped and
  IV-reclassified coefficients).

- tests: add regression tests for all of the above (residual dof, CI labels,
  continuous-slope FE dof, small-scale collinearity, reinsert_omitted!).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Author: matthieugomez

src/FixedEffectModel.jl

@@ -256,9 +256,12 @@ function StatsAPI.coeftable(m::FixedEffectModel; level = 0.95)
         coefnms = coefnms[newindex]
     end
     tt = cc ./ se
+    # Label the confidence-interval columns with the requested level, not a hard-coded 95%.
+    pct = 100 * level
+    pctstr = isinteger(pct) ? string(Int(pct)) : string(pct)
     CoefTable(
         hcat(cc, se, tt, fdistccdf.(Ref(1), Ref(StatsAPI.dof_residual(m)), abs2.(tt)), conf_int[:, 1:2]),
-        ["Estimate","Std. Error","t-stat", "Pr(>|t|)", "Lower 95%", "Upper 95%" ],
+        ["Estimate","Std. Error","t-stat", "Pr(>|t|)", "Lower $(pctstr)%", "Upper $(pctstr)%" ],
         ["$(coefnms[i])" for i = 1:length(cc)], 4)
 end
 

src/fit.jl

@@ -325,9 +325,14 @@ function StatsAPI.fit(::Type{FixedEffectModel},
     # Compute Fstat
     F = Fstat(coef, matrix_vcov, has_intercept)
 
-    # dof_ is the number of estimated coefficients beyond the intercept.
+    # dof_ is the number of estimated coefficients beyond the intercept (F-stat numerator).
     dof_ = size(X, 2) - has_intercept
-    dof_tstat_ = max(1, Vcov.dof_residual(vcov_data, vcov_method) - (has_intercept || has_fe_intercept))
+    # Residual dof for t-tests/p-values/CIs. Vcov.dof_residual already returns the final
+    # value: nobs - size(X,2) - dof_fes for simple/robust (size(X,2) already counts the
+    # intercept column, or the constant is inside dof_fes when absorbed by a fixed effect),
+    # 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)
 
     # Compute Fstat of First Stage

src/utils/basecol.jl

@@ -5,7 +5,8 @@ Generalized inverse via sweep operator.
 Returns minus the symmetric inverse (negated so callers negate back).
 """
 function invsym!(X::Symmetric; has_intercept = false, setzeros = false, diagonal = 1:size(X, 2))
-    tols = max.(diag(X), 1)
+    # Per-column reference scale for the numerically-zero pivot test below.
+    tols = sweep_tolerances(X)
     buffer = zeros(size(X, 1))
     for j in diagonal
         d = X[j,j]
@@ -17,17 +18,31 @@ function invsym!(X::Symmetric; has_intercept = false, setzeros = false, diagonal
             copy!(buffer, view(X, :, j))
             Symmetric(BLAS.syrk!('U', 'N', -1/d, buffer, one(eltype(X)), X.data))
             rmul!(buffer, 1 / d)
-            @views copy!(X.data[1:j-1,j], buffer[1:j-1])        
-            @views copy!(X.data[j, j+1:end], buffer[j+1:end])   
+            @views copy!(X.data[1:j-1,j], buffer[1:j-1])
+            @views copy!(X.data[j, j+1:end], buffer[j+1:end])
             X[j,j] = - 1 / d
         end
         if setzeros && has_intercept && j == 1
-            tols = max.(diag(X), 1)
+            tols = sweep_tolerances(X)
         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)
 
@@ -197,28 +212,19 @@ function reinsert_omitted!(
     perm::Union{Nothing, Vector{Int}}        # permutation from IV reclassification, or nothing
 )
     if !all(basis_coef)
+        kept = findall(basis_coef)
         newcoef = zeros(length(basis_coef))
+        newcoef[kept] = coef
+
         newmatrix_vcov = fill(NaN, (length(basis_coef), length(basis_coef)))
-        newindex = [searchsortedfirst(cumsum(basis_coef), i) for i in 1:length(coef)]
-        for i in eachindex(newindex)
-            newcoef[newindex[i]] = coef[i]
-            for j in eachindex(newindex)
-                newmatrix_vcov[newindex[i], newindex[j]] = matrix_vcov[i, j]
-            end
-        end
+        newmatrix_vcov[kept, kept] = Matrix(matrix_vcov)
         coef = newcoef
         matrix_vcov = Symmetric(newmatrix_vcov)
     end
     if perm !== nothing
-        _invperm = invperm(perm)
-        coef = coef[_invperm]
-        newmatrix_vcov = zeros(size(matrix_vcov))
-        for i in 1:size(newmatrix_vcov, 1)
-            for j in 1:size(newmatrix_vcov, 1)
-                newmatrix_vcov[i, j] = matrix_vcov[_invperm[i], _invperm[j]]
-            end
-        end
-        matrix_vcov = Symmetric(newmatrix_vcov)
+        invp = invperm(perm)
+        coef = coef[invp]
+        matrix_vcov = Symmetric(Matrix(matrix_vcov)[invp, invp])
     end
     return coef, matrix_vcov
 end

src/utils/fixedeffects.jl

@@ -52,8 +52,20 @@ end
 
 function nunique(fe::FixedEffect)
     seen = falses(fe.n)
-    @inbounds for ref in fe.refs
-        seen[ref] = true
+    if fe.interaction isa UnitWeights
+        @inbounds for ref in fe.refs
+            seen[ref] = true
+        end
+    else
+        # Continuous-slope fixed effect (e.g. fe(id)&x): a group whose interaction is
+        # identically zero contributes no column to the design and absorbs no degree of
+        # freedom (the demeaning sets its scale to 0, SolverCPU.scale!). Count only groups
+        # with at least one nonzero interaction value, so dof_fes is not overstated.
+        @inbounds for i in eachindex(fe.refs, fe.interaction)
+            if !iszero(fe.interaction[i])
+                seen[fe.refs[i]] = true
+            end
+        end
     end
     count(seen)
 end

test/collinearity.jl

@@ -1,4 +1,5 @@
-using DataFrames, CSV, FixedEffectModels, Random, Statistics, Test
+using DataFrames, CSV, FixedEffectModels, Random, Statistics, Test, LinearAlgebra
+using FixedEffectModels: reinsert_omitted!
 
 @testset "collinearity_with_fixedeffects" begin
 
@@ -44,3 +45,37 @@ using DataFrames, CSV, FixedEffectModels, Random, Statistics, Test
   rr = reg(df, @formula(x1 ~ x2))
   @test all(!isnan, stderror(rr))
 end
+
+@testset "reinsert omitted" begin
+  coef = [10.0, 20.0]
+  vcov = Symmetric([1.0 0.2; 0.2 4.0])
+  basis_coef = BitVector([true, false, true])
+  perm = [2, 1, 3]
+
+  newcoef, newvcov = reinsert_omitted!(coef, vcov, basis_coef, perm)
+  newvcov_matrix = Matrix(newvcov)
+
+  @test newcoef == [0.0, 10.0, 20.0]
+  @test isnan(newvcov_matrix[1, 1])
+  @test isnan(newvcov_matrix[1, 2])
+  @test isnan(newvcov_matrix[1, 3])
+  @test isnan(newvcov_matrix[2, 1])
+  @test newvcov_matrix[2, 2] == 1.0
+  @test newvcov_matrix[2, 3] == 0.2
+  @test isnan(newvcov_matrix[3, 1])
+  @test newvcov_matrix[3, 2] == 0.2
+  @test newvcov_matrix[3, 3] == 4.0
+end
+
+@testset "small-scale independent regressor is kept" begin
+  # The rank test is scale-invariant: a genuinely independent regressor whose total
+  # sum of squares is below sqrt(eps) used to be dropped as collinear (NaN stderror).
+  df = DataFrame(CSV.File(joinpath(dirname(pathof(FixedEffectModels)), "../dataset/Cigar.csv")))
+  df.tiny = df.Price .* 3e-8
+  m = reg(df, @formula(Sales ~ tiny))
+  mP = reg(df, @formula(Sales ~ Price))
+  @test !isnan(stderror(m)[2])
+  # scaling a regressor by c scales its coefficient by 1/c and leaves the t-stat invariant
+  @test coef(m)[2] * 3e-8 ≈ coef(mP)[2] rtol = 1e-4
+  @test coef(m)[2] / stderror(m)[2] ≈ coef(mP)[2] / stderror(mP)[2] rtol = 1e-4
+end

test/fit.jl

@@ -778,3 +778,33 @@ end
 	df_ivinf = DataFrame(y = [1.0, 2.0, 3.0], x = [0.0, 1.0, 2.0], z = [1.0, Inf, 3.0])
 	@test_throws ArgumentError reg(df_ivinf, @formula(y ~ (x ~ z)))
 end
+
+@testset "residual degrees of freedom and CI labels" begin
+	df = DataFrame(CSV.File(joinpath(dirname(pathof(FixedEffectModels)), "../dataset/Cigar.csv")))
+
+	# The intercept is counted exactly once: residual dof is N - K, not N - K - 1.
+	@test dof_residual(reg(df, @formula(Sales ~ Price))) == 1378
+	@test dof_residual(reg(df, @formula(Sales ~ 0 + Price))) == 1379
+	@test dof_residual(reg(df, @formula(Sales ~ CPI + (Price ~ Pimin)))) == 1377
+	# A fixed effect absorbs the constant once (dof = N - slopes - #levels).
+	@test dof_residual(reg(df, @formula(Sales ~ Price + fe(State)))) == 1380 - 1 - length(unique(df.State))
+	# Cluster-robust dof is G - 1 (clusters minus one), not G - 2.
+	@test dof_residual(reg(df, @formula(Sales ~ Price), Vcov.cluster(:State))) == length(unique(df.State)) - 1
+
+	# coeftable labels the confidence-interval columns with the requested level.
+	m = reg(df, @formula(Sales ~ Price))
+	@test coeftable(m).colnms[end-1:end] == ["Lower 95%", "Upper 95%"]
+	@test coeftable(m; level = 0.90).colnms[end-1:end] == ["Lower 90%", "Upper 90%"]
+	@test coeftable(m; level = 0.99).colnms[end-1:end] == ["Lower 99%", "Upper 99%"]
+end
+
+@testset "continuous-slope FE degrees of freedom" begin
+	# A group whose continuous interaction is identically zero contributes no column
+	# to the design and must not be counted in dof_fes.
+	df = DataFrame(id = repeat(1:3, inner = 4),
+	               z = Float64[3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8],
+	               y = Float64[5, 2, 7, 3, 1, 8, 4, 9, 6, 2, 5, 3])
+	df.x = Float64.(df.id .!= 1)          # interaction identically zero for id == 1
+	m = reg(df, @formula(y ~ z + fe(id)&x))
+	@test FixedEffectModels.dof_fes(m) == 2
+end