[avehtari]

This PR was assisted by Claude.

Summary

For brms::cumulative(), rstanarm::stan_polr() and brms::categorical() reference models, the augmented-data projection fits each projected draw with MASS::polr() (fit_cumul()) or nnet::multinom() (fit_categ()). Both build the model frame of the augmented data (C copies of the n observations) and minimize the negative log-likelihood with BFGS; polr() evaluates its objective and gradient in R and finds start values with glm.fit(). After each projection, init_submodl() predicts each draw's fit with predict() (which builds the model frame again) and computes the cross-entropy with a loop over observations and draws. In profiles of cv_varsel(validate_search = TRUE), the fitter took 87 % of the time for an ordinal model with 6 categories (N 800, 11 predictors) and 82 % for a categorical model with 4 categories (N 400, 20 predictors).

This PR has two commits.

1. Ordinal submodels.
   • cumul_native_multi(), which divmin_augdat() calls before the draw-wise loop for cumulative-link submodels without group-level terms, builds the model matrix once and writes the augmented weights of each draw as an n x C matrix of category weights. The new C++ function cumul_fit_multi_c() then fits all draws.
   • The objective and the latent distributions are transcribed from MASS::polr(): its negative log-likelihood, including the truncation of the latent-scale arguments to [-100, 100], and its pfun / dfun. All links of polr() are supported.
   • The fitter is new: Newton steps with the analytic gradient and Hessian, damped if the Hessian is not positive definite, and a backtracking line search that keeps the thresholds increasing. It starts from the model without predictors. It computes the upper-category probabilities from the survival function, so that they do not round to zero; for the cloglog link, polr() itself fails in such cases with "initial value in 'vmmin' is not finite".
   • cumul_native_multi() returns polr objects with the elements that predict.polr(), coef(), as.matrix.projection() and the convergence check use.
   • subprd_augdat() first calls the new subprd_augdat_direct(): for polr fits (and multinom fits) of the same submodel, it builds the model matrix for newdata once and computes the linear predictors from the coefficients ($\zeta_j - \eta$ for polr, $X \beta_k$ for multinom) instead of one predict() call per draw followed by a round trip through the probabilities.
   • family$ce_ptwise() of the augmented-data families sums over the categories with one rowSums() call.
2. Categorical submodels. categ_native_multi() does the same for multinomial-logit submodels without group-level terms, with the new C++ function categ_fit_multi_c(). It minimizes the objective of nnet::multinom() (no weight decay, as in fit_categ(); the objective is convex) with Newton steps with the analytic gradient and Hessian, from the intercepts of the model without predictors. nnet::multinom() itself uses BFGS. categ_native_multi() returns copies of one multinom object from nnet::multinom(maxit = 0) per projection, with the weights, objective, deviance, AIC, fitted values, residuals and convergence code of each draw, so predict.multinom(), coef() and as.matrix.projection() work as before. With rang = 0, nnet::multinom() draws no random numbers, so the random number stream is unchanged.

In both cases, divmin_augdat() uses the draw-wise fitter as before if ... contains arguments that fit_cumul() / fit_categ() pass on to polr() / optim() / multinom(), if the model frame has offsets or missing values, or if the augmented data do not have the expected layout. The internal option projpred.augdat_native = FALSE switches the native fitters off.

Results

The fitters change from BFGS to Newton iterations to convergence, so the submodel fits agree with those of polr() and multinom() up to their convergence tolerance. Against polr() with reltol = 1e-15, coefficients and thresholds agree within 2e-6; against multinom(), the objective was at most 9e-8 lower and the predicted probabilities agreed within 4e-6. The direct predictions differ from the probability-based ones of master by at most 2e-13 (relative), and the cross-entropy is bitwise identical. In the bench below, the rankings, the suggested sizes and the CV ranking proportions are identical to master, and the largest relative difference in the ELPD estimates is 4.5e-6.

Bench: winequality-red, random subsample N = 800, 11 standardized predictors, quality with 6 categories (brms cumulative("logit"), rstanarm stan_polr()); simulated data with 4 categories, N = 400, 20 correlated predictors (brms categorical()); forward search; vs = varsel() (nterms_max 11 or 15), cvt10 = cv_varsel(validate_search = TRUE, nloo = 10, refit_prj = TRUE) (nterms_max 8 or 10); 1 CPU per run, all versions in the same job on one node (CSC Roihu, AMD EPYC 9965). Times in seconds.

model	call	master	this PR	speedup
brms cumulative()	vs	451.5	27.4	16.5
brms cumulative()	cvt10	3712.3	159.7	23.2
rstanarm stan_polr()	vs	437.6	23.4	18.7
rstanarm stan_polr()	cvt10	3577.1	152.7	23.4
brms categorical()	vs	212.9	27.2	7.8
brms categorical()	cvt10	1194.8	148.6	8.0

After this PR, the C++ fits take most of the time for ordinal models; for categorical models, brms's inverse link (inv_link_categorical(), called via family$linkinv) takes about a third.

Tests

The test suite passes with NOT_CRAN=true: 261 tests, 55 887 expectations, 0 failures, 0 errors, the same 4 warnings as master. Tested with this change applied to master 73cf0ee (projpred 2.11.0.9000), R 4.6.1, MASS 7.3-65, nnet 7.3-20, rstanarm 2.32.2, brms 2.23.0, on CSC Roihu (Linux, AMD EPYC).

New tests:

• tests/testthat/test_cumul_native.R: cumul_fit_multi_c() gives the polr() estimates (tolerance 1e-5) and a deviance not above that of polr() for the logistic, probit, loglog and cauchit links; for cloglog, the numeric gradient is zero at the fit; several weight matrices in one call give the same fits as one at a time. cumul_native_multi() returns polr objects with the names of polr(), the same predicted probabilities, a converged state, and a named empty coefficient vector without predictors; NULL for polr() arguments in ..., with the option set to FALSE, and for permuted augmented data.
• tests/testthat/test_augdat_predict.R: subprd_augdat_direct() gives the linear predictors and dimension names of the probability-based path for polr (logistic, probit, cloglog; factor predictor; no predictors) and multinom fits; NULL for fits of different submodels, mixed classes and other fits; ce_ptwise() gives exactly the values of the loop, also with zero probabilities.
• tests/testthat/test_categ_native.R: categ_fit_multi_c() gives the estimates of a tightly converged multinom() (tolerance 1e-5) and an objective not above it; categ_native_multi() returns multinom objects with the names, levels, objective, predicted probabilities and fitted values of multinom(), also without predictors; NULL for multinom() arguments in ..., with the option set to FALSE, and for permuted augmented data.

The new tests fail on master (the functions do not exist there).

Overlap with other PRs

• Leave the left-out observation out of the fold projections in validated LOO CV #580: this PR and Leave the left-out observation out of the fold projections in validated LOO CV #580 change the same lines in divmin_augdat(). The second PR to merge needs a rebase. The resolution keeps the removal of the observations with weight 0 from Leave the left-out observation out of the fold projections in validated LOO CV #580 first. The native paths of this PR then get the data without these observations. Both native paths return early, so each of them must set the attribute projpred_lost_lvls on outdmin before return(outdmin), as divmin_augdat() does at its end. Leave the left-out observation out of the fold projections in validated LOO CV #580 also multiplies the augmented-data weights (projpred_ws_aug) by the observation weights in proj_to_submodl(). The native paths use projpred_ws_aug as their only weights, so they need no change for this part.

[avehtari]

I added three commits to this PR. They fix two problems of the native cumulative-link fitter.

Start values (5f57373): the fitter starts the thresholds at the link quantiles of the cumulative weight shares. With the cauchit link, a lowest or highest category with less than about 0.3 % of the weight gives a start threshold outside [-100, 100]. The objective is then infinite at the start values. cl_fit_one() stopped with the C++ error max(): object has no elements (the Hessian was empty). In one test, the data have 200 observations and 0.2 % of the weight in the lowest category. This error occurred in that test with the arguments that cumul_native_multi() passes (offset 0, thresh = 1e-12, maxit = 100). The fitter now limits the start thresholds to [-50, 50], so each category has a positive probability at the start. In that test, the fit converges with a deviance of 396.55; polr() gives 397.48. For the other links, the start thresholds are within +-19, so their fits do not change. Cauchit fits with a lowest or highest cumulative weight share between about 0.3 % and 0.6 % now start at -50 or 50; their results can change within the convergence tolerance.

Non-finite input (fb3f43a, ecb6f00): with a non-finite predictor value, the linear predictor at the start values is Inf * 0 = NaN. std::min(100, NaN) returns 100, so cl_eval() treated the NaN as a value at the bound. The fit ran with NaN derivatives (with Armadillo chol() warnings) and returned a finite objective. cumul_native_multi() now sends data with non-finite predictor values or weights to polr(), as before this PR. cumul_fit_multi_c() now returns non-converged fits (convergence code 1, infinite objective) for a non-finite model matrix or offset, and cl_eval() returns Inf for a non-finite linear predictor. If the objective at the start values is infinite for another reason, cl_fit_one() now returns convergence code 1 instead of the C++ error.

Two new tests in test_cumul_native.R cover the sparse cauchit category, an infinite objective at the start values and non-finite input. Without the fixes, both tests stop with the C++ error. Tested on Triton with R 4.5.3 and testthat 3.3.2 (augdat, cumul_native and div_minimizer tests: 405 expectations, 0 failures).

[avehtari]

Thanks, I committed both suggestions. I added one commit on top.

categ_fit_multi_c() now checks the model matrix and the weight matrix before the fits. cumul_fit_multi_c() already has the same check for the model matrix. With a non-finite predictor value, the linear predictor at the start values contains Inf * 0. Not every BLAS gives NaN for this product, so the start objective can be finite. Then the new early return in ml_fit_one() does not apply. For a non-finite model matrix, categ_fit_multi_c() now returns non-converged results (convergence code 1, objective Inf, coefficients 0) and does not call ml_fit_one(). For a non-finite weight matrix, it does the same for that draw only.

Before these changes, categ_native_multi() stopped with an error for an Inf weight or predictor value: the error came from its nnet::multinom() template call. Now it returns NULL, and fit_categ() calls nnet::multinom() as before this PR.

The new test in test_categ_native.R covers categ_fit_multi_c() and categ_native_multi() with an Inf predictor value and an Inf weight. With R/divergence_minimizers.R and src/categfit.cpp from before your suggestions, it gives 4 failures and 1 error. With the new commit, test_categ_native.R passes (36 expectations, R 4.6.1, testthat 3.3.2).