Use ucminf() with forward differences in the clmm() submodel fits

NEWS.md

@@ -12,6 +12,7 @@ If you read this from a place other than <https://mc-stan.org/projpred/news/inde
 * Faster clustering of the reference model draws for a large number of observations: `kmeans()` runs on the coordinates of the draws in their low-dimensional span, which gives the same clusters. For example, `varsel()` with 50 000 observations, 4 000 draws and 20 clusters is about 4.5 times faster. (GitHub: #576)
 * Faster projections onto submodels without group-level or additive terms: the projected draws (or clusters) share one model matrix, and for the Gaussian family with the identity link, one linear solve gives the submodel fits for all draws. For example, `varsel()` with forward search for a Gaussian model with 50 observations and 200 predictor terms is about 7 times faster. (GitHub: #574)
 * Faster computation of the predictive variance of the reference model for the Gaussian and Student-t families. This matters for a large number of observations and many projected draws. For example, `project()` with 10 000 observations and 4 000 projected draws is about 7 times faster. (GitHub: #575)
+* Faster `ordinal::clmm()` submodel fits with a single group-level intercept: projpred now uses `ucminf()` with forward differences, and refits with `nlminb()` if needed. For example, `varsel()` for a cumulative model with 300 observations is about 1.5 to 2 times faster. (GitHub: #596)
 * Faster `lme4::glmer()` submodel fits: projpred now uses BOBYQA in both optimization stages, which changes the fits only within the optimizer tolerance. For example, `varsel()` for a Bernoulli model with 300 observations and two grouping factors is about 2 times faster. (GitHub: #592)
 * Faster projections onto ordinal submodels without group-level terms (`brms::cumulative()`, `rstanarm::stan_polr()`): one C++ call fits all projected draws, with results equal to those of `MASS::polr()` up to its convergence tolerance. For example, `cv_varsel()` with `validate_search = TRUE` is about 23 times faster for a model with 6 categories. (GitHub: #583)
 * Faster projections onto categorical submodels without group-level terms (`brms::categorical()`): one C++ call fits all projected draws, with results equal to those of `nnet::multinom()` up to its convergence tolerance. For example, `cv_varsel()` with `validate_search = TRUE` is about 8 times faster for a model with 4 categories. (GitHub: #583)

R/divergence_minimizers.R

@@ -1839,6 +1839,15 @@ fit_cumul_mlvl <- function(formula, data, family, weights, ...) {
       methods::formalArgs(stats::nlminb),
       methods::formalArgs(stats::optim))
   )]
+  # clmm() has no analytic gradient, so its default optimizer nlminb() spends
+  # most objective evaluations on finite-difference gradients.  ucminf() with
+  # forward differences needs fewer evaluations and gives the same fits within
+  # the optimizer tolerance.  clmm() can use ucminf() only for a single
+  # `(1 | <group>)` term.  Use it then, unless the caller chose an optimizer:
+  bars <- reformulas::findbars(formula)
+  use_ucminf <- is.null(dot_args[["method"]]) &&
+    is.null(dot_args[["control"]]) && length(bars) == 1 &&
+    identical(bars[[1]][[2]], 1)
   # Since ordinal::clmm() sometimes seems to have issues with formulas lacking
   # an intercept, add it here (note that this is a quite "hacky" solution;
   # perhaps reformulate() could be used instead):
@@ -1850,20 +1859,55 @@ fit_cumul_mlvl <- function(formula, data, family, weights, ...) {
     link_nm <- "probit"
   }
   # Call the submodel fitter:
-  fitobj <- do.call(ordinal::clmm, c(
-    list(formula = formula,
-         data = data,
-         weights = quote(projpred_internal_w_aug),
-         contrasts = NULL,
-         Hess = FALSE,
-         model = FALSE,
-         link = link_nm),
-    dot_args
-  ))
-  # Needed for the ordinal:::predict.clm() workaround (the value `"negative"` is
-  # the default, see `?ordinal::clm.control`):
-  fitobj$control$sign.location <- "negative"
-  return(fitobj)
+  fit_clmm <- function(dot_args) {
+    fitobj <- do.call(ordinal::clmm, c(
+      list(formula = formula,
+           data = data,
+           weights = quote(projpred_internal_w_aug),
+           contrasts = NULL,
+           Hess = FALSE,
+           model = FALSE,
+           link = link_nm),
+      dot_args
+    ))
+    # Needed for the ordinal:::predict.clm() workaround (the value `"negative"`
+    # is the default, see `?ordinal::clm.control`):
+    fitobj$control$sign.location <- "negative"
+    return(fitobj)
+  }
+  if (use_ucminf) {
+    # If ucminf() fails or does not converge, refit with the nlminb() default:
+    return(fit_with_fallback(
+      function() fit_clmm(c(dot_args, list(method = "ucminf",
+                                           grad = "forward"))),
+      function() fit_clmm(dot_args),
+      accept = check_conv_s
+    ))
+  }
+  return(fit_clmm(dot_args))
+}
+
+# Return `fit_first()` if it gives no error and `accept()` returns `TRUE` for
+# its result; then signal its warnings and messages again.  Otherwise return
+# `fit_second()`, without the warnings and messages of `fit_first()`.
+fit_with_fallback <- function(fit_first, fit_second, accept) {
+  conds <- list()
+  keep_cond <- function(cnd, restart) {
+    conds[[length(conds) + 1L]] <<- cnd
+    invokeRestart(restart)
+  }
+  res <- try(withCallingHandlers(
+    fit_first(),
+    warning = function(w) keep_cond(w, "muffleWarning"),
+    message = function(m) keep_cond(m, "muffleMessage")
+  ), silent = TRUE)
+  if (!inherits(res, "try-error") && isTRUE(accept(res))) {
+    for (cnd in conds) {
+      if (inherits(cnd, "warning")) warning(cnd) else message(cnd)
+    }
+    return(res)
+  }
+  return(fit_second())
 }
 
 # Use nnet::multinom() to fit submodels for the brms::categorical() family:
@@ -2159,6 +2203,16 @@ check_conv_s <- function(fit_s) {
   } else if (inherits(fit_s, "multinom")) {
     return(fit_s$convergence == 0)
   } else if (inherits(fit_s, "clmm")) {
+    if (is.null(fit_s$optRes$objective)) {
+      # No `objective` (which nlminb() returns), so ucminf() ran.  Its codes:
+      # 1 (small gradient), 2 (small step) and 4 (zero step in the line
+      # search, which forward differences can cause near the optimum).
+      # `gradient` is the central-difference gradient after the fit; the
+      # limit is absolute, as `grtol` of ucminf():
+      return(fit_s$optRes$convergence %in% c(1, 2, 4) &&
+               all(is.finite(fit_s$gradient)) &&
+               max(abs(fit_s$gradient)) < 1e-3)
+    }
     return(fit_s$optRes$convergence == 0)
   } else if (inherits(fit_s, "polr")) {
     return(fit_s$convergence == 0)

R/projpred-package.R

@@ -81,7 +81,12 @@
 #'     + For the augmented-data projection: [ordinal::clmm()] (the returned
 #'     object inherits from class `clmm`) for the [brms::cumulative()] family,
 #'     [mclogit::mblogit()] (the returned object inherits from class `mmblogit`)
-#'     for the [brms::categorical()] family.
+#'     for the [brms::categorical()] family. For [ordinal::clmm()] with a single
+#'     group-level term `(1 | <group>)`, \pkg{projpred} uses
+#'     `method = "ucminf"` and `grad = "forward"` (see
+#'     [ordinal::clmm.control()]). If this fit fails or does not converge,
+#'     \pkg{projpred} refits with the default settings of [ordinal::clmm()]. An
+#'     argument `method` or `control` in `...` replaces this default.
 #' * Submodel without multilevel but additive terms: [mgcv::gam()] (the returned
 #' object inherits from class `gam`).
 #' * Submodel with multilevel and additive terms: [gamm4::gamm4()] (within

tests/testthat/test_fit_cumul_mlvl.R

@@ -0,0 +1,111 @@
+# The optimizer of the ordinal::clmm() submodel fits and their convergence
+# check. These tests need no reference model fits, so they also run on CRAN.
+
+clmm_test_data <- function() {
+  set.seed(5678)
+  dat <- data.frame(x = rnorm(150), g = gl(15, 10))
+  eta <- 0.6 * dat$x + rnorm(15, sd = 0.7)[dat$g] + rlogis(150)
+  dat$y <- factor(cut(eta, c(-Inf, -1, 0, 1, Inf), labels = FALSE),
+                  levels = 1:4, ordered = TRUE)
+  dat
+}
+
+fit_clmm_test <- function(dat, formula = y ~ x + (1 | g), ...) {
+  fit_cumul_mlvl(formula = formula, data = dat,
+                 family = list(link = "logit"), weights = rep(1, nrow(dat)),
+                 ...)
+}
+
+test_that("fit_cumul_mlvl() uses ucminf() with forward differences", {
+  skip_if_not_installed("ordinal")
+  dat <- clmm_test_data()
+  fit <- fit_clmm_test(dat)
+  expect_s3_class(fit, "clmm")
+  expect_identical(fit$control$method, "ucminf")
+  expect_identical(fit$control$optCtrl$grad, "forward")
+  expect_true(check_conv_s(fit))
+  # A `method` or `control` from the caller replaces the default:
+  fit_nlminb <- fit_clmm_test(dat, method = "nlminb")
+  expect_identical(fit_nlminb$control$method, "nlminb")
+  expect_true(check_conv_s(fit_nlminb))
+  fit_ctrl <- fit_clmm_test(dat, control = ordinal::clmm.control())
+  expect_identical(fit_ctrl$control$method, "nlminb")
+  # Both optimizers give the same fit within the optimizer tolerance:
+  expect_equal(fit$logLik, fit_nlminb$logLik, tolerance = 1e-6)
+  expect_equal(fit$coefficients, fit_nlminb$coefficients, tolerance = 1e-3)
+})
+
+test_that("fit_cumul_mlvl() uses nlminb() if clmm() cannot use ucminf()", {
+  skip_if_not_installed("ordinal")
+  # clmm() can use ucminf() only for a single `(1 | <group>)` term:
+  expect_no_warning(
+    fit <- fit_clmm_test(clmm_test_data(), formula = y ~ x + (x | g)),
+    message = "ucminf"
+  )
+  expect_identical(fit$control$method, "nlminb")
+  expect_false(is.null(fit$optRes$objective))
+})
+
+test_that("fit_cumul_mlvl() refits with nlminb() if ucminf() does not converge", {
+  skip_if_not_installed("ordinal")
+  local_mocked_bindings(check_conv_s = function(fit_s) {
+    !identical(fit_s$control$method, "ucminf")
+  })
+  fit <- fit_clmm_test(clmm_test_data())
+  expect_identical(fit$control$method, "nlminb")
+})
+
+test_that("fit_with_fallback() keeps or drops the first fit and its conditions", {
+  fit_first <- function() {
+    warning("first warning")
+    message("first message")
+    1
+  }
+  fit_second <- function() 2
+  # Accepted: the result and the conditions of the first fit.
+  expect_message(
+    expect_warning(
+      res <- fit_with_fallback(fit_first, fit_second, function(x) TRUE),
+      "first warning"
+    ),
+    "first message"
+  )
+  expect_identical(res, 1)
+  # Not accepted: the second fit, without the conditions of the first.
+  expect_no_condition(
+    res <- fit_with_fallback(fit_first, fit_second, function(x) FALSE)
+  )
+  expect_identical(res, 2)
+  # An error in the first fit:
+  expect_no_condition(
+    res <- fit_with_fallback(function() stop("failed"), fit_second,
+                             function(x) TRUE)
+  )
+  expect_identical(res, 2)
+})
+
+test_that("check_conv_s() accepts ucminf() fits by code and gradient", {
+  skip_if_not_installed("ordinal")
+  fit <- fit_clmm_test(clmm_test_data())
+  fit$gradient[] <- 1e-5
+  for (code in c(1, 2, 4)) {
+    fit$optRes$convergence <- code
+    expect_true(check_conv_s(fit), info = paste("code", code))
+  }
+  # Evaluation limit or no start:
+  for (code in c(3, -2)) {
+    fit$optRes$convergence <- code
+    expect_false(check_conv_s(fit), info = paste("code", code))
+  }
+  # A zero step far from the optimum:
+  fit$optRes$convergence <- 4
+  fit$gradient[1] <- 0.5
+  expect_false(check_conv_s(fit))
+  fit$gradient[1] <- NA
+  expect_false(check_conv_s(fit))
+  # An nlminb() fit is judged by its own code, also if `method` is "ucminf"
+  # (clmm() uses nlminb() for models that ucminf() cannot fit):
+  fit_nlminb <- fit_clmm_test(clmm_test_data(), method = "nlminb")
+  fit_nlminb$control$method <- "ucminf"
+  expect_true(check_conv_s(fit_nlminb))
+})