diff --git a/NEWS.md b/NEWS.md index 9fc8844..aa616b5 100644 --- a/NEWS.md +++ b/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 diff --git a/R/data.processing.R b/R/data.processing.R index 6c1dd87..b00ca87 100644 --- a/R/data.processing.R +++ b/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) } diff --git a/R/matrix_utils.R b/R/matrix_utils.R index 1b5e728..90d9dbd 100644 --- a/R/matrix_utils.R +++ b/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)) diff --git a/R/ttsfm.R b/R/ttsfm.R index 88cc719..c525bbb 100644 --- a/R/ttsfm.R +++ b/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) diff --git a/man/ttsfm.Rd b/man/ttsfm.Rd index db0fb06..7961fa7 100644 --- a/man/ttsfm.Rd +++ b/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{ diff --git a/tests/testthat/test-a53-missing-se-report.R b/tests/testthat/test-a53-missing-se-report.R new file mode 100644 index 0000000..0f17f47 --- /dev/null +++ b/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 + ) +})