Name the standard errors ttsfm() could not compute, instead of saying nothing (A53 item 2)

NEWS.md

@@ -2,6 +2,27 @@
 
 ## Bug fixes
 
+* **`ttsfm()` said nothing when it could not compute a standard error, and
+  errored outright when the Hessian was exactly singular.** On a fit whose
+  Hessian is indefinite the negative-variance coordinates lose their standard
+  error, which is correct -- no standard error exists for them there -- but
+  the fit reported `NaN` in those cells and no warning of any kind. Measured on
+  `data_gen_cs(N = 300, rand = 3002, sig_u = .5, sig_v = .2, cons = 5,
+  beta1 = 1.5, beta2 = 2, a = 5, mu = .1)` fitted with `model_name = "TTHN"`:
+  428 seconds, one negative eigenvalue (`-0.00312` against a largest of
+  `30578.7`), two of six standard errors `NaN`, zero warnings. Seed 3049 of the
+  same design gives three of six, also silently.
+
+  `ttsfm()` now names the parameters whose standard error could not be
+  computed, reports the Hessian's smallest and largest eigenvalue, and says
+  explicitly that an indefinite Hessian here is *either* a point that is not a
+  maximum *or* a boundary solution -- where the Hessian in the full parameter
+  space need not be positive definite -- rather than asserting one reading. The
+  result is recorded in `$se_unresolved` and printed by `print()` and
+  `summary()`, since a warning does not survive `saveRDS()`. `solve()` on a
+  singular Hessian is also caught now: it used to throw away a fit that had
+  already taken minutes.
+
 * **`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/data.processing.R

@@ -416,6 +416,21 @@ print.sfareg <- function(x, ...) {
       " both the split and any\n  level from coef() as unidentified in this",
       " sample. See ?psfm.\n", sep = "")
   }
+  ## A standard error that could not be computed at all, named for the same
+  ## reason as the flags above: a warning does not survive saveRDS().
+  if (!is.null(x$se_unresolved)) {
+    cat("NOTE: no standard error could be computed for ",
+      paste(x$se_unresolved$parameters, collapse = ", "), " (",
+      x$se_unresolved$n_missing, " of ", x$se_unresolved$n_par,
+      " parameters)",
+      if (isTRUE(x$se_unresolved$indefinite)) {
+        ": the Hessian at the reported\n  point is not positive definite"
+      } else {
+        ""
+      },
+      ". Either the point is not a maximum or it is a\n  boundary solution --",
+      " the two readings differ. Inspect $se_unresolved.\n", sep = "")
+  }
   invisible(NULL)
 }
 

R/matrix_utils.R

@@ -99,6 +99,66 @@
 ## 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).
+## 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 = "") {
+  bad <- which(!is.finite(se))
+  if (!length(bad)) {
+    return(NULL)
+  }
+  if (length(nms) != length(se)) nms <- paste0("[", seq_along(se), "]")
+  ev <- if (is.matrix(hessian) && all(is.finite(hessian))) {
+    tryCatch(eigen(hessian, symmetric = TRUE, only.values = TRUE)$values,
+      error = function(e) NULL
+    )
+  } else {
+    NULL
+  }
+  info <- list(
+    parameters = nms[bad], n_missing = length(bad), n_par = length(se),
+    min_eigenvalue = if (is.null(ev)) NA_real_ else min(ev),
+    max_eigenvalue = if (is.null(ev)) NA_real_ else max(ev),
+    indefinite = if (is.null(ev)) NA else any(ev < 0)
+  )
+  ## The conditioning, when it could be read, and otherwise nothing: a claim
+  ## about eigenvalues that were never computed is worse than silence.
+  cond <- if (isTRUE(info$indefinite)) {
+    paste0(
+      " The Hessian at the reported point is not positive definite -- its ",
+      "smallest eigenvalue is ", signif(info$min_eigenvalue, 3), " against a ",
+      "largest of ", signif(info$max_eigenvalue, 3), " -- so the estimated ",
+      "variance of those parameters is negative and no standard error exists ",
+      "for them there."
+    )
+  } else if (is.null(ev)) {
+    paste0(
+      " The Hessian could not be read, so the estimated variance of those ",
+      "parameters is unavailable rather than merely imprecise."
+    )
+  } else {
+    paste0(
+      " The Hessian's smallest eigenvalue is ", signif(info$min_eigenvalue, 3),
+      " against a largest of ", signif(info$max_eigenvalue, 3), ", so it is ",
+      "positive definite but close enough to singular that inverting it lost ",
+      "those entries."
+    )
+  }
+  info$message <- paste0(
+    label, "no standard error could be computed for ",
+    paste(info$parameters, collapse = ", "), " (", info$n_missing, " of ",
+    info$n_par, " parameters).", cond,
+    " This is EITHER a point that is not a maximum OR a boundary solution, ",
+    "where the Hessian in the full parameter space need not be positive ",
+    "definite; the two readings differ and are not distinguished by this ",
+    "check. The standard errors that ARE reported are conditional on this ",
+    "point. For a scale that has collapsed onto its lower bound see ",
+    "tt_boundary_report(); for a covariance that needs no Hessian, refit with ",
+    "keep_objective = TRUE and use vcov(type = \"bhhh\"). Inspect ",
+    "$se_unresolved."
+  )
+  info
+}
+
 .log_nr_g <- function(z, K = 80L) {
   z <- as.numeric(z)
   out <- rep(NA_real_, length(z))

R/ttsfm.R

@@ -510,10 +510,14 @@ ttsfm <- function(formula,
     }
 
     if (optHessian == TRUE) {
+      ## solve() on a singular Hessian errors outright, which threw away a fit
+      ## that had already taken minutes; report it as missing instead (A53).
       st_err <- if (isTRUE(as.numeric(sum(colMeans(opt$hessian))) == 0)) {
         rep(NA, length(opt$par))
       } else {
-        suppressWarnings(sqrt(diag(solve(opt$hessian))))
+        tryCatch(suppressWarnings(sqrt(diag(solve(opt$hessian)))),
+          error = function(e) rep(NA_real_, length(opt$par))
+        )
       }
     }
     ## Natural scale and distinct names; see .tt_report() above.
@@ -606,9 +610,23 @@ ttsfm <- function(formula,
 
     metric.ne.res <- metrics.ne(p = opt$par)
 
+    ## A fit that could not report a standard error used to say nothing about
+    ## it, however long it had taken to get there (A53).
+    .sm <- if (optHessian == TRUE) {
+      .se_missing_report(out["st_err", ], colnames(out), opt$hessian,
+        paste0("ttsfm(model_name = \"", model_name, "\"): ")
+      )
+    } else {
+      NULL
+    }
+    if (!is.null(.sm)) warning(.sm$message, call. = FALSE)
+
     results <- list(t(out), c(opt), End.Time, start_v, model_name, formula, out["par", ], out["st_err", ], out["t-val", ], metric.ne.res, call)
     class(results) <- "sfareg"
     names(results) <- c("out", "opt", "total_time", "start_v", "model_name", "formula", "coefficients", "std.errors", "t.values", "metrics", "call")
+    ## Appended, not rebuilt: re-listing eleven values against eleven names to
+    ## add one field is how a misalignment gets introduced (see sfm.R).
+    if (!is.null(.sm)) results$se_unresolved <- .sm
     ## Rows actually used; bread() scales by this.
     results$nobs <- length(as.numeric(Y))
     results$anchor_resid <- .sfa_anchor_resid(data_i_vars, Y)
@@ -819,10 +837,14 @@ ttsfm <- function(formula,
     }
 
     if (optHessian == TRUE) {
+      ## solve() on a singular Hessian errors outright, which threw away a fit
+      ## that had already taken minutes; report it as missing instead (A53).
       st_err <- if (isTRUE(as.numeric(sum(colMeans(opt$hessian))) == 0)) {
         rep(NA, length(opt$par))
       } else {
-        suppressWarnings(sqrt(diag(solve(opt$hessian))))
+        tryCatch(suppressWarnings(sqrt(diag(solve(opt$hessian)))),
+          error = function(e) rep(NA_real_, length(opt$par))
+        )
       }
     }
     ## Natural scale and distinct names; see .tt_report() above.
@@ -921,9 +943,23 @@ ttsfm <- function(formula,
 
     metric.hn.res <- tryCatch(metrics.hn(p = opt$par), error = function(e) NULL)
 
+    ## A fit that could not report a standard error used to say nothing about
+    ## it, however long it had taken to get there (A53).
+    .sm <- if (optHessian == TRUE) {
+      .se_missing_report(out["st_err", ], colnames(out), opt$hessian,
+        paste0("ttsfm(model_name = \"", model_name, "\"): ")
+      )
+    } else {
+      NULL
+    }
+    if (!is.null(.sm)) warning(.sm$message, call. = FALSE)
+
     results <- list(t(out), c(opt), End.Time, start_v, model_name, formula, out["par", ], out["st_err", ], out["t-val", ], metric.hn.res, call)
     class(results) <- "sfareg"
     names(results) <- c("out", "opt", "total_time", "start_v", "model_name", "formula", "coefficients", "std.errors", "t.values", "metrics", "call")
+    ## Appended, not rebuilt: re-listing eleven values against eleven names to
+    ## add one field is how a misalignment gets introduced (see sfm.R).
+    if (!is.null(.sm)) results$se_unresolved <- .sm
     ## Rows actually used; bread() scales by this.
     results$nobs <- length(as.numeric(Y))
     results$anchor_resid <- .sfa_anchor_resid(data_i_vars, Y)

man/ttsfm.Rd

@@ -91,6 +91,7 @@ Where \eqn{\sigma_u} or \eqn{\sigma_w} carries determinants there is no single s
 \item{call}{The matched call.}
 \item{nobs}{The number of observations used in the fit.}
 \item{objective}{The likelihood closure, when \code{keep_objective = TRUE} (\code{"TTNE"} and \code{"TTHN"} only).}
+\item{se_unresolved}{Present, for \code{"TTNE"} and \code{"TTHN"} fitted with \code{optHessian = TRUE}, only when at least one reported standard error could not be computed. A list naming those parameters (\code{parameters}, \code{n_missing} of \code{n_par}), the Hessian's \code{min_eigenvalue} and \code{max_eigenvalue}, whether it is \code{indefinite} (\code{NA} when the Hessian could not be read), and the \code{message} the fit warned with. An indefinite Hessian means the estimated variance of those parameters is negative, so no standard error exists for them at that point; it is \emph{either} a point that is not a maximum \emph{or} a boundary solution, where the Hessian in the full parameter space need not be positive definite, and this component does not distinguish the two -- see \code{\link{tt_boundary_report}} for a collapsed scale. The standard errors that are reported are conditional on the same point. Printed by \code{print()} and \code{summary()}, because a warning does not survive \code{saveRDS()}.}
 }
 \examples{
 \donttest{

tests/testthat/test-a53-missing-se-report.R

@@ -0,0 +1,106 @@
+## A53 item 2. A ttsfm() fit whose Hessian is indefinite loses the standard
+## errors of the negative-variance coordinates and used to say nothing about
+## it, after minutes of fitting. Driven on hand-built matrices rather than
+## through a fit: the fits that reach this state sit on flat ridges, which is
+## what does not reproduce across platforms. The live reproduction is the slow
+## test at the bottom.
+
+mk_h <- function(ev) {
+  set.seed(11)
+  Q <- qr.Q(qr(matrix(stats::rnorm(length(ev)^2), length(ev), length(ev))))
+  Q %*% diag(ev) %*% t(Q)
+}
+
+test_that("nothing is reported when every standard error is finite", {
+  expect_null(sfa:::.se_missing_report(c(a = 1, b = 2), c("a", "b"), diag(2)))
+})
+
+test_that("the missing parameters are named, with the conditioning", {
+  h <- mk_h(c(30578.7, 300, 100, 10, 1, -0.00311802))
+  se <- c(5.03, 1.48, 1.99, NaN, 0.51, NaN)
+  nms <- c("(Intercept)", "x1", "x2", "sigma_v", "sigma_u", "sigma_w")
+  r <- sfa:::.se_missing_report(se, nms, h, "ttsfm(model_name = \"TTHN\"): ")
+  expect_equal(r$parameters, c("sigma_v", "sigma_w"))
+  expect_equal(r$n_missing, 2L)
+  expect_equal(r$n_par, 6L)
+  expect_true(r$indefinite)
+  expect_equal(r$min_eigenvalue, -0.00311802, tolerance = 1e-9)
+  expect_equal(r$max_eigenvalue, 30578.7, tolerance = 1e-9)
+  expect_match(r$message, "^ttsfm\\(model_name = \"TTHN\"\\): ")
+  expect_match(r$message, "sigma_v, sigma_w", fixed = TRUE)
+  expect_match(r$message, "not positive definite", fixed = TRUE)
+  expect_match(r$message, "-0.00312", fixed = TRUE)
+})
+
+test_that("the message does not decide between the two readings", {
+  ## A boundary optimum has no reason to have a positive definite Hessian in
+  ## the full parameter space, so "not a maximum" must not be asserted alone.
+  r <- sfa:::.se_missing_report(c(1, NaN), c("a", "b"), mk_h(c(10, -1)))
+  expect_match(r$message, "EITHER a point that is not a maximum OR a boundary",
+    fixed = TRUE
+  )
+  expect_match(r$message, "tt_boundary_report()", fixed = TRUE)
+})
+
+test_that("an unreadable Hessian is reported as unavailable, not as indefinite", {
+  r <- sfa:::.se_missing_report(c(1, NA_real_), c("a", "b"), NULL)
+  expect_true(is.na(r$indefinite))
+  expect_true(is.na(r$min_eigenvalue))
+  expect_match(r$message, "could not be read", fixed = TRUE)
+  expect_false(grepl("not positive definite", r$message, fixed = TRUE))
+  r2 <- sfa:::.se_missing_report(c(1, NA_real_), c("a", "b"),
+    matrix(c(1, NA, NA, 1), 2, 2)
+  )
+  expect_true(is.na(r2$indefinite))
+})
+
+test_that("a positive definite but near-singular Hessian gets its own wording", {
+  r <- sfa:::.se_missing_report(c(1, NaN), c("a", "b"), mk_h(c(1e12, 1e-14)))
+  expect_false(r$indefinite)
+  expect_match(r$message, "close enough to singular", fixed = TRUE)
+  expect_false(grepl("not positive definite", r$message, fixed = TRUE))
+})
+
+test_that("unnamed standard errors fall back to positional labels", {
+  r <- sfa:::.se_missing_report(c(1, NaN), character(0), diag(c(1, -1)))
+  expect_equal(r$parameters, "[2]")
+})
+
+test_that("print() reports the flag, because a warning does not survive saveRDS", {
+  x <- structure(
+    list(
+      out = matrix(1, 1, 1), formula = y ~ x, total_time = 1,
+      opt = list(value = 1), model_name = "TTHN",
+      se_unresolved = list(parameters = c("sigma_v", "sigma_w"), n_missing = 2L,
+                           n_par = 6L, indefinite = TRUE)
+    ),
+    class = "sfareg"
+  )
+  out <- paste(capture.output(print(x)), collapse = " ")
+  expect_match(out, "no standard error could be computed for sigma_v, sigma_w",
+    fixed = TRUE
+  )
+  expect_match(out, "2 of 6 parameters", fixed = TRUE)
+  expect_match(out, "not positive definite", fixed = TRUE)
+})
+
+test_that("a TTHN fit with an indefinite Hessian warns and carries the flag", {
+  skip_on_cran()
+  skip_if_not(nzchar(Sys.getenv("SFA_TEST_SLOW")),
+    "one TTHN fit is ~400s; set SFA_TEST_SLOW to run"
+  )
+  ## Seed 3002 of A53's design: the fit's own Hessian has one negative
+  ## eigenvalue (-0.00312) and two of six standard errors come back NaN.
+  d <- data_gen_cs(N = 300, rand = 3002, sig_u = .5, sig_v = .2, cons = 5,
+                   beta1 = 1.5, beta2 = 2, a = 5, mu = .1)
+  expect_warning(
+    f <- ttsfm(y_tthn ~ x1 + x2, model_name = "TTHN", data = d),
+    "no standard error could be computed"
+  )
+  expect_false(is.null(f$se_unresolved))
+  expect_true(f$se_unresolved$indefinite)
+  expect_equal(f$se_unresolved$n_missing, sum(!is.finite(f$out[, "st_err"])))
+  expect_match(paste(capture.output(print(f)), collapse = " "),
+    "no standard error could be computed", fixed = TRUE
+  )
+})