vcov(): recommend type = "bhhh" only when it is actually available

[davidhbernstein]

Found while measuring test coverage, not from a report.

The mechanism

vcov(type = "sandwich") and vcov(type = "clustered") both need the Hessian for their bread. When that is unavailable — the fit carries none, or it is singular — both paths advised the user to switch to type = "bhhh", on the correct observation that BHHH needs no Hessian.

That reasoning is true but not sufficient. BHHH still has to invert crossprod(G), and a fit that reaches these errors is frequently flat in a parameter — which makes the outer product singular too. The advice was issued unconditionally, so on exactly those fits it sent the user to a second error.

Reproduction, sfm(model_name = "NTN") on data_gen_cs(N = 150, rand = 2):

vcov(type = "sandwich") : "... the Hessian is singular, so the bread is
                           undefined. type = "bhhh" needs no Hessian and
                           is defined here."
vcov(type = "bhhh")     : "... the outer product of gradients is singular."

Both paths now probe solve(crossprod(G)) first. BHHH is recommended only when it succeeds; when it does not, the message says so, and either way it still names the parameter responsible through .sfa_dead_parameters().

Why that fit is degenerate — verified, not assumed

NTN has no interior maximum in lambda on this sample. Profiled by continuation from the reported solution (seeding each lambda from the previous optimum, rather than from a fixed start, which otherwise fails to find sfm()'s own point at large lambda):

lambda	profile logLik
1e1	−190.3374
1e2	−188.5223
1e5	−188.1011
9.72e5 (sfm()'s)	−188.1003
1e9 … 1e15	−188.1001876 (flat to 7 dp)

Monotone increasing, converging to a finite limit it never attains. So the likelihood genuinely is flat in lambda at the reported fit, and the diagnostic's "the likelihood is flat in lambda" is independently correct rather than self-confirming.

Supporting numbers at that fit: the score column for lambda is 1.24e-10 against 1.15e+07 for the largest, and crossprod(G) has condition number 6.8e32.

How it was found, and why it had survived

Test coverage measured with covr over the full suite: 87.07% overall, but .sfa_dead_parameters() — the function that supplies the parameter name in these messages — had no coverage at all, 57 of 57 expressions never executed. The entire diagnostic path for these two covariances had never run in the suite, which is why the inconsistency between the advice and reality was never observed. The test added here is the first thing to execute it.

Tests

tests/testthat/test-vcov-bhhh-advice.R, 17 assertions, two blocks:

1. An invariant, not a string match. The message may claim BHHH is available only when vcov(type = "bhhh") actually succeeds — both observed in the same run and compared. It therefore cannot pass because the wording drifted, and cannot fail because the wording was reworded.
2. The healthy path pinned, since the availability probe is new work on a hot code path: hessian, bhhh, sandwich and clustered all still return a full-rank p × p covariance on an ordinary NHN fit.

Evidence

• Fails on main, passes here — checked against c94721a specifically: 4 failures there, 0 on this branch.
• Full suite on this branch (installed build, NOT_CRAN=true): 0 failures, 5272 passing, 43 warnings, 8 skips, 42.4 min. The warnings and skips are pre-existing; this test file contributes none of them.
• The four plm::pdata.frame: already a pdata.frame errors documented as expected in local runs did not occur in this run.
• CI on five platforms.

Scope

A55 (PR #17, merged) owns .sfa_dead_parameters() and its logic is unchanged here — this only decides whether the surrounding message recommends BHHH, and passes the diagnostic through on one additional path.

🤖 Generated with Claude Code