[avehtari]

This PR was assisted by Claude.

Summary

This PR builds on #574 (one model matrix for all draws). For non-Gaussian GLM submodels (and Gaussian submodels with a non-identity link), divmin() still fits each projected draw separately: fit_glm_ridge_callback() calls glm_ridge(), which standardizes the shared model matrix again and runs IRLS in C++ (glm_ridge_c()). Each IRLS iteration and each line-search step calls the R function pseudo_data() (with family$linkinv(), family$mu.eta(), family$variance() and family$deviance()), and each draw-wise fit runs inside two capture.output() calls (in fit_glm_ridge_callback() for the Rcout warnings of the C++ code, and in capt_mssgs_warns() in divmin()). A profile of cv_varsel(validate_search = TRUE) for a Poisson model showed 32 % of the time in pseudo_data() and about 8 % in textConnection() and sink().

This PR makes three changes:

1. Messages without output capture. The C++ fitter returns its warning messages (as a sixth list element) instead of printing them with Rcout. glm_ridge() gets the argument print_msgs (default TRUE, so the output of glm_elnet() does not change). fit_glm_ridge_callback() uses print_msgs = FALSE and turns the messages into one warning as before. For this fitter, divmin() collects the messages and warnings with condition handlers (new capt_conds(), based on withCallingHandlers()). The other submodel fitters (lme4, mgcv, ordinal, ...) keep capt_mssgs_warns().
2. Pseudo-observations in C++. For the gaussian, binomial and poisson families with the identity, log, logit, probit, cloglog or cauchit link, the C++ code computes the pseudo-observations, the loss and the gradient itself (PseudoObsNative in src/glmfun.cpp). The link functions are transcribed from R's stats package (make.link(), and for the logit link logit_linkinv() and logit_mu_eta() in src/library/stats/src/family.c), including their thresholds. The variances and deviances are those of the stats families and of extend_family(), and the code sums the deviance in long double as R's sum() does. The IRLS loop is the existing glm_ridge() of src/glmfun.cpp; this PR only moves the pseudo-observations behind an interface (PseudoObs) with an R-callback and a native implementation. native_pseudo_codes() selects it only if family$linkinv, family$mu.eta and family$variance have the bodies of the stats family with the same name and link (cached per family and link); other families (for example Student_t() or a modified family object) keep the R callback. brms reference models use the stats family objects, so they also get the native code. The internal option projpred.glm_native = FALSE switches this off.
3. All draws in one C++ call. glm_ridge_native_multi(), which divmin() calls after the Gaussian multi-draw solve of Faster submodel projections: one model matrix and one Gaussian solve for all draws #574, standardizes the shared model matrix once and calls glm_ridge_native_multi_c(), which runs the IRLS fit of glm_ridge() for each draw from the same starting values. The subfit objects have the same elements as those of fit_glm_ridge_callback(), and the per-draw messages go to the same warn_prj_drawwise() and check_conv(). divmin() uses the draw-wise path as before for other families and links, an intercept-only submodel, unusable observation weights, and if ... contains a formal argument of glm_ridge() other than family, weights, qa_updates_max and ls_iter_max.

Results

The results are the same as on master. In the unit tests, the native pseudo-observations give bit-identical coefficients and weights for 7 family/link pairs on Linux (tested with tolerance 1e-12), and the multi-draw fits equal the draw-wise fits, also for non-converged fits. 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 1.9e-16.

Bench: rstanarm reference models, forward search, 1 CPU per run, all versions in the same job on one node (CSC Roihu, AMD EPYC 9965). Models: collinear (Poisson, N 200, P 4, from the projpred case studies), a simulated Bernoulli data set with logit and with probit link (N 400, P 30), and two models for the student data (binomial with 3 trials, N 407, P 26; Poisson, N 382, P 25). vs = varsel(nterms_max = 15) (collinear: 4), cvt50 = cv_varsel(validate_search = TRUE, nloo = 50, nterms_max = 10, refit_prj = TRUE), cvt = the same with full LOO. Times in seconds.

model	call	master	#574	this PR (with #574)	speedup
collinear (Poisson)	vs	3.4	2.2	1.3	2.6
sim. Bernoulli, logit	vs	25.7	13.3	6.3	4.1
sim. Bernoulli, probit	vs	31.0	18.3	7.8	4.0
student, binomial	vs	23.5	14.4	6.0	3.9
student, Poisson	vs	28.7	18.0	5.9	4.9
collinear (Poisson)	cvt	596.4	369.2	117.7	5.1
collinear (Poisson)	cvt50	170.6	102.9	28.5	6.0
sim. Bernoulli, logit	cvt50	728.9	479.4	169.2	4.3
sim. Bernoulli, probit	cvt50	853.2	595.3	219.0	3.9
student, binomial	cvt50	716.0	389.6	154.3	4.6
student, Poisson	cvt50	802.4	519.7	159.6	5.0

The times of cv_varsel(validate_search = FALSE) change as those of varsel().

Tests

The test suite passes with NOT_CRAN=true: 262 tests, 55 921 expectations, 0 failures, 0 errors, the same 4 warnings as master. Tested with this change applied to master 73cf0ee (projpred 2.11.0.9000) together with #574, R 4.6.1, rstanarm 2.32.2, brms 2.23.0, on CSC Roihu (Linux, AMD EPYC).

New tests in tests/testthat/test_glm_ridge.R:

• the warning messages of glm_ridge() are returned in msgs, printed only with print_msgs = TRUE, and turned into a warning by fit_glm_ridge_callback();
• the native pseudo-observations give the same fits as the R callback for 7 family/link pairs, with and without intercept and predictors, with observation weights and offsets; native_pseudo_codes() returns NULL for Gamma(), poisson("sqrt"), a modified linkinv and with the option set to FALSE;
• glm_ridge_native_multi() gives the same subfit objects and per-draw messages as the draw-wise fitter for binomial (logit, probit), poisson and gaussian (log link), with a factor predictor and observation weights, also with qa_updates_max = 2 (non-converged fits); it returns NULL for Gamma() and if offset is given;
• capt_conds() collects messages and multi-line warnings and passes errors on.

The new tests fail on master (the functions and arguments 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 the draw-wise loop of divmin(). The second PR to merge needs a rebase. The resolution keeps both changes: capt_fun(soutdmin <- do.call(sdivmin, c(list(...), dot_args))). The GLM multi-draw path of this PR returns before the removal of the observations with weight 0 that Leave the left-out observation out of the fold projections in validated LOO CV #580 adds. Leave the left-out observation out of the fold projections in validated LOO CV #580 removes these observations only for multilevel and additive submodels, and the C++ IRLS gives an observation with weight 0 a working weight of 0, as pseudo_data() does, so the GLM path needs no change.
• Keep the data-dependent bases #573: the multi-draw subfit objects of this PR need terms = shared_x$terms. Keep the data-dependent bases #573 adds this element to the draw-wise subfit objects.