Skip to content

Report a negative variance estimate as NA, not as a zero standard error (A64) - #62

Merged
davidhbernstein merged 3 commits into
mainfrom
fix/negative-variance-se
Oct 27, 2026
Merged

davidhbernstein merged 3 commits into
mainfrom
fix/negative-variance-se

Conversation

@davidhbernstein

@davidhbernstein davidhbernstein commented Oct 23, 2026 •

Copy link
Copy Markdown
Owner

A negative variance estimate was reported as a standard error of exactly zero

Standard errors were taken as sqrt(pmax(diag(V), 0)) at 21 places in R/.
The clamp means a parameter whose estimated variance came back negative --
which is exactly what an indefinite Hessian gives -- was reported as
st_err = 0 and t-val = par/0 = Inf: an infinitely precise estimate for the
one direction the covariance machinery has nothing to say about. The oldest of
those sites dates to CRAN 1.0.4 (2026-02-04).

Reachable on the first seed tried, on this machine

psfm(model_name = "GTRE_FML") on A55's own design -- data_gen_p(t = 5, N = 40, rand = 1, sig_u = 1, sig_v = 0.3, sig_r = 0.2, sig_h = 0.4, cons = 0.5, beta1 = 0.5, beta2 = 0.5, eta = 0.1), 3.1 seconds to fit:

              (Intercept)        x1         x2        sigr       sigv       sigh      sigu
par              0.413756  0.542644   0.475742  5.2403e-08   0.409347   0.567744  0.841880
diag(solve(H))  1.006e-05 -2.83e-06  8.535e-06     4.0e-18  -3.869e-04 -3.762e-04 1.403e-05
st_err           0.003171  0.000000   0.002921     2.0e-09   0.000000   0.000000  0.003745
t-val             130.474       Inf    162.843     26.2017        Inf        Inf   224.800

Three of seven parameters, one an ordinary slope (x1 = 0.5426), reported with
a standard error of zero. Over seeds 1-10 of that design: 1 shows this, 4
return all-NA standard errors because solve() fails outright (already
handled), 5 are ordinary. It is not rare arithmetic -- it is what an indefinite
Hessian always produces.

The change

.sfa_se_from_var() (R/matrix_utils.R) reports a negative or non-finite
variance estimate as NA and takes the square root otherwise. All 21 sites
call it: psfm (10), panel_classical and uhet_test (2 each), and one each
in sfm, esfm, het, kss, selsfm_kts, robust_divergence and
matrix_utils. On seed 1 the three negative coordinates become NA, the
other four standard errors are unchanged to every digit
, and no Inf
t-value survives. This is what the unclamped entry points (ttsfm, zsfm,
lcsfm, copsfm, ivsfm, selsfm) already did by taking the square root of
the negative value directly.

vcov()'s singular-Hessian fallback builds a diagonal from the reported
standard errors (R/sfareg_methods.R:209), so it now carries NA for those
parameters instead of a zero variance. Its !all(is.na(...)) guard is
unaffected.

Deliberately left alone

  • Quantities that are not a parameter's standard error. sfm()'s per-firm
    posterior spread sd_exp_u_hat is sqrt(pmax(E[exp(-2u)] - E[exp(-u)]^2, 0))
    by quadrature, where a small negative is rounding on a genuinely near-zero
    variance and 0 is the right report; changing it would blank the column for
    well-behaved firms. Likewise the 1e-3 curvature floor and the
    double.xmin density floor, which are not variances.
  • An exactly-zero variance, which keeps reporting 0. The defect is the
    clamp manufacturing zeros out of negatives.
  • sigr in the table above, at 5.24e-08 with variance 4e-18, standard
    error 2e-09 and t-value 26.2 -- a collapsed parameter reported as
    comfortably significant. That is A55's open territory (name the parameter
    that died), it survives this change untouched, and it is arguably more
    misleading than the Inf because it looks like an ordinary number.

Tests

tests/testthat/test-a64-negative-variance-se.R, 14 assertions: the helper's
contract, that a t-value built on it is NA rather than Inf, and the
end-to-end GTRE_FML fit above. On main the file gives 7 failures
(4 errors for the missing helper, plus the three end-to-end assertions:
se == 0 count, infinite-t count, and NA on exactly the negative
coordinates); on this branch 14 pass, 0 fail on macOS.

The end-to-end case is platform-dependent, and the first CI run proved it.
It originally asserted expect_gt(sum(v < 0), 0) -- that this design really
does reach a negative variance estimate. It does here (3 of 7 parameters) and
on none of CI's five platforms, so all five jobs failed on that one
assertion and nothing else (FAIL 1 each, 44m-2h06m) while the fix itself was
correct. Which point the 3-stage minimiser lands on, and the sign of a small
eigenvalue there, is BLAS-dependent arithmetic; the fixed seed makes the data
identical everywhere, not the fit. That premise is now gated on
skip_if_not(any(v < 0), ...) with the solve() wrapped in tryCatch, so the
clause skips off this Mac. The invariant stays unconditional -- no finite
standard error of exactly 0, no infinite t-value -- and the helper's four
deterministic unit tests, whose input is constructed rather than optimised to,
carry the proof of the fix.

Full local suite against the merged branch build, after #63 landed and this
branch took main (macOS, load_package = "installed", NOT_CRAN=true,
TESTTHAT_CPUS=2): pass 5531, fail 0, warn 48, skip 9. For the record, CI's
warning counts on the pre-merge branch ranged 182-225 across the five platforms
against identical code, which is why 48 is reported rather than compared.

The merge with main, and one thing it fixes

NEWS.md conflicted as expected -- both entries kept, A53's first since it is
already on main. R/matrix_utils.R conflicted because both branches inserted
their helper at the same place: between .log_nr_g's three-line comment and
.log_nr_g itself, so each had silently orphaned that comment onto its own
function, and main carries the same wart from #63. The merge puts
.sfa_se_from_var() and .se_missing_report() above the comment and
re-attaches it, so landing this also fixes that on main.
There is no same-invocation main baseline for the warning count this round,
so 48 is reported rather than compared; the change cannot add warnings, since
.sfa_se_from_var() emits none and sqrt() of NA is silent where
sqrt() of a negative was not.

Found while scoping A53 item 5, which claimed sqrt(diag(solve(h))) "destroys
every one of them together" and prescribed per-parameter reporting. R's
sqrt() is already element-wise (sqrt(c(4, -1, 9)) is 2 NaN 3), so that
item is refuted and recorded as such; the clamped sibling expression is the
real defect.

🤖 Generated with Claude Code

davidhbernstein and others added 3 commits October 22, 2026 21:44
Standard errors were taken as sqrt(pmax(diag(V), 0)) at 21 places, so a
parameter whose estimated variance came back negative -- what an indefinite
Hessian gives -- was reported as st_err 0 and t-val par/0 = Inf: an infinitely
precise estimate for the one direction the covariance machinery has nothing to
say about. The oldest such site dates to CRAN 1.0.4.

Reachable on the first seed tried. psfm("GTRE_FML") on A55's design
(data_gen_p(t = 5, N = 40, rand = 1, ...)) gives three of seven parameters a
negative variance estimate -- x1 -2.83e-06, sigv -3.87e-04, sigh -3.76e-04 --
and reported all three with a standard error of 0 and an infinite t-value,
including the slope x1 = 0.5426.

.sfa_se_from_var() reports a negative or non-finite variance estimate as NA and
takes the square root otherwise; all 21 sites call it. The other four standard
errors on that fit are unchanged to every digit. This matches what the
unclamped entry points already did by taking the square root directly.

Left alone: quantities that are not a parameter's standard error, where a small
negative is rounding on a genuinely near-zero variance and 0 is the right
report -- sfm()'s per-firm posterior spread, the curvature floor, the density
floor; and an exactly-zero variance, since the defect is the clamp
manufacturing zeros out of negatives.

Found while scoping A53 item 5, which claimed sqrt(diag(solve(h))) loses every
negative entry together. R's sqrt() is element-wise, so that item is refuted;
the clamped sibling expression is the real defect. horserace/FUNCTIONALITY_GAPS.md A64.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The end-to-end case asserted expect_gt(sum(v < 0), 0) -- that the
GTRE_FML fit on this design really does land on a point whose Hessian
inverse has a negative diagonal entry. It does on macOS arm64, where the
defect was found (3 of 7 parameters: x1, sigv, sigh), and on none of
CI's five platforms, so all five failed on that one assertion and
nothing else.

The optimizer's landing point is not what the fix changes, so the
premise does not belong in an expectation. The invariant that does hold
everywhere -- no finite standard error of exactly 0, no infinite t-value
-- stays unconditional; the clause about which coordinates lose their
standard error is now gated on this platform actually producing a
negative estimate, and solve() is wrapped so an exactly singular Hessian
skips rather than errors. The helper's own contract, four deterministic
tests above, is what pins the fix down.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Two conflicts, both routine, and one of them worth a note.

NEWS.md: both branches added an entry under the top "## Bug fixes".
Both kept, A53's first since it is already on main.

R/matrix_utils.R: both branches inserted their new helper at the same
place -- between .log_nr_g's three-line comment and .log_nr_g itself --
so each had silently orphaned that comment onto its own function, and
main carries the same wart from #63. Resolved by putting
.sfa_se_from_var() and .se_missing_report() above the comment and
re-attaching it to .log_nr_g, which this merge therefore also fixes on
main.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
@davidhbernstein
davidhbernstein merged commit 5163c44 into main Oct 27, 2026
9 of 10 checks passed
davidhbernstein added a commit that referenced this pull request Oct 27, 2026
The test hard-coded the landing point of a single fit on macOS arm64: the
names of the two parameters that lose a standard error, and the count
2 of 6. Which point the optimizer stops at, and the sign of a small
eigenvalue there, is BLAS-dependent -- exactly what made #62 fail on all
five CI platforms a few days ago, and exactly what notes from that round
said not to do again.

Every expectation is now derived from the fit in front of it: the
reported parameters must be the ones whose standard error is non-finite,
n_missing must be that count, n_par the length of the vector, and
"indefinite" must hold iff a negative eigenvalue was read. The regime is
a skip condition, not an assertion, in both directions -- the first test
skips where every standard error is finite, since there is then nothing
for the report to say, and the false-positive guard skips where the
N = 300 fit is not healthy, since its claim is about a healthy fit.
expect_silent() is gone too: it demanded total silence from a fit, which
is a stronger claim than the test needs.

12 of 12 pass against the branch build and 2 of 4 fail against main, so
it still discriminates.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
davidhbernstein added a commit that referenced this pull request Oct 27, 2026
… (A65) (#66)

* Name the standard errors the other five entry points could not compute (A65)

zsfm(), lcsfm(), copsfm(), ivsfm() and selsfm() each took their standard
errors as suppressWarnings(sqrt(diag(solve(opt$hessian)))). On an
indefinite Hessian the negative-variance coordinates come back NaN, and
suppressWarnings() discarded R's own "NaNs produced" -- the only signal
there was. Nothing set a flag, and the convergence code read 0.

Reproduced on zsfm("ZISF"), data_gen_cs(N = 25, rand = 8, sig_u = 1,
sig_v = .3, cons = .5, beta1 = .5, beta2 = .5, mu = .5, a = 5), 0.5s to
fit: eigenvalues 3.48e6 ... -1474, two of six standard errors NaN -- one
of them the slope x2, whose true value is 0.5 -- conv 0, no warnings, no
flags, and print()/summary() said nothing beyond printing NaN.

All five now call .se_missing_report(), the model-agnostic helper #63
added for ttsfm(), on the REPORTED standard errors and names, so
copsfm()'s and ivsfm()'s delta-method scales carry through correctly.
$se_unresolved is appended to the results list rather than re-listed
against its name vector. solve() is wrapped too: the
sum(colMeans(hessian)) == 0 test caught only an all-zero Hessian, so a
singular but non-zero one errored and discarded a finished fit.

The test file fails on main (3 failures: no warning, no flag) and passes
17 of 17 here, and includes a false-positive guard -- the same design at
N = 300 must stay silent with no flag.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

* Derive A65's expectations from the fit instead of pinning one platform's

The test hard-coded the landing point of a single fit on macOS arm64: the
names of the two parameters that lose a standard error, and the count
2 of 6. Which point the optimizer stops at, and the sign of a small
eigenvalue there, is BLAS-dependent -- exactly what made #62 fail on all
five CI platforms a few days ago, and exactly what notes from that round
said not to do again.

Every expectation is now derived from the fit in front of it: the
reported parameters must be the ones whose standard error is non-finite,
n_missing must be that count, n_par the length of the vector, and
"indefinite" must hold iff a negative eigenvalue was read. The regime is
a skip condition, not an assertion, in both directions -- the first test
skips where every standard error is finite, since there is then nothing
for the report to say, and the false-positive guard skips where the
N = 300 fit is not healthy, since its claim is about a healthy fit.
expect_silent() is gone too: it demanded total silence from a fit, which
is a stronger claim than the test needs.

12 of 12 pass against the branch build and 2 of 4 fail against main, so
it still discriminates.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

* NEWS: conv 0 on the A65 reproduction is machine-specific

The entry reported "conv 0, and zero warnings" as the headline. The zero
warnings is the defect and holds; conv 0 does not travel. On all five CI
platforms the same fit reports convergence 52, so there the convergence
code carries a signal and only the missing standard errors are
unexplained. On this Mac it reports 0, which is the fully silent case.

The CI run also establishes something the local run could not: the
warning fired on macOS, ubuntu devel, ubuntu release and windows -- the
only failing expectation there was the convergence code -- so the
indefinite-Hessian regime is reached on four of the five platforms.
ubuntu oldrel-1 does not reach it, which is why the rewritten test skips
rather than asserts.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

* Extend A65 to psfm(), which has a seventh site and a third spelling

psfm() carries the same two defects at seven Hessian-diagonal sites,
covering GTRE and TRE (which share one), GTRE_Z, TRE_Z, GTRE_FML,
TFE_WMLE, TFE and FD. Four of the seven left the inverse unguarded, so a
singular but non-zero Hessian discarded a finished PANEL fit -- the most
expensive kind this package produces. Their guard tests
any(hessian == 0), which catches a zero ENTRY rather than singularity.

The count recorded two days ago was six, and it was wrong. FD spells the
same quantity sqrt(diag(qr.solve(opt$hessian))) -- qr.solve, and
unclamped -- which makes it invisible to A64's sweep, to A65's first
pass, and to the .sfa_se_from_var(diag(solve(...))) pattern that found
the other six. It surfaced only because a replace-all asserted 4
occurrences and got 5; the sweep was then redone on $hessian itself,
spelling-agnostically, which is what should have happened first.

Deliberately excluded after that sweep: psfm()'s two other Hessian
inversions (the sandwich H_inv, the delta-method vc), which already
guard solve() and take their standard errors through a Jacobian, so the
report's eigenvalue narrative does not map onto them as directly.

$se_unresolved is appended at eight assemblies -- GTRE and TRE share a
site but keep separate results lists -- and the insertion was scoped to
each enclosing if (model_name ...) block so it could not leak into
GTRE_SEQ1/SSFE/PL80/GTRE_SEQ2, where .sm does not exist.

Two instances came out of the existing suite, which is why the local
warning count moved:

  * psfm("FD") at test-psfm.R:23 is GENUINE -- 3 of 5 standard errors
    missing, smallest eigenvalue -0.00119 against a largest of 17400. So
    the site every sweep missed is reachable in the package's own tests,
    and had been reporting three unexplained NAs there all along.
  * lcsfm("LCM_CN") at test-a52-opt-shape.R:91 was a FALSE POSITIVE of
    mine. That call passes optHessian = FALSE, so there is no Hessian by
    request; LCM_CN is the one branch whose standard-error test does not
    mention optHessian, and so the one place the guard was missing. Now
    guarded.

Test uses A64's own design, psfm("GTRE_FML") on data_gen_p(t = 5,
N = 40, rand = 1, ...) -- 3 seconds, three of seven parameters with a
negative variance estimate -- gated on the regime from the start, since
A64's end-to-end case had to be gated after all five CI platforms
disagreed with this Mac.

Full local suite on the final code: pass 5556, fail 0, warn 49, skip 9,
against a 48-warning baseline -- the one addition being the genuine FD
report.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

---------

Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant