[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.