From 6b44016c36d103d62dd8d056714596f787a88be5 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Thu, 22 Oct 2026 23:16:10 +0300 Subject: [PATCH 1/5] Share one model matrix across the draw-wise GLM submodel fits --- R/divergence_minimizers.R | 190 ++++++++++++++++++++++++++++++-------- 1 file changed, 151 insertions(+), 39 deletions(-) diff --git a/R/divergence_minimizers.R b/R/divergence_minimizers.R index 36345eda8..7e6a3855f 100644 --- a/R/divergence_minimizers.R +++ b/R/divergence_minimizers.R @@ -77,6 +77,15 @@ divmin <- function( rep(projpred_formulas_no_random, length(formulas)) ) } + # The draw-wise formulas of a non-multilevel non-additive submodel differ only + # in the response, so fit_glm_ridge_callback() can use one model matrix for + # all draws: + if (identical(sdivmin, fit_glm_ridge_callback) && length(formulas) > 1) { + shared_x <- shared_model_matrix(formulas, data = list(...)$data) + if (!is.null(shared_x)) { + sdivmin <- glm_ridge_callback_shared(shared_x) + } + } if (is.null(getOption("projpred.parallel_proj_trigger")) && !is.null(getOption("projpred.prll_prj_trigger"))) { @@ -173,28 +182,22 @@ divmin <- function( # submodels: fit_glm_ridge_callback <- function(formula, data, projpred_var = matrix(nrow = nrow(data)), - regul = 1e-4, thresh_conv = 1e-7, ...) { + regul = 1e-4, thresh_conv = 1e-7, + projpred_shared_x = NULL, ...) { # Preparations: - fr <- model.frame(formula, data = data, drop.unused.levels = TRUE) - da_classes <- attr(attr(fr, "terms"), "dataClasses") - nms_chr_fac <- names(da_classes)[da_classes %in% c("character", "factor")] - resp_nm <- all.vars(attr(fr, "terms"))[attr(attr(fr, "terms"), "response")] - nms_chr_fac <- setdiff(nms_chr_fac, resp_nm) - if (length(nms_chr_fac) > 0) { - xlvls <- lapply(setNames(nm = nms_chr_fac), function(nm_chr_fac) { - levels(as.factor(fr[[nm_chr_fac]])) - }) + if (!is.null(projpred_shared_x)) { + # Model matrix built once by shared_model_matrix(); only the response + # differs between the draws: + x <- projpred_shared_x$x + xlvls <- projpred_shared_x$xlvls + y <- data[[all.vars(formula[[2]])]] + names(y) <- projpred_shared_x$y_names } else { - xlvls <- NULL + xy <- glm_model_matrix(formula, data = data) + x <- xy$x + xlvls <- xy$xlvls + y <- xy$y } - # TODO: In the following model.matrix() call, allow user-specified contrasts - # to be passed to argument `contrasts.arg`. The `contrasts.arg` default - # (`NULL`) uses `options("contrasts")` internally, but it might be more - # convenient to let users specify contrasts directly. At that occasion, - # contrasts should also be tested thoroughly (not done until now). - x <- model.matrix(formula, data = fr) - x <- x[, colnames(x) != "(Intercept)", drop = FALSE] - y <- model.response(fr) # Exclude arguments from `...` which cannot be passed to glm_ridge(): dot_args <- list(...) dot_args <- dot_args[intersect( @@ -221,6 +224,72 @@ fit_glm_ridge_callback <- function(formula, data, return(sub) } +# Model frame preparations for fit_glm_ridge_callback(): predictor model matrix +# (without intercept column), levels of character and factor predictors, and +# response. +glm_model_matrix <- function(formula, data) { + fr <- model.frame(formula, data = data, drop.unused.levels = TRUE) + da_classes <- attr(attr(fr, "terms"), "dataClasses") + nms_chr_fac <- names(da_classes)[da_classes %in% c("character", "factor")] + resp_nm <- all.vars(attr(fr, "terms"))[attr(attr(fr, "terms"), "response")] + nms_chr_fac <- setdiff(nms_chr_fac, resp_nm) + if (length(nms_chr_fac) > 0) { + xlvls <- lapply(setNames(nm = nms_chr_fac), function(nm_chr_fac) { + levels(as.factor(fr[[nm_chr_fac]])) + }) + } else { + xlvls <- NULL + } + # TODO: In the following model.matrix() call, allow user-specified contrasts + # to be passed to argument `contrasts.arg`. The `contrasts.arg` default + # (`NULL`) uses `options("contrasts")` internally, but it might be more + # convenient to let users specify contrasts directly. At that occasion, + # contrasts should also be tested thoroughly (not done until now). + x <- model.matrix(formula, data = fr) + x <- x[, colnames(x) != "(Intercept)", drop = FALSE] + return(list(x = x, xlvls = xlvls, y = model.response(fr), n_fr = nrow(fr))) +} + +# Model matrix shared by the draw-wise formulas `formulas` (which have to differ +# only in the response). Returns `NULL` if sharing is not possible: different +# right-hand sides, a response that is not a plain column of `data`, or missing +# values (which would make the model frames differ between the draws). +shared_model_matrix <- function(formulas, data) { + if (!is.data.frame(data)) { + return(NULL) + } + rhs_1 <- formulas[[1]][[3]] + if (!all(vapply(formulas, function(fml) { + length(fml) == 3 && identical(fml[[3]], rhs_1) + }, logical(1)))) { + return(NULL) + } + resp_nms <- vapply(formulas, function(fml) { + if (is.name(fml[[2]])) as.character(fml[[2]]) else NA_character_ + }, character(1)) + if (anyNA(resp_nms) || !all(resp_nms %in% names(data)) || + anyNA(data[resp_nms])) { + return(NULL) + } + xy <- glm_model_matrix(formulas[[1]], data = data) + if (xy$n_fr != nrow(data)) { + return(NULL) + } + # The response names (as from model.response()) are also shared, so the + # draw-wise responses do not each get their own copy: + return(list(x = xy$x, xlvls = xy$xlvls, y_names = names(xy$y))) +} + +# fit_glm_ridge_callback() with a fixed shared model matrix. The enclosing +# environment holds only `shared_x` (relevant for the parallel case in divmin(), +# where the function is exported). +glm_ridge_callback_shared <- function(shared_x) { + force(shared_x) + function(...) { + fit_glm_ridge_callback(..., projpred_shared_x = shared_x) + } +} + # Alternative to fit_glm_ridge_callback() (may be used via global option # `projpred.glm_fitter`): fit_glm_callback <- function(formula, family, ...) { @@ -1159,6 +1228,12 @@ check_conv_s <- function(fit_s) { # Prediction functions for submodels -------------------------------------- subprd <- function(fits, newdata) { + if (!is.null(newdata) && length(fits) > 1 && inherits(fits[[1]], "subfit")) { + prd <- subprd_subfit_shared(fits, newdata) + if (!is.null(prd)) { + return(prd) + } + } prd_list <- lapply(fits, function(fit) { is_glmmPQL <- inherits(fit, "glmmPQL") is_glmm <- inherits(fit, c("lmerMod", "glmerMod")) @@ -1290,28 +1365,65 @@ predict.subfit <- function(object, newdata = NULL, ...) { if (is.null(beta)) { return(as.matrix(rep(alpha, nrow(newdata)))) } else { - # TODO: In the following model.matrix() call, allow user-specified - # contrasts to be passed to argument `contrasts.arg`. The `contrasts.arg` - # default (`NULL`) uses `options("contrasts")` internally, but it might be - # more convenient to let users specify contrasts directly. At that - # occasion, contrasts should also be tested thoroughly (not done until - # now). - x <- model.matrix(delete.response(terms(object$formula)), data = newdata, - xlev = object$xlvls) - # Need to ensure that the columns of `x` match the coefficients (note that - # `rownames(beta)` is the same as `colnames(object$x)`, see the tests): - nms_to_use <- c("(Intercept)", colnames(object$x)) - nms_to_use <- reorder_ias(nms_to_use, colnames(x)) - if (all(!is.na(nms_to_use)) && all(nms_to_use %in% colnames(x))) { - x <- x[, nms_to_use, drop = FALSE] - } else { - stop("The column names of the new model matrix don't match the ", - "column names of the original model matrix. Please notify the ", - "package maintainer.") - } - return(x %*% rbind(alpha, beta)) + return(subfit_newdata_x(object, newdata) %*% rbind(alpha, beta)) + } + } +} + +# Model matrix (with intercept column) of a `subfit` object for new data, with +# the columns in the order of the coefficients. +subfit_newdata_x <- function(object, newdata) { + # TODO: In the following model.matrix() call, allow user-specified + # contrasts to be passed to argument `contrasts.arg`. The `contrasts.arg` + # default (`NULL`) uses `options("contrasts")` internally, but it might be + # more convenient to let users specify contrasts directly. At that + # occasion, contrasts should also be tested thoroughly (not done until + # now). + x <- model.matrix(delete.response(terms(object$formula)), data = newdata, + xlev = object$xlvls) + # Need to ensure that the columns of `x` match the coefficients (note that + # `rownames(beta)` is the same as `colnames(object$x)`, see the tests): + nms_to_use <- c("(Intercept)", colnames(object$x)) + nms_to_use <- reorder_ias(nms_to_use, colnames(x)) + if (all(!is.na(nms_to_use)) && all(nms_to_use %in% colnames(x))) { + x <- x[, nms_to_use, drop = FALSE] + } else { + stop("The column names of the new model matrix don't match the ", + "column names of the original model matrix. Please notify the ", + "package maintainer.") + } + return(x) +} + +# Predictions of several `subfit` objects with the same predictor terms (the +# draw-wise fits of one submodel) for new data: one model matrix, one matrix +# product. Returns `NULL` if the fits do not share their predictor terms. +subprd_subfit_shared <- function(fits, newdata) { + fit_1 <- fits[[1]] + rhs_1 <- fit_1$formula[[3]] + nms_1 <- rownames(fit_1$beta) + is_shared <- all(vapply(fits, function(fit) { + inherits(fit, "subfit") && identical(fit$formula[[3]], rhs_1) && + identical(rownames(fit$beta), nms_1) && identical(fit$xlvls, fit_1$xlvls) + }, logical(1))) + if (!is_shared) { + return(NULL) + } + alphas <- vapply(fits, function(fit) as.numeric(fit$alpha), numeric(1)) + if (is.null(fit_1$beta) || NROW(fit_1$beta) == 0) { + if (is.null(fit_1$beta)) { + # As predict.subfit() for `beta = NULL`: + return(matrix(alphas, nrow = nrow(newdata), ncol = length(fits), + byrow = TRUE)) } + return(NULL) + } + betas <- vapply(fits, function(fit) as.numeric(fit$beta), + numeric(length(nms_1))) + if (!is.matrix(betas)) { + betas <- matrix(betas, nrow = 1) } + return(subfit_newdata_x(fit_1, newdata) %*% rbind(alphas, betas)) } #' @noRd From 3325739fecc23fcd8a994bfe4346dc36ba494ff5 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sat, 24 Oct 2026 11:23:15 +0300 Subject: [PATCH 2/5] Solve Gaussian identity-link submodel projections for all draws at once --- NEWS.md | 2 + R/divergence_minimizers.R | 71 +++++++++++- tests/testthat/test_submodel_fit_shared.R | 125 ++++++++++++++++++++++ 3 files changed, 196 insertions(+), 2 deletions(-) create mode 100644 tests/testthat/test_submodel_fit_shared.R diff --git a/NEWS.md b/NEWS.md index 94ffe21cd..9ae4aed91 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,6 +7,8 @@ If you read this from a place other than ) + ## Minor changes ## Bug fixes and maintenance diff --git a/R/divergence_minimizers.R b/R/divergence_minimizers.R index 7e6a3855f..eab961494 100644 --- a/R/divergence_minimizers.R +++ b/R/divergence_minimizers.R @@ -79,10 +79,18 @@ divmin <- function( } # The draw-wise formulas of a non-multilevel non-additive submodel differ only # in the response, so fit_glm_ridge_callback() can use one model matrix for - # all draws: - if (identical(sdivmin, fit_glm_ridge_callback) && length(formulas) > 1) { + # all draws (and the joint Gaussian solve can fit them, also for a single + # draw): + if (identical(sdivmin, fit_glm_ridge_callback) && length(formulas) >= 1) { shared_x <- shared_model_matrix(formulas, data = list(...)$data) if (!is.null(shared_x)) { + # Gaussian submodels with the identity link: one solve for all draws + # (this fitter cannot throw messages or warnings in this case): + outdmin <- glm_ridge_gauss_multi(formulas, shared_x = shared_x, + dot_args = list(...)) + if (!is.null(outdmin)) { + return(outdmin) + } sdivmin <- glm_ridge_callback_shared(shared_x) } } @@ -290,6 +298,65 @@ glm_ridge_callback_shared <- function(shared_x) { } } +# Fits all draw-wise submodels of fit_glm_ridge_callback() at once, for a +# Gaussian submodel with the identity link. For this family, glm_ridge() solves +# one weighted ridge least-squares problem in one exact step. The draws share +# the model matrix, the observation weights, the standardization and the +# penalty. Thus, one factorization with one right-hand side per draw gives the +# same solutions. Returns `NULL` if one of the checks below fails; then divmin() +# uses the draw-wise path. +glm_ridge_gauss_multi <- function(formulas, shared_x, dot_args) { + family <- dot_args$family + if (is.null(family) || family$family != "gaussian" || + family$link != "identity") { + return(NULL) + } + # fit_glm_ridge_callback() passes only the formal arguments of glm_ridge() on + # and ignores the others. Of these, this function uses only `family`, + # `weights` and `thresh` (not used for this family): + args_glm_ridge <- intersect(names(dot_args), methods::formalArgs(glm_ridge)) + if (!all(args_glm_ridge %in% c("family", "weights", "thresh"))) { + return(NULL) + } + x <- shared_x$x + if (ncol(x) == 0) { + return(NULL) + } + data <- dot_args$data + weights <- dot_args$weights %||% rep(1.0, nrow(x)) + regul <- dot_args$regul %||% 1e-4 + if (length(weights) != nrow(x) || anyNA(weights)) { + return(NULL) + } + resp_nms <- vapply(formulas, function(fml) as.character(fml[[2]]), + character(1)) + y <- as.matrix(data[resp_nms]) + # As in glm_ridge(): + transf <- standardization(x, center = TRUE, scale = TRUE, weights = weights) + if (any(transf$scale == 0)) { + return(NULL) + } + x_std <- cbind(1, t((t(x) - transf$shift) / transf$scale)) + sqrt_w <- sqrt(weights) + x_w <- x_std * sqrt_w + regmat <- diag(c(0, rep(regul, ncol(x))), nrow = ncol(x_std)) + coefs <- solve(crossprod(x_w) + regmat, crossprod(x_w, y * sqrt_w)) + beta <- coefs[-1, , drop = FALSE] / transf$scale + beta0 <- coefs[1, ] - colSums(transf$shift * beta) + # As from glm_ridge_c() (double, also for integer observation weights): + w <- matrix(as.numeric(weights), ncol = 1) + return(lapply(seq_along(formulas), function(s) { + beta_s <- beta[, s, drop = FALSE] + dimnames(beta_s) <- list(colnames(x), NULL) + y_s <- y[, s] + names(y_s) <- shared_x$y_names + structure(list(alpha = unname(beta0[s]), beta = beta_s, w = w, + formula = formulas[[s]], x = x, y = y_s, + xlvls = shared_x$xlvls), + class = "subfit") + })) +} + # Alternative to fit_glm_ridge_callback() (may be used via global option # `projpred.glm_fitter`): fit_glm_callback <- function(formula, family, ...) { diff --git a/tests/testthat/test_submodel_fit_shared.R b/tests/testthat/test_submodel_fit_shared.R new file mode 100644 index 000000000..d5560f0d8 --- /dev/null +++ b/tests/testthat/test_submodel_fit_shared.R @@ -0,0 +1,125 @@ +context("shared model matrix and joint Gaussian solve for submodel fits") + +# Data with one response column per projected draw, as divmin() gets them, and +# the draw-wise formulas. The formulas include a factor and an interaction. +shared_fit_data <- function(n = 40, S = 5, seed = 4721) { + set.seed(seed) + dat <- data.frame(x1 = rnorm(n), x2 = rnorm(n), + f = factor(sample(c("a", "b", "c"), n, replace = TRUE))) + for (s in seq_len(S)) { + dat[[paste0("y", s)]] <- 0.5 * dat$x1 - 0.3 * dat$x2 + (dat$f == "b") + + rnorm(n, sd = 0.3) + } + fmls <- lapply(seq_len(S), function(s) { + as.formula(paste0("y", s, " ~ x1 + x2 + f + x1:x2")) + }) + list(data = dat, formulas = fmls, w = runif(n, 0.5, 2)) +} + +test_that(paste( + "the joint Gaussian solve gives the same fits as the draw-wise", + "fit_glm_ridge_callback()" +), { + d <- shared_fit_data() + shared_x <- shared_model_matrix(d$formulas, data = d$data) + expect_false(is.null(shared_x)) + fits_multi <- glm_ridge_gauss_multi( + d$formulas, shared_x = shared_x, + dot_args = list(data = d$data, family = gaussian(), weights = d$w) + ) + expect_length(fits_multi, length(d$formulas)) + fits_loop <- lapply(d$formulas, fit_glm_ridge_callback, data = d$data, + family = gaussian(), weights = d$w) + for (s in seq_along(d$formulas)) { + expect_s3_class(fits_multi[[s]], "subfit") + expect_identical(names(fits_multi[[s]]), names(fits_loop[[s]])) + expect_equal(unname(fits_multi[[s]]$alpha), unname(fits_loop[[s]]$alpha), + tolerance = 1e-10) + expect_equal(fits_multi[[s]]$beta, fits_loop[[s]]$beta, tolerance = 1e-10) + expect_identical(fits_multi[[s]]$x, fits_loop[[s]]$x) + expect_identical(fits_multi[[s]]$y, fits_loop[[s]]$y) + expect_identical(fits_multi[[s]]$xlvls, fits_loop[[s]]$xlvls) + expect_equal(fits_multi[[s]]$w, fits_loop[[s]]$w) + } +}) + +test_that("the joint Gaussian solve falls back for other families", { + d <- shared_fit_data() + shared_x <- shared_model_matrix(d$formulas, data = d$data) + expect_null(glm_ridge_gauss_multi( + d$formulas, shared_x = shared_x, + dot_args = list(data = d$data, family = gaussian(link = "log")) + )) + expect_null(glm_ridge_gauss_multi( + d$formulas, shared_x = shared_x, + dot_args = list(data = d$data, family = poisson()) + )) +}) + +test_that("shared_model_matrix() returns NULL if the model matrix differs", { + d <- shared_fit_data() + # Different right-hand sides: + fmls_rhs <- d$formulas + fmls_rhs[[2]] <- y2 ~ x1 + x2 + expect_null(shared_model_matrix(fmls_rhs, data = d$data)) + # A response that is not a plain column of the data: + fmls_resp <- d$formulas + fmls_resp[[1]] <- log(abs(y1)) ~ x1 + x2 + f + x1:x2 + expect_null(shared_model_matrix(fmls_resp, data = d$data)) + # A missing predictor value: + dat_na <- d$data + dat_na$x2[3] <- NA + expect_null(shared_model_matrix(d$formulas, data = dat_na)) +}) + +test_that(paste( + "the joint predictions for new data equal the draw-wise predictions of", + "predict.subfit()" +), { + d <- shared_fit_data() + shared_x <- shared_model_matrix(d$formulas, data = d$data) + fits <- glm_ridge_gauss_multi( + d$formulas, shared_x = shared_x, + dot_args = list(data = d$data, family = gaussian(), weights = d$w) + ) + nd <- d$data[c(2, 5, 7, 11), ] + prd_loop <- do.call(cbind, lapply(fits, predict, newdata = nd)) + prd_shared <- subprd_subfit_shared(fits, nd) + expect_equal(unname(prd_shared), unname(prd_loop), tolerance = 1e-12) + expect_equal(unname(subprd(fits, nd)), unname(prd_loop), tolerance = 1e-12) +}) + +test_that("project() uses the joint Gaussian solve for a Gaussian submodel", { + d <- shared_fit_data(n = 60) + dat <- d$data[, c("x1", "x2", "f")] + dat$y <- d$data$y1 + fml <- y ~ x1 + x2 + f + mm <- model.matrix(fml, dat) + ndraws <- 100 + beta <- drop(solve(crossprod(mm), crossprod(mm, dat$y))) + + matrix(rnorm(ncol(mm) * ndraws, sd = 0.05), ncol(mm)) + fit <- list(beta = beta, data = dat) + ref_predfun <- function(fit, newdata = NULL) { + nd <- if (is.null(newdata)) fit$data else newdata + model.matrix(~ x1 + x2 + f, nd) %*% fit$beta + } + refm <- suppressWarnings(init_refmodel( + fit, data = dat, formula = fml, family = gaussian(), + ref_predfun = ref_predfun, dis = rep(0.3, ndraws) + )) + calls <- new.env() + calls$n <- 0L + ns <- environment(glm_ridge_gauss_multi) + suppressMessages(trace( + "glm_ridge_gauss_multi", where = ns, print = FALSE, + exit = bquote(if (!is.null(returnValue(NULL))) { + assign("n", get("n", envir = .(calls)) + 1L, envir = .(calls)) + }) + )) + on.exit(suppressMessages(untrace("glm_ridge_gauss_multi", where = ns)), + add = TRUE) + prj <- project(refm, predictor_terms = c("x1", "f"), nclusters = 5, + seed = 3906, verbose = FALSE) + expect_gt(calls$n, 0L) + expect_length(prj$outdmin, 5) +}) From 3eebda0c2f7ffc19036bbde801fbb3a7675fe7c9 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sat, 24 Oct 2026 13:11:26 +0300 Subject: [PATCH 3/5] Add the PR number to the NEWS entry --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 9ae4aed91..96fc28373 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,7 +7,7 @@ If you read this from a place other than ) +* 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) ## Minor changes From 9fabbc5640196760a33ec68119334bcc5c28c1a7 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Mon, 26 Oct 2026 10:36:04 +0200 Subject: [PATCH 4/5] Faster GLM submodel fits: C++ pseudo-observations, all draws in one C++ call and no output capture --- NEWS.md | 2 + R/RcppExports.R | 8 + R/divergence_minimizers.R | 127 ++++++++++++-- R/glmfun.R | 99 +++++++++-- R/misc.R | 20 +++ src/RcppExports.cpp | 52 ++++++ src/glmfun.cpp | 283 +++++++++++++++++++++++++++++--- tests/testthat/test_glm_ridge.R | 136 +++++++++++++++ 8 files changed, 678 insertions(+), 49 deletions(-) diff --git a/NEWS.md b/NEWS.md index 96fc28373..92da0a6d0 100644 --- a/NEWS.md +++ b/NEWS.md @@ -9,6 +9,8 @@ If you read this from a place other than 0) { warning(paste(out_capt, collapse = "\n")) } @@ -357,6 +376,84 @@ glm_ridge_gauss_multi <- function(formulas, shared_x, dot_args) { })) } +# Fits all draw-wise submodels of fit_glm_ridge_callback() in one C++ call, for +# the families and links with native pseudo-observations (see +# native_pseudo_codes()). The draws share the model matrix, the observation +# weights, the standardization and the penalty, so this sets up glm_ridge() +# once. Each draw gets the IRLS iterations of glm_ridge() and the same result. +# Attribute `mssgs_warns_capts` holds the messages of each draw, as divmin() +# collects them from fit_glm_ridge_callback(). Returns `NULL` for other cases; +# divmin() then fits the draws one by one. +glm_ridge_native_multi <- function(formulas, shared_x, dot_args) { + family <- dot_args$family + if (is.null(family)) { + return(NULL) + } + family <- extend_family(family) + codes <- native_pseudo_codes(family) + if (is.null(codes)) { + return(NULL) + } + # fit_glm_ridge_callback() passes only the formal arguments of glm_ridge() on. + # Of these, the following are reproduced here: + args_glm_ridge <- intersect(names(dot_args), methods::formalArgs(glm_ridge)) + if (!all(args_glm_ridge %in% c("family", "weights", "qa_updates_max", + "ls_iter_max"))) { + return(NULL) + } + x <- shared_x$x + if (ncol(x) == 0) { + return(NULL) + } + data <- dot_args$data + weights <- dot_args$weights %||% rep(1.0, nrow(x)) + if (length(weights) != nrow(x) || anyNA(weights)) { + return(NULL) + } + regul <- dot_args$regul %||% 1e-4 + thresh <- dot_args$thresh_conv %||% 1e-7 + # As in glm_ridge(): + if (family$family == "gaussian" && family$link == "identity") { + qa_updates_max <- 1 + ls_iter_max <- 1 + } else { + qa_updates_max <- dot_args$qa_updates_max %||% 100 + ls_iter_max <- dot_args$ls_iter_max %||% 30 + } + resp_nms <- vapply(formulas, function(fml) as.character(fml[[2]]), + character(1)) + y <- as.matrix(data[resp_nms]) + transf <- standardization(x, center = TRUE, scale = TRUE, weights = weights) + penalty <- rep(1.0, ncol(x)) + penalty[transf$scale == 0] <- Inf + transf$scale[transf$scale == 0] <- 1 + x_std <- t((t(x) - transf$shift) / transf$scale) + out <- glm_ridge_native_multi_c( + x_std, y, as.numeric(weights), rep(0.0, nrow(x)), codes[1], codes[2], + regul, TRUE, penalty, rep(0, ncol(x) + 1), weights, thresh, + qa_updates_max, ls_iter_max + ) + beta <- out$beta / transf$scale + beta0 <- as.vector(out$beta0) - colSums(transf$shift * beta) + outdmin <- lapply(seq_along(formulas), function(s) { + beta_s <- beta[, s, drop = FALSE] + dimnames(beta_s) <- list(colnames(x), NULL) + y_s <- y[, s] + names(y_s) <- shared_x$y_names + structure(list(alpha = beta0[s], beta = beta_s, + w = out$w[, s, drop = FALSE], formula = formulas[[s]], + x = x, y = y_s, xlvls = shared_x$xlvls), + class = "subfit") + }) + # As fit_glm_ridge_callback() turns the messages into one warning and + # capt_conds() splits it into lines: + attr(outdmin, "mssgs_warns_capts") <- lapply(out$msgs, function(msgs_s) { + msgs_s <- unique(grep("[Ww]arning|bug", msgs_s, value = TRUE)) + setdiff(as.character(unlist(strsplit(msgs_s, "\n", fixed = TRUE))), "") + }) + return(outdmin) +} + # Alternative to fit_glm_ridge_callback() (may be used via global option # `projpred.glm_fitter`): fit_glm_callback <- function(formula, family, ...) { diff --git a/R/glmfun.R b/R/glmfun.R index 044678b9a..dcfe0a47a 100644 --- a/R/glmfun.R +++ b/R/glmfun.R @@ -207,10 +207,12 @@ glm_ridge <- function(x, y, family = gaussian(), lambda = 0, thresh = 1e-7, qa_updates_max = NULL, weights = NULL, offset = NULL, obsvar = 0, intercept = TRUE, penalty = NULL, normalize = TRUE, beta_init = NULL, beta0_init = NULL, - ls_iter_max = 30) { + ls_iter_max = 30, print_msgs = TRUE) { # # Fits GLM with ridge penalty on the regression coefficients. # Does not handle any dispersion parameters. + # The warning messages of the C++ fitter are returned in element `msgs`. With + # `print_msgs = TRUE`, they are also printed to the standard output. # if (is.null(x)) { x <- matrix(ncol = 0, nrow = length(y)) @@ -251,14 +253,13 @@ glm_ridge <- function(x, y, family = gaussian(), lambda = 0, thresh = 1e-7, # constant predictor) x <- matrix(rep(1, length(y)), ncol = 1) w0 <- weights - pseudo_obs <- function(f, wprev) - pseudo_data(f, y, family, offset = offset, weights = weights, - obsvar = obsvar, wprev = wprev) - out <- glm_ridge_c(x, pseudo_obs, lambda, FALSE, 0, beta_start, w0, - thresh, qa_updates_max, ls_iter_max) + out <- glm_ridge_fit_c(x, y, family, weights, offset, obsvar, lambda, + FALSE, 0, beta_start, w0, thresh, qa_updates_max, + ls_iter_max) + print_ridge_msgs(out[[6]], print_msgs) return(list(beta = matrix(integer(length = 0)), beta0 = as.vector(out[[1]]), w = out[[3]], loss = out[[4]], - qa_updates = out[[5]])) + qa_updates = out[[5]], msgs = out[[6]])) } else { # null model with no predictors and no intercept pseudo_obs <- function(f, wprev) @@ -266,7 +267,7 @@ glm_ridge <- function(x, y, family = gaussian(), lambda = 0, thresh = 1e-7, obsvar = obsvar, wprev = wprev) pobs <- pseudo_obs(rep(0, length(y)), weights) return(list(beta = matrix(integer(length = 0)), beta0 = 0, w = pobs$wobs, - qa_updates = 0)) + qa_updates = 0, msgs = character())) } } @@ -284,21 +285,91 @@ glm_ridge <- function(x, y, family = gaussian(), lambda = 0, thresh = 1e-7, # compute the solution w0 <- weights - pseudo_obs <- function(f, wprev) - pseudo_data(f, y, family, offset = offset, weights = weights, - obsvar = obsvar, wprev = wprev) - out <- glm_ridge_c(x, pseudo_obs, lambda, intercept, penalty, beta_start, w0, - thresh, qa_updates_max, ls_iter_max) + out <- glm_ridge_fit_c(x, y, family, weights, offset, obsvar, lambda, + intercept, penalty, beta_start, w0, thresh, + qa_updates_max, ls_iter_max) beta <- out[[1]] beta0 <- as.vector(out[[2]]) w <- out[[3]] loss <- out[[4]] + print_ridge_msgs(out[[6]], print_msgs) # return the intercept and the coefficients on the original scale beta_orig <- beta / transf$scale beta0_orig <- beta0 - sum(transf$shift * beta_orig) out <- nlist(beta = beta_orig, beta0 = beta0_orig, w, - qa_updates = out[[5]]) + qa_updates = out[[5]], msgs = out[[6]]) return(out) } + +# Calls the C++ fitter: glm_ridge_native_c() if it supports the family and the +# link (see native_pseudo_codes()), otherwise glm_ridge_c() with pseudo_data() +# as the R callback. Both give the same results. +glm_ridge_fit_c <- function(x, y, family, weights, offset, obsvar, lambda, + intercept, penalty, beta_start, w0, thresh, + qa_updates_max, ls_iter_max) { + codes <- native_pseudo_codes(family) + if (!is.null(codes)) { + return(glm_ridge_native_c( + x, as.numeric(y), as.numeric(weights), as.numeric(offset), codes[1], + codes[2], lambda, intercept, penalty, beta_start, w0, thresh, + qa_updates_max, ls_iter_max + )) + } + pseudo_obs <- function(f, wprev) { + pseudo_data(f, y, family, offset = offset, weights = weights, + obsvar = obsvar, wprev = wprev) + } + return(glm_ridge_c(x, pseudo_obs, lambda, intercept, penalty, beta_start, + w0, thresh, qa_updates_max, ls_iter_max)) +} + +# Family and link codes for glm_ridge_native_c(), or `NULL` if the +# pseudo-observations have to be computed by pseudo_data() in R. The native +# code reproduces the family objects of the stats package, so other family +# objects with the same names are not covered. Global option +# `projpred.glm_native` (default `TRUE`) switches the native code off. +native_pseudo_codes <- function(family) { + if (!getOption("projpred.glm_native", TRUE)) { + return(NULL) + } + fam_code <- match(family$family, c("gaussian", "binomial", "poisson")) + link_code <- match(family$link, c("identity", "log", "logit", "probit", + "cloglog", "cauchit")) + if (is.na(fam_code) || is.na(link_code)) { + return(NULL) + } + # The function bodies of the stats family, cached per family and link: + key <- paste(family$family, family$link) + ref_bodies <- native_pseudo_cache[[key]] + if (is.null(ref_bodies)) { + stats_fam <- get(family$family, envir = asNamespace("stats"), + mode = "function") + stats_fam <- tryCatch(stats_fam(link = family$link), + error = function(e) NULL) + if (is.null(stats_fam)) { + ref_bodies <- FALSE + } else { + ref_bodies <- lapply(stats_fam[c("linkinv", "mu.eta", "variance")], body) + } + assign(key, ref_bodies, envir = native_pseudo_cache) + } + if (isFALSE(ref_bodies) || + !identical(body(family$linkinv), ref_bodies$linkinv) || + !identical(body(family$mu.eta), ref_bodies$mu.eta) || + !identical(body(family$variance), ref_bodies$variance)) { + return(NULL) + } + return(c(fam_code, link_code)) +} +native_pseudo_cache <- new.env(parent = emptyenv()) + +# Prints the warning messages of glm_ridge_c() to the standard output, one line +# per message: +print_ridge_msgs <- function(msgs, print_msgs) { + if (print_msgs && length(msgs) > 0) { + cat(paste0(msgs, "\n"), sep = "") + } + return(invisible()) +} diff --git a/R/misc.R b/R/misc.R index 28ee8b9c9..014b05093 100644 --- a/R/misc.R +++ b/R/misc.R @@ -690,6 +690,26 @@ capt_mssgs_warns <- function(expr) { return(mssgs_warns_capt) } +# Like capt_mssgs_warns(), but collects the messages and warnings by condition +# handlers. Faster, but does not capture text that is written directly to the +# message stream, so use it only for fitters that throw R conditions. +capt_conds <- function(expr) { + capt <- character() + withCallingHandlers( + expr, + warning = function(w) { + capt <<- c(capt, conditionMessage(w)) + invokeRestart("muffleWarning") + }, + message = function(m) { + capt <<- c(capt, sub("\n$", "", conditionMessage(m))) + invokeRestart("muffleMessage") + } + ) + capt <- as.character(unlist(strsplit(capt, "\n", fixed = TRUE))) + return(setdiff(capt, "")) +} + # Parse the argument containing the observation weights (`wobs` or `weights`) # for the $ppd() functions used by proj_predict(): parse_wobs_ppd <- function(wobs, n_obs) { diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index c85e42fe8..26bc3a105 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -56,10 +56,62 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// glm_ridge_native_c +List glm_ridge_native_c(arma::mat x, arma::vec y, arma::vec weights, arma::vec offset, int family_code, int link_code, double lambda, bool intercept, arma::vec penalty, arma::vec beta_init, arma::vec w_init, double thresh, int qa_updates_max, int ls_iter_max, bool debug); +RcppExport SEXP _projpred_glm_ridge_native_c(SEXP xSEXP, SEXP ySEXP, SEXP weightsSEXP, SEXP offsetSEXP, SEXP family_codeSEXP, SEXP link_codeSEXP, SEXP lambdaSEXP, SEXP interceptSEXP, SEXP penaltySEXP, SEXP beta_initSEXP, SEXP w_initSEXP, SEXP threshSEXP, SEXP qa_updates_maxSEXP, SEXP ls_iter_maxSEXP, SEXP debugSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::mat >::type x(xSEXP); + Rcpp::traits::input_parameter< arma::vec >::type y(ySEXP); + Rcpp::traits::input_parameter< arma::vec >::type weights(weightsSEXP); + Rcpp::traits::input_parameter< arma::vec >::type offset(offsetSEXP); + Rcpp::traits::input_parameter< int >::type family_code(family_codeSEXP); + Rcpp::traits::input_parameter< int >::type link_code(link_codeSEXP); + Rcpp::traits::input_parameter< double >::type lambda(lambdaSEXP); + Rcpp::traits::input_parameter< bool >::type intercept(interceptSEXP); + Rcpp::traits::input_parameter< arma::vec >::type penalty(penaltySEXP); + Rcpp::traits::input_parameter< arma::vec >::type beta_init(beta_initSEXP); + Rcpp::traits::input_parameter< arma::vec >::type w_init(w_initSEXP); + Rcpp::traits::input_parameter< double >::type thresh(threshSEXP); + Rcpp::traits::input_parameter< int >::type qa_updates_max(qa_updates_maxSEXP); + Rcpp::traits::input_parameter< int >::type ls_iter_max(ls_iter_maxSEXP); + Rcpp::traits::input_parameter< bool >::type debug(debugSEXP); + rcpp_result_gen = Rcpp::wrap(glm_ridge_native_c(x, y, weights, offset, family_code, link_code, lambda, intercept, penalty, beta_init, w_init, thresh, qa_updates_max, ls_iter_max, debug)); + return rcpp_result_gen; +END_RCPP +} +// glm_ridge_native_multi_c +List glm_ridge_native_multi_c(arma::mat x, arma::mat y, arma::vec weights, arma::vec offset, int family_code, int link_code, double lambda, bool intercept, arma::vec penalty, arma::vec beta_init, arma::vec w_init, double thresh, int qa_updates_max, int ls_iter_max, bool debug); +RcppExport SEXP _projpred_glm_ridge_native_multi_c(SEXP xSEXP, SEXP ySEXP, SEXP weightsSEXP, SEXP offsetSEXP, SEXP family_codeSEXP, SEXP link_codeSEXP, SEXP lambdaSEXP, SEXP interceptSEXP, SEXP penaltySEXP, SEXP beta_initSEXP, SEXP w_initSEXP, SEXP threshSEXP, SEXP qa_updates_maxSEXP, SEXP ls_iter_maxSEXP, SEXP debugSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::mat >::type x(xSEXP); + Rcpp::traits::input_parameter< arma::mat >::type y(ySEXP); + Rcpp::traits::input_parameter< arma::vec >::type weights(weightsSEXP); + Rcpp::traits::input_parameter< arma::vec >::type offset(offsetSEXP); + Rcpp::traits::input_parameter< int >::type family_code(family_codeSEXP); + Rcpp::traits::input_parameter< int >::type link_code(link_codeSEXP); + Rcpp::traits::input_parameter< double >::type lambda(lambdaSEXP); + Rcpp::traits::input_parameter< bool >::type intercept(interceptSEXP); + Rcpp::traits::input_parameter< arma::vec >::type penalty(penaltySEXP); + Rcpp::traits::input_parameter< arma::vec >::type beta_init(beta_initSEXP); + Rcpp::traits::input_parameter< arma::vec >::type w_init(w_initSEXP); + Rcpp::traits::input_parameter< double >::type thresh(threshSEXP); + Rcpp::traits::input_parameter< int >::type qa_updates_max(qa_updates_maxSEXP); + Rcpp::traits::input_parameter< int >::type ls_iter_max(ls_iter_maxSEXP); + Rcpp::traits::input_parameter< bool >::type debug(debugSEXP); + rcpp_result_gen = Rcpp::wrap(glm_ridge_native_multi_c(x, y, weights, offset, family_code, link_code, lambda, intercept, penalty, beta_init, w_init, thresh, qa_updates_max, ls_iter_max, debug)); + return rcpp_result_gen; +END_RCPP +} static const R_CallMethodDef CallEntries[] = { {"_projpred_glm_elnet_c", (DL_FUNC) &_projpred_glm_elnet_c, 14}, {"_projpred_glm_ridge_c", (DL_FUNC) &_projpred_glm_ridge_c, 11}, + {"_projpred_glm_ridge_native_c", (DL_FUNC) &_projpred_glm_ridge_native_c, 15}, + {"_projpred_glm_ridge_native_multi_c", (DL_FUNC) &_projpred_glm_ridge_native_multi_c, 15}, {NULL, NULL, 0} }; diff --git a/src/glmfun.cpp b/src/glmfun.cpp index 2d2c72385..8491666e6 100644 --- a/src/glmfun.cpp +++ b/src/glmfun.cpp @@ -3,6 +3,8 @@ #include #include #include +#include +#include //[[Rcpp::depends(RcppArmadillo)]] @@ -291,6 +293,141 @@ List glm_elnet_c(arma::mat x, // input matrix +/** + * Pseudo-observations of the quadratic approximation to the loss function at + * the latent values f (see pseudo_data() in R): locations z, weights w, the + * pointwise gradient of the loss w.r.t. f, and the loss. + */ +struct PseudoObs { + virtual ~PseudoObs() {} + virtual void eval(const vec& f, const vec& wprev, vec& z, vec& w, vec& grad, + double& loss) = 0; +}; + +// Pseudo-observations from an R function (any family supported by +// pseudo_data()). +struct PseudoObsR : public PseudoObs { + Function fn; + PseudoObsR(Function fn_) : fn(fn_) {} + void eval(const vec& f, const vec& wprev, vec& z, vec& w, vec& grad, + double& loss) { + List obs = fn(f, wprev); + z = as(obs["z"]); + w = as(obs["wobs"]); + grad = as(obs["grad"]); + loss = (double) obs["loss"]; + } +}; + +// Family and link codes of PseudoObsNative (see native_pseudo_codes() in R). +enum { FAM_GAUSSIAN = 1, FAM_BINOMIAL = 2, FAM_POISSON = 3 }; +enum { LINK_IDENTITY = 1, LINK_LOG = 2, LINK_LOGIT = 3, LINK_PROBIT = 4, + LINK_CLOGLOG = 5, LINK_CAUCHIT = 6 }; + +// pmax() and pmin() of R for a scalar bound (NaN propagates): +static inline double r_pmax(double x, double bound) { + return (std::isnan(x) || x >= bound) ? x : bound; +} +static inline double r_pmin(double x, double bound) { + return (std::isnan(x) || x <= bound) ? x : bound; +} + +// Pseudo-observations computed in C++. linkinv() and mu.eta() transcribe the +// link functions of R's stats package (make.link(), and for the logit link +// logit_linkinv() and logit_mu_eta() in src/library/stats/src/family.c of R, +// GPL-2 | GPL-3). The variances and deviances are those of the stats families +// and of extend_family(): +struct PseudoObsNative : public PseudoObs { + vec y, weights, offset; + int fam, link; + PseudoObsNative(const vec& y_, const vec& weights_, const vec& offset_, + int fam_, int link_) : + y(y_), weights(weights_), offset(offset_), fam(fam_), link(link_) {} + + // linkinv() and mu.eta() of the link: + void link_funs(double eta, double& mu, double& dmu) const { + const double eps = std::numeric_limits::epsilon(); + switch (link) { + case LINK_IDENTITY: + mu = eta; + dmu = 1.0; + break; + case LINK_LOG: + mu = r_pmax(std::exp(eta), eps); + dmu = mu; + break; + case LINK_LOGIT: { + // As logit_linkinv() and logit_mu_eta() in R's stats/src/family.c: + const double thr = 30.0; + double tmp = (eta < -thr) ? eps : ((eta > thr) ? 1.0 / eps : std::exp(eta)); + mu = tmp / (1.0 + tmp); + double opexp = 1.0 + std::exp(eta); + dmu = (eta > thr || eta < -thr) ? eps : std::exp(eta) / (opexp * opexp); + break; + } + case LINK_PROBIT: { + const double thr = -R::qnorm(eps, 0.0, 1.0, 1, 0); + double eta_c = r_pmin(r_pmax(eta, -thr), thr); + mu = R::pnorm(eta_c, 0.0, 1.0, 1, 0); + dmu = r_pmax(R::dnorm(eta, 0.0, 1.0, 0), eps); + break; + } + case LINK_CLOGLOG: { + mu = r_pmax(r_pmin(-std::expm1(-std::exp(eta)), 1.0 - eps), eps); + double eta_c = r_pmin(eta, 700.0); + dmu = r_pmax(std::exp(eta_c) * std::exp(-std::exp(eta_c)), eps); + break; + } + case LINK_CAUCHIT: { + const double thr = -R::qcauchy(eps, 0.0, 1.0, 1, 0); + double eta_c = r_pmin(r_pmax(eta, -thr), thr); + mu = R::pcauchy(eta_c, 0.0, 1.0, 1, 0); + dmu = r_pmax(R::dcauchy(eta, 0.0, 1.0, 0), eps); + break; + } + default: + Rcpp::stop("Unknown link code."); + } + } + + void eval(const vec& f, const vec& wprev, vec& z, vec& w, vec& grad, + double& loss) { + size_t n = f.n_elem; + z.set_size(n); + w.set_size(n); + grad.set_size(n); + // R's sum() accumulates in long double: + long double dev_sum = 0.0; + double mu, dmu, var, dev; + for (size_t i = 0; i < n; ++i) { + link_funs(f(i) + offset(i), mu, dmu); + z(i) = f(i) + (y(i) - mu) / dmu; + switch (fam) { + case FAM_GAUSSIAN: + var = 1.0; + dev = -2.0 * weights(i) * (-0.5 * ((y(i) - mu) * (y(i) - mu)) - 0.0); + break; + case FAM_BINOMIAL: + var = mu * (1.0 - mu); + dev = -2.0 * (weights(i) * (y(i) * std::log(mu) + + (1.0 - y(i)) * std::log(1.0 - mu))); + break; + case FAM_POISSON: + var = mu; + dev = -2.0 * (weights(i) * (y(i) * std::log(mu) - mu)); + break; + default: + Rcpp::stop("Unknown family code."); + } + w(i) = (weights(i) * (dmu * dmu)) / var; + grad(i) = (-w(i)) * (z(i) - f(i)); + dev_sum += dev; + } + loss = 0.5 * (double) dev_sum; + } +}; + + /** Internal function that gives the output in c++ way (writes into allocated memory). * See glm_ridge_c for a wrapper that is being called from R. */ @@ -298,8 +435,9 @@ void glm_ridge( vec& beta, // output: regression coefficients (contains int double& loss, // output: value of the loss function vec& w, // output: weights of the pseudo-gaussian observations at the optimum (needed for Student-t model) int& qau, // output: number of quadratic approximation updates + std::vector& msgs, // output: warning messages arma::mat x, - Function pseudo_obs, + PseudoObs& pseudo_obs, double lambda, bool intercept, arma::vec penalty, // relative penalties for the variables @@ -334,11 +472,11 @@ void glm_ridge( vec& beta, // output: regression coefficients (contains int mat regmat = lambda*diagmat(penalty);//eye(D,D); // regularization matrix // initial quadratic approximation - List obs = pseudo_obs(f,w); - vec z = as(obs["z"]); - w = as(obs["wobs"]); - grad_f = as(obs["grad"]); - double loss_initial = ((double) obs["loss"]) + elnet_penalty(beta, lambda, 0, penalty); + vec z(n); + double loss_obs; + pseudo_obs.eval(f, w, z, w, grad_f, loss_obs); + vec z_ls(n), w_ls(n), grad_ls(n); // pseudo-observations in the line search + double loss_initial = loss_obs + elnet_penalty(beta, lambda, 0, penalty); double loss_old = loss_initial; // will be updated iteratively loss = loss_initial; // will be updated iteratively double tol = thresh*fabs(loss_initial); // threshold for convergence @@ -370,8 +508,8 @@ void glm_ridge( vec& beta, // output: regression coefficients (contains int t = b*t; f = x*(beta+t*dbeta); - obs = pseudo_obs(f,w); - loss = ((double) obs["loss"]) + elnet_penalty(beta+t*dbeta, lambda, 0, penalty); + pseudo_obs.eval(f, w, z_ls, w_ls, grad_ls, loss_obs); + loss = loss_obs + elnet_penalty(beta+t*dbeta, lambda, 0, penalty); ++ls_iter; if (std::isnan(loss)) @@ -381,9 +519,11 @@ void glm_ridge( vec& beta, // output: regression coefficients (contains int if (loss < loss_old - a*t*decrement ) break; } else { - Rcpp::Rcout << "The search direction is not a descent direction "; - Rcpp::Rcout << "(newton decrement = " << decrement << ", should be positive), "; - Rcpp::Rcout << ", this is likely a bug. Please report to the developers." << '\n'; + std::ostringstream msg; + msg << "The search direction is not a descent direction "; + msg << "(newton decrement = " << decrement << ", should be positive), "; + msg << ", this is likely a bug. Please report to the developers."; + msgs.push_back(msg.str()); } } @@ -396,15 +536,15 @@ void glm_ridge( vec& beta, // output: regression coefficients (contains int // Rcpp::Rcout << "tol = " << tol << "\n"; // Rcpp::Rcout << "decrement = " << decrement << "\n"; // Rcpp::Rcout << "\n\n"; - Rcpp::Rcout << "glm_ridge warning: maximum number of line search iterations reached. The optimization can be ill-behaved.\n"; + msgs.push_back("glm_ridge warning: maximum number of line search iterations reached. The optimization can be ill-behaved."); break; } // update the solution beta = beta + t*dbeta; - z = as(obs["z"]); - w = as(obs["wobs"]); - grad_f = as(obs["grad"]); + z = z_ls; + w = w_ls; + grad_f = grad_ls; loss_old = loss; ++qau; } @@ -413,8 +553,10 @@ void glm_ridge( vec& beta, // output: regression coefficients (contains int if (decrement/fabs(loss_initial) > 100*tol) { // warn the user if the max number of iterations is reached and we are relatively far // (two orders of magnitude) from the given convergence threshold - Rcpp::Rcout << "glm_ridge warning: maximum number of quadratic approximation updates reached, within "; - Rcpp::Rcout << decrement << " from optimum (tolerance = " << thresh << ").\n"; + std::ostringstream msg; + msg << "glm_ridge warning: maximum number of quadratic approximation updates reached, within "; + msg << decrement << " from optimum (tolerance = " << thresh << ")."; + msgs.push_back(msg.str()); } } } @@ -445,14 +587,115 @@ List glm_ridge_c( arma::mat x, vec w = w_init; int qau; double loss; - glm_ridge(beta, loss, w, qau, x, pseudo_obs, lambda, intercept, penalty, thresh, qa_updates_max, ls_iter_max, debug); + std::vector msgs; + PseudoObsR po(pseudo_obs); + glm_ridge(beta, loss, w, qau, msgs, x, po, lambda, intercept, penalty, thresh, qa_updates_max, ls_iter_max, debug); if (intercept) - return List::create(vec(beta.tail(D-1)), beta(0), w, loss, qau); + return List::create(vec(beta.tail(D-1)), beta(0), w, loss, qau, wrap(msgs)); // return List::create(vec(beta.tail(D-1)), beta(0), w, qau); else - return List::create(beta, 0.0, w, loss, qau); + return List::create(beta, 0.0, w, loss, qau, wrap(msgs)); // return List::create(beta, 0.0, w, qau); } +/** + * As glm_ridge_c, but with the pseudo-observations computed in C++ for the + * family and link given by their codes (see native_pseudo_codes() in R). + */ +// [[Rcpp::export]] +List glm_ridge_native_c( arma::mat x, + arma::vec y, + arma::vec weights, + arma::vec offset, + int family_code, + int link_code, + double lambda, + bool intercept, + arma::vec penalty, + arma::vec beta_init, + arma::vec w_init, + double thresh, + int qa_updates_max, + int ls_iter_max=100, + bool debug=false) +{ + int D = x.n_cols; + if (intercept) + D++; + + vec beta = beta_init; + vec w = w_init; + int qau; + double loss; + std::vector msgs; + PseudoObsNative po(y, weights, offset, family_code, link_code); + glm_ridge(beta, loss, w, qau, msgs, x, po, lambda, intercept, penalty, thresh, qa_updates_max, ls_iter_max, debug); + + if (intercept) + return List::create(vec(beta.tail(D-1)), beta(0), w, loss, qau, wrap(msgs)); + else + return List::create(beta, 0.0, w, loss, qau, wrap(msgs)); +} + + +/** + * glm_ridge_native_c for several responses (the columns of y) with the same + * model matrix, observation weights, offsets and penalty, as for the draw-wise + * submodels of one projection. Each column is fitted separately, from the same + * initial values. + */ +// [[Rcpp::export]] +List glm_ridge_native_multi_c( arma::mat x, + arma::mat y, + arma::vec weights, + arma::vec offset, + int family_code, + int link_code, + double lambda, + bool intercept, + arma::vec penalty, + arma::vec beta_init, + arma::vec w_init, + double thresh, + int qa_updates_max, + int ls_iter_max=100, + bool debug=false) +{ + int D = x.n_cols; + if (intercept) + D++; + int S = y.n_cols; + int n = y.n_rows; + + mat beta_out(intercept ? D-1 : D, S); + vec beta0_out(S); + mat w_out(n, S); + vec loss_out(S); + IntegerVector qau_out(S); + List msgs_out(S); + for (int s = 0; s < S; ++s) { + vec beta = beta_init; + vec w = w_init; + int qau; + double loss; + std::vector msgs; + PseudoObsNative po(y.col(s), weights, offset, family_code, link_code); + glm_ridge(beta, loss, w, qau, msgs, x, po, lambda, intercept, penalty, thresh, qa_updates_max, ls_iter_max, debug); + if (intercept) { + beta_out.col(s) = beta.tail(D-1); + beta0_out(s) = beta(0); + } else { + beta_out.col(s) = beta; + beta0_out(s) = 0.0; + } + w_out.col(s) = w; + loss_out(s) = loss; + qau_out[s] = qau; + msgs_out[s] = wrap(msgs); + } + return List::create(Named("beta") = beta_out, Named("beta0") = beta0_out, + Named("w") = w_out, Named("loss") = loss_out, + Named("qa_updates") = qau_out, Named("msgs") = msgs_out); +} diff --git a/tests/testthat/test_glm_ridge.R b/tests/testthat/test_glm_ridge.R index c7ccca48c..9e6415031 100644 --- a/tests/testthat/test_glm_ridge.R +++ b/tests/testthat/test_glm_ridge.R @@ -271,6 +271,142 @@ test_that("glm_ridge: poisson, log-link, no intercept, lambda = 0", { # tolerance = tol) # }) +test_that("glm_ridge: warning messages are returned and printed", { + fam <- extend_family(binomial()) + x_w <- x[, 1:2] + y_w <- rbinom(n, 1, fam$linkinv(4 * x_w[, 1] - 3 * x_w[, 2])) + # Too few quadratic approximation updates: + expect_output( + fit_print <- glm_ridge(x_w, y_w, family = fam, qa_updates_max = 2), + "glm_ridge warning: maximum number of quadratic approximation updates" + ) + expect_silent( + fit_quiet <- glm_ridge(x_w, y_w, family = fam, qa_updates_max = 2, + print_msgs = FALSE) + ) + expect_length(fit_quiet$msgs, 1) + expect_match(fit_quiet$msgs, + "^glm_ridge warning: maximum number of quadratic approximation") + expect_identical(fit_print[setdiff(names(fit_print), "msgs")], + fit_quiet[setdiff(names(fit_quiet), "msgs")]) + # Converged fit: + fit_conv <- glm_ridge(x_w, y_w, family = fam, print_msgs = FALSE) + expect_identical(fit_conv$msgs, character()) + # The submodel fitter turns the messages into one R warning: + dat_w <- data.frame(y = y_w, x1 = x_w[, 1], x2 = x_w[, 2]) + expect_warning( + fit_glm_ridge_callback(y ~ x1 + x2, data = dat_w, family = fam, + qa_updates_max = 2), + "maximum number of quadratic approximation updates" + ) +}) + +test_that("glm_ridge: native pseudo-observations equal the R callback", { + with_native <- function(native, expr) { + opts_old <- options(projpred.glm_native = native) + on.exit(options(opts_old)) + expr + } + eta_n <- drop(x[, 1:4] %*% c(0.8, -0.6, 0.4, 0.2)) + fams <- list( + list(binomial("logit"), rbinom(n, weights, plogis(eta_n)) / weights), + list(binomial("probit"), rbinom(n, weights, pnorm(eta_n)) / weights), + list(binomial("cloglog"), + rbinom(n, weights, 1 - exp(-exp(eta_n - 0.5))) / weights), + list(binomial("cauchit"), rbinom(n, weights, pcauchy(eta_n)) / weights), + list(poisson("log"), rpois(n, exp(eta_n / 2))), + list(gaussian("log"), exp(eta_n / 3) + rnorm(n, 0, 0.1)), + list(gaussian("identity"), eta_n + rnorm(n)) + ) + for (fam_y in fams) { + fam <- extend_family(fam_y[[1]]) + y_n <- fam_y[[2]] + expect_false(is.null(native_pseudo_codes(fam))) + for (intcpt in c(TRUE, FALSE)) { + for (x_n in list(x_tr, NULL)) { + if (is.null(x_n) && !intcpt) next + fit_args <- list(x = x_n, y = y_n, family = fam, lambda = 1e-4, + weights = weights, offset = offset, + intercept = intcpt) + fit_native <- with_native(TRUE, do.call(glm_ridge, fit_args)) + fit_r <- with_native(FALSE, do.call(glm_ridge, fit_args)) + info_str <- paste(fam$family, fam$link, intcpt, is.null(x_n)) + expect_equal(fit_native, fit_r, tolerance = 1e-12, info = info_str) + } + } + } + # Families without native pseudo-observations: + expect_null(native_pseudo_codes(extend_family(Gamma()))) + expect_null(native_pseudo_codes(extend_family(poisson("sqrt")))) + fam_mod <- extend_family(binomial()) + fam_mod$linkinv <- function(eta) 1 / (1 + exp(-eta)) + expect_null(native_pseudo_codes(fam_mod)) + expect_null(with_native(FALSE, + native_pseudo_codes(extend_family(binomial())))) +}) + +test_that(paste( + "glm_ridge_native_multi() gives the same submodel fits as the draw-wise", + "fitter" +), { + n_draws <- 4 + eta_d <- x[, 1:3] %*% matrix(rnorm(3 * n_draws, 0, 0.6), 3, n_draws) + fams_d <- list(binomial("logit"), binomial("probit"), poisson("log"), + gaussian("log")) + for (fam_d in fams_d) { + fam <- extend_family(fam_d) + mu_d <- fam$linkinv(eta_d + if (fam$family == "gaussian") 1 else 0) + dat_d <- data.frame(x[, 1:4], xf = factor(rep(c("a", "b"), length.out = n))) + resp_nms <- paste0(".y", seq_len(n_draws)) + dat_d[resp_nms] <- mu_d + fmls <- lapply(resp_nms, function(resp_nm) { + as.formula(paste(resp_nm, "~ X1 + X2 + X3 + X4 + xf")) + }) + shared_x <- shared_model_matrix(fmls, data = dat_d) + for (qa_max in list(NULL, 2)) { + dot_args <- list(data = dat_d, family = fam, weights = weights) + if (!is.null(qa_max)) dot_args$qa_updates_max <- qa_max + out_multi <- glm_ridge_native_multi(fmls, shared_x = shared_x, + dot_args = dot_args) + sdivmin_shared <- glm_ridge_callback_shared(shared_x) + capts_drawwise <- list() + out_drawwise <- lapply(seq_along(fmls), function(s) { + capts_drawwise[[s]] <<- capt_conds( + fit_s <- do.call(sdivmin_shared, c(list(formula = fmls[[s]]), + dot_args)) + ) + fit_s + }) + info_str <- paste(fam$family, fam$link, is.null(qa_max)) + expect_identical(attr(out_multi, "mssgs_warns_capts"), capts_drawwise, + info = info_str) + attr(out_multi, "mssgs_warns_capts") <- NULL + expect_equal(out_multi, out_drawwise, tolerance = 1e-12, info = info_str) + if (!is.null(qa_max) && fam$family != "gaussian") { + expect_true(any(lengths(capts_drawwise) > 0), info = info_str) + } + } + } + # Cases that use the draw-wise fitter: + dot_args <- list(data = dat_d, family = extend_family(Gamma())) + expect_null(glm_ridge_native_multi(fmls, shared_x, dot_args)) + dot_args <- list(data = dat_d, family = extend_family(binomial()), + offset = rep(0, n)) + expect_null(glm_ridge_native_multi(fmls, shared_x, dot_args)) +}) + +test_that("capt_conds() collects messages and warnings", { + expect_silent(capt <- capt_conds({ + message("a message") + warning("line 1\nline 2") + warning("line 1\nline 2") + invisible(1) + })) + expect_identical(capt, c("a message", "line 1", "line 2")) + expect_identical(capt_conds(invisible(1)), character()) + expect_error(capt_conds(stop("an error")), "an error") +}) + if (exists("rng_old")) assign(".Random.seed", rng_old, envir = .GlobalEnv) # Clean up the workspace: rm(list = setdiff(ls(), ls_bu)) From 92f05858019c2d81043113c9d00cc7e1d3f58c93 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Mon, 26 Oct 2026 20:10:24 +0200 Subject: [PATCH 5/5] Add the PR number to the NEWS entry --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 92da0a6d0..5bdbeb947 100644 --- a/NEWS.md +++ b/NEWS.md @@ -9,7 +9,7 @@ If you read this from a place other than