diff --git a/CHANGELOG.md b/CHANGELOG.md
index fdcf8482..40a0432b 100644
--- a/CHANGELOG.md
+++ b/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).
diff --git a/Project.toml b/Project.toml
index e4f039a1..9508abe0 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"
@@ -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"
 
 
diff --git a/src/FixedEffectModel.jl b/src/FixedEffectModel.jl
index f56a427d..be248661 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
 
@@ -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
diff --git a/src/fit.jl b/src/fit.jl
index 2d097e1b..815d1ec3 100644
--- a/src/fit.jl
+++ b/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
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/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 50800963..dfcde3fa 100644
--- a/src/utils/fixedeffects.jl
+++ b/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
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 e4ac7b75..e2d97fd6 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,96 @@ 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
+
+@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 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