Fix some discrepancies with reghdfe/ivreghdfe

CHANGELOG.md

@@ -1,5 +1,17 @@
 # Changelog
 
+## 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 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.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
 
 - 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).

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"
@@ -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"
@@ -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.5"
 julia = "1.10"
 
 

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
 
@@ -138,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

src/fit.jl

@@ -238,14 +238,26 @@ 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
+
+    # 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 .-= mean(x, weights)
+            # 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 +269,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
@@ -323,21 +337,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
@@ -346,7 +352,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
@@ -356,15 +362,16 @@ 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
     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 +410,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

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()

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)
 

src/utils/fixedeffects.jl

@@ -153,3 +153,27 @@ function isnested(fe::FixedEffect, prefs)
     end
     return true
 end
+
+
+##############################################################################
+##
+## Degrees of freedom absorbed by fixed effects
+##
+##############################################################################
+
+# 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 !(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 nothing else is absorbed
+    if isempty(kept) && !isempty(fes)
+        out += 1
+    end
+    return out
+end

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)

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