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

NEWS.md

@@ -23,6 +23,22 @@
   singular Hessian is also caught now: it used to throw away a fit that had
   already taken minutes.
 
+* **A negative variance estimate was reported as a standard error of
+  exactly zero, and a t-value of `Inf`.** Standard errors were taken as
+  `sqrt(pmax(diag(V), 0))` at 21 places, so a parameter whose estimated
+  variance came back negative -- which is what an indefinite Hessian gives --
+  was reported as infinitely precise rather than as unidentified. On
+  `psfm(model_name = "GTRE_FML")` with `data_gen_p(t = 5, N = 40, rand = 1)`,
+  three of seven parameters have a negative variance estimate (`x1`
+  -2.83e-06, `sigv` -3.87e-04, `sigh` -3.76e-04) and all three were reported
+  with `st_err` 0 and `t-val` `Inf`, including the slope `x1 = 0.5426`. Such a
+  parameter now reports `NA`, which is what the other entry points already did
+  by taking the square root of the negative value directly. The remaining
+  standard errors on that fit are unchanged, and quantities that are not a
+  parameter's standard error -- the per-firm posterior spread in `sfm()`, where
+  a small negative is quadrature rounding on a genuinely near-zero variance --
+  still report zero.
+
 * **`sfm()` could report a fit as converged, with standard errors, at a point
   that is not a maximum.** On a flat ridge `optim()` can stop with
   convergence code 0 short of the maximum. In an RTMB scoping study on

R/esfm.R

@@ -154,7 +154,7 @@ esfm <- function(formula, data, dist = c("halfnormal", "exponential"),
   se <- rep(NA_real_, length(par))
   if (!is.null(H)) {
     V <- tryCatch(solve(H), error = function(e) NULL)
-    if (!is.null(V)) se <- sqrt(pmax(diag(V), 0))
+    if (!is.null(V)) se <- .sfa_se_from_var(diag(V))
   }
   out <- cbind(par = par, st_err = se, "t-val" = par / se)
 

R/het.R

@@ -250,7 +250,7 @@
   st_err <- rep(NA_real_, length(opt$par))
   if (isTRUE(optHessian) && !is.null(opt$hessian)) {
     vc <- tryCatch(suppressWarnings(solve(opt$hessian)), error = function(e) NULL)
-    if (!is.null(vc)) st_err <- suppressWarnings(sqrt(pmax(diag(vc), 0)))
+    if (!is.null(vc)) st_err <- suppressWarnings(.sfa_se_from_var(diag(vc)))
   }
 
   out <- matrix(NA_real_, nrow = 3L, ncol = length(opt$par),

R/kss.R

@@ -201,7 +201,7 @@
   df <- n * Tn - K - (n * L + Tn * L - L^2)
   s2 <- sum(resid^2) / max(df, 1L)
   vc <- tryCatch(s2 * solve(crossprod(Xflat)), error = function(e) NULL)
-  se <- if (is.null(vc)) rep(NA_real_, K) else sqrt(pmax(diag(vc), 0))
+  se <- if (is.null(vc)) rep(NA_real_, K) else .sfa_se_from_var(diag(vc))
 
   ## Back to the original row order.
   ## Everything above works on the (firm x period) matrix; `ord` maps it back

R/matrix_utils.R

@@ -96,9 +96,13 @@
   }
 }
 
-## log(phi(z) - z * Phi(-z)), the normal-Rayleigh density's one-sided factor.
-## Direct below z = 3; above, the two terms cancel, so use the continued
-## fraction for the Mills ratio (see notes/code_history/matrix_utils.md, A26).
+## sqrt() of a variance estimate, with a NEGATIVE estimate reported as NA.
+## pmax(v, 0) reported it as 0 -- an infinitely precise estimate (A64).
+.sfa_se_from_var <- function(v) {
+  v[!is.finite(v) | v < 0] <- NA_real_
+  sqrt(v)
+}
+
 ## Why a reported standard error is missing, in the fit's own names. An
 ## indefinite Hessian puts negative entries on the inverse diagonal (A53).
 .se_missing_report <- function(se, nms, hessian = NULL, label = "") {
@@ -159,6 +163,9 @@
   info
 }
 
+## log(phi(z) - z * Phi(-z)), the normal-Rayleigh density's one-sided factor.
+## Direct below z = 3; above, the two terms cancel, so use the continued
+## fraction for the Mills ratio (see notes/code_history/matrix_utils.md, A26).
 .log_nr_g <- function(z, K = 80L) {
   z <- as.numeric(z)
   out <- rep(NA_real_, length(z))
@@ -1508,7 +1515,7 @@
   s2 <- sum(e^2) / (n - k)
   R <- qr.R(fit$qr)
   V <- tryCatch(chol2inv(R) * s2, error = function(err) matrix(NA_real_, k, k))
-  se_b <- sqrt(pmax(diag(V), 0))
+  se_b <- .sfa_se_from_var(diag(V))
   if (!is.na(intercept_col) && intercept_col >= 1 && intercept_col <= length(se_b)) {
     se_b[intercept_col] <- NA_real_
   }

R/panel_classical.R

@@ -95,7 +95,7 @@
   s2 <- sum(fit$residuals^2) / max(df, 1L)
   XtX <- crossprod(MX)
   vc <- tryCatch(s2 * solve(XtX), error = function(e) NULL)
-  se <- if (is.null(vc)) rep(NA_real_, K) else sqrt(pmax(diag(vc), 0))
+  se <- if (is.null(vc)) rep(NA_real_, K) else .sfa_se_from_var(diag(vc))
 
   ## The firm trajectories, recovered from the residual of the frontier.
   r <- y - as.numeric(X %*% beta)
@@ -205,7 +205,7 @@
   df <- n - K - (nlevels(id) + Tn - 1L)
   s2 <- sum(resid^2) / max(df, 1L)
   vc <- tryCatch(s2 * solve(crossprod(X)), error = function(e) NULL)
-  se <- if (is.null(vc)) rep(NA_real_, K) else sqrt(pmax(diag(vc), 0))
+  se <- if (is.null(vc)) rep(NA_real_, K) else .sfa_se_from_var(diag(vc))
   names(b) <- colnames(X)
 
   list(beta = b, se = se, df = df, alpha_it = as.numeric(alpha_it),

R/psfm.R

@@ -648,7 +648,7 @@ psfm <- function(formula,
     st_err <- if (isTRUE(any(opt$hessian == 0)) | optHessian == FALSE) {
       rep(NA, length(opt$par))
     } else {
-      suppressWarnings(sqrt(pmax(diag(solve(opt$hessian)), 0)))
+      suppressWarnings(.sfa_se_from_var(diag(solve(opt$hessian))))
     }
     t_val <- opt$par / st_err
     out[1, ] <- opt$par
@@ -1344,7 +1344,7 @@ psfm <- function(formula,
         st_err <- rep(NA, length(opt$par))
       } else {
         st_err <- tryCatch(
-          suppressWarnings(sqrt(pmax(diag(solve(opt$hessian)), 0))),
+          suppressWarnings(.sfa_se_from_var(diag(solve(opt$hessian)))),
           error = function(e) rep(NA_real_, length(opt$par))
         )
       }
@@ -1462,7 +1462,7 @@ psfm <- function(formula,
               warning("OPG matrix is singular; OPG standard errors set to NA.", call. = FALSE)
               rep(NA_real_, length(opt$par))
             } else {
-              sqrt(pmax(diag(OPG_vcov), 0))
+              .sfa_se_from_var(diag(OPG_vcov))
             }
           },
           error = function(e) {
@@ -1490,7 +1490,7 @@ psfm <- function(formula,
                 rep(NA_real_, length(opt$par))
               } else {
                 sw_vcov <- H_inv %*% OPG_meat %*% H_inv
-                sqrt(pmax(diag(sw_vcov), 0))
+                .sfa_se_from_var(diag(sw_vcov))
               }
             }
           },
@@ -1784,7 +1784,7 @@ psfm <- function(formula,
     st_err <- if (isTRUE(any(opt$hessian == 0) | optHessian == FALSE)) {
       rep(NA, length(opt$par))
     } else {
-      suppressWarnings(sqrt(pmax(diag(solve(opt$hessian)), 0)))
+      suppressWarnings(.sfa_se_from_var(diag(solve(opt$hessian))))
     }
     t_val <- opt$par / st_err
     out[1, ] <- opt$par
@@ -1997,7 +1997,7 @@ psfm <- function(formula,
     st_err <- if (isTRUE(any(opt$hessian == 0)) | optHessian == FALSE) {
       rep(NA, length(opt$par))
     } else {
-      tryCatch(suppressWarnings(sqrt(pmax(diag(solve(opt$hessian)), 0))),
+      tryCatch(suppressWarnings(.sfa_se_from_var(diag(solve(opt$hessian)))),
         error = function(e) rep(NA_real_, length(opt$par))
       )
     }
@@ -2196,7 +2196,7 @@ psfm <- function(formula,
     st_err <- if (isTRUE(any(opt$hessian == 0)) | optHessian == FALSE) {
       rep(NA, length(opt$par))
     } else {
-      suppressWarnings(sqrt(pmax(diag(solve(opt$hessian)), 0)))
+      suppressWarnings(.sfa_se_from_var(diag(solve(opt$hessian))))
     }
     t_val <- opt$par / st_err
     out[1, ] <- opt$par
@@ -2371,7 +2371,7 @@ psfm <- function(formula,
     st_err <- if (isTRUE(any(opt$hessian == 0)) | optHessian == FALSE) {
       rep(NA, length(opt$par))
     } else {
-      tryCatch(suppressWarnings(sqrt(pmax(diag(solve(opt$hessian)), 0))),
+      tryCatch(suppressWarnings(.sfa_se_from_var(diag(solve(opt$hessian)))),
         error = function(e) rep(NA_real_, length(opt$par))
       )
     }
@@ -2635,7 +2635,7 @@ psfm <- function(formula,
       if (is.null(jac) || is.null(vc)) {
         c(NA_real_, NA_real_)
       } else {
-        suppressWarnings(sqrt(pmax(diag(jac %*% vc %*% t(jac)), 0)))
+        suppressWarnings(.sfa_se_from_var(diag(jac %*% vc %*% t(jac))))
       }
     }
 
@@ -2908,7 +2908,7 @@ psfm <- function(formula,
       if (is.null(jac_rep) || is.null(vc)) {
         rep(NA, length(par_names))
       } else {
-        suppressWarnings(sqrt(pmax(diag(jac_rep %*% vc %*% t(jac_rep)), 0)))
+        suppressWarnings(.sfa_se_from_var(diag(jac_rep %*% vc %*% t(jac_rep))))
       }
     }
     out[2, ] <- st_err

R/robust_divergence.R

@@ -87,5 +87,5 @@
     return(na_result)
   }
 
-  sqrt(pmax(diag(V), 0))
+  .sfa_se_from_var(diag(V))
 }

R/selsfm_kts.R

@@ -180,7 +180,7 @@
   nm <- c(paste0("t0.", colnames(X)), paste0("t1.", colnames(X)),
     "sigma_v0", "sigma_v1", "sigma_u0", "sigma_u1",
     paste0("choice.", colnames(Z)), "delta")
-  se <- tryCatch(sqrt(pmax(diag(solve(opt$hessian)), 0)),
+  se <- tryCatch(.sfa_se_from_var(diag(solve(opt$hessian))),
     error = function(e) rep(NA_real_, length(th)))
   ## The four scales are estimated as logs and reported on the natural scale,
   ## so their standard errors move with them: d(exp b)/db = exp(b).

R/sfm.R

@@ -1360,7 +1360,7 @@ sfm <- function(formula,
         rep(NA, length(opt$par))
       } else {
         vc <- tryCatch(suppressWarnings(solve(opt$hessian)), error = function(e) NULL)
-        if (is.null(vc)) rep(NA, length(opt$par)) else suppressWarnings(sqrt(pmax(diag(vc), 0)))
+        if (is.null(vc)) rep(NA, length(opt$par)) else suppressWarnings(.sfa_se_from_var(diag(vc)))
       }
     }
 

R/uhet_test.R

@@ -197,8 +197,8 @@ uhet_test <- function(object, z, data = NULL) {
     stringsAsFactors = FALSE
   )
   attr(out, "gamma") <- gam[keep]
-  attr(out, "se_corrected") <- sqrt(pmax(diag(Vc)[keep], 0))
-  attr(out, "se_naive") <- sqrt(pmax(diag(Vn)[keep], 0))
+  attr(out, "se_corrected") <- .sfa_se_from_var(diag(Vc)[keep])
+  attr(out, "se_naive") <- .sfa_se_from_var(diag(Vn)[keep])
   rownames(out) <- NULL
   out
 }

tests/testthat/test-a64-negative-variance-se.R

@@ -0,0 +1,63 @@
+## A64. A negative variance estimate was clamped by pmax(v, 0) and reported as
+## a standard error of exactly 0, with t = par/0 = Inf -- an infinitely precise
+## estimate for a direction whose estimated variance was negative. The helper
+## reports such a parameter as NA instead, and the end-to-end case is the fit
+## the defect was found on.
+
+test_that("a negative variance estimate is NA, not zero", {
+  v <- c(a = 4, b = -0.000387, c = 0.25)
+  se <- sfa:::.sfa_se_from_var(v)
+  expect_equal(se[["a"]], 2)
+  expect_equal(se[["c"]], 0.5)
+  expect_true(is.na(se[["b"]]))
+  ## The old expression returned 0 here, and par/0 is Inf.
+  expect_false(isTRUE(se[["b"]] == 0))
+  expect_identical(names(se), names(v))
+})
+
+test_that("non-finite variance entries are NA, and the rest survive", {
+  se <- sfa:::.sfa_se_from_var(c(1, NA_real_, NaN, Inf, -Inf, 9))
+  expect_equal(se[c(1, 6)], c(1, 3))
+  expect_true(all(is.na(se[2:5])))
+})
+
+test_that("a zero variance is left at zero, which is a separate question", {
+  ## Only the manufactured zeros are the defect; an exact zero that the
+  ## algebra really produced is reported as it stands.
+  expect_equal(sfa:::.sfa_se_from_var(0), 0)
+})
+
+test_that("a t-value built on the result is NA rather than Inf", {
+  par <- c(x1 = 0.542644)
+  se <- sfa:::.sfa_se_from_var(c(x1 = -3e-06))
+  expect_true(is.na(par / se))
+  expect_false(is.infinite(par / se))
+})
+
+test_that("psfm('GTRE_FML') reports no zero standard error on the A55 design", {
+  skip_on_cran()
+  ## Seed 1, N = 40, T = 5. Whether the optimizer lands on a point whose
+  ## Hessian inverse carries a negative diagonal entry is platform-dependent:
+  ## three parameters do on macOS arm64, where the defect was found, and none
+  ## do on any of CI's five platforms. So the invariant below is unconditional
+  ## and the negative-variance clause is not; the helper's own contract above
+  ## is what pins the fix down deterministically.
+  d <- as.data.frame(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))
+  f <- suppressWarnings(psfm(y_gtre ~ x1 + x2, model_name = "GTRE_FML",
+                             data = d, individual = "name"))
+  se <- f$out[, "st_err"]
+  tv <- f$out[, "t-val"]
+  expect_equal(sum(is.finite(se) & se == 0), 0L)
+  expect_equal(sum(is.infinite(tv)), 0L)
+  ## Where the variance estimates really are negative, the fit is not being
+  ## asked for a different Hessian, only to report this one honestly.
+  v <- tryCatch(diag(solve(f$opt$hessian)), error = function(e) NULL)
+  skip_if_not(!is.null(v) && any(v < 0),
+    "no negative variance estimate on this platform"
+  )
+  ## Exactly the negative-variance coordinates lose their standard error.
+  expect_true(all(is.na(se[v < 0])))
+  expect_true(all(is.finite(se[v > 0])))
+})