diff --git a/NEWS.md b/NEWS.md index 17f663e2e..4c495e5da 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,6 +7,8 @@ If you read this from a place other than 0 || length(trms$additive_terms) > 0) { + return(NULL) + } + tt <- terms(as.formula(paste("~", paste(terms_all, collapse = " + ")))) + data <- refmodel$fetch_data() + mf <- model.frame(tt, data = data, drop.unused.levels = TRUE) + if (nrow(mf) != nrow(data)) { + return(NULL) + } + fac <- attr(tt, "factors") + is_categ <- vapply(rownames(fac), function(v) { + is.factor(mf[[v]]) || is.character(mf[[v]]) || is.logical(mf[[v]]) + }, logical(1)) + if (any(attr(tt, "order") > 1 & + colSums(fac[is_categ, , drop = FALSE]) > 0)) { + return(NULL) + } + x <- model.matrix(tt, data = mf) + assign_x <- attr(x, "assign") + keep <- assign_x > 0 + x <- x[, keep, drop = FALSE] + assign_x <- assign_x[keep] + cols_by_term <- split(seq_along(assign_x), + norm_label_gauss(attr(tt, "term.labels"))[assign_x]) + wobs <- refmodel$wobs + sum_w <- sum(wobs) + n <- nrow(x) + x_means <- colSums(x * wobs) / sum_w + x_c <- sweep(x, 2, x_means) + g_mat <- crossprod(x_c, x_c * wobs) + sd_x <- sqrt(n / (n - 1) * diag(g_mat) / sum_w) + if (n < 2 || any(!is.finite(sd_x)) || any(sd_x == 0)) { + return(NULL) + } + lambda_x <- regul * sd_x^2 + a_full <- g_mat + diag(a_full) <- diag(a_full) + lambda_x + # For new data, the terms of the model frame keep the data-dependent bases + # of the training data (attribute "predvars", for example of poly()): + return(nlist(tt_pred = attr(mf, "terms"), xlev = .getXlevels(tt, mf), keep, + x_c, x_means, cols_by_term, g_mat, a_full, lambda_x, wobs, + sum_w)) +} + +# The term labels of a term string (see norm_label_gauss()). +labels_gauss <- function(term_string) { + if (identical(term_string, "1")) { + return(character(0)) + } + return(unique(norm_label_gauss(trimws( + strsplit(term_string, "+", fixed = TRUE)[[1]] + )))) +} + +# The performance evaluation of perf_eval() along the predictor ranking for +# Gaussian submodels, from cross-products (see the top of this file): for each +# size in `nterms`, the projected dispersion, the cross-entropy and the +# predictions for the test data, without the draw-wise submodel fits. Returns +# the list `out_by_size` of perf_eval() (elements `sub_summary`, `ce` and +# `dis_sub`), or `NULL` if a condition does not hold. +perf_eval_gauss <- function(search_path, nterms, refmodel, p_ref, + refmodel_fulldata, indices_test, newdata_test, + offset_test, wobs_test, y_test, y_oscale_test, + divmin_args) { + regul <- gauss_crossprod_regul(refmodel, divmin_args) + ranking <- search_path$predictor_ranking + if (is.null(regul) || length(nterms) == 0 || any(nterms > length(ranking)) || + any(offset_test != 0)) { + return(NULL) + } + des <- gauss_crossprod_design(refmodel, c("1", ranking), regul) + if (is.null(des)) { + return(NULL) + } + data_test <- newdata_test %||% refmodel_fulldata$fetch_data(obs = indices_test) + x_test <- tryCatch({ + mf_test <- model.frame(des$tt_pred, data = data_test, xlev = des$xlev, + na.action = na.pass) + model.matrix(des$tt_pred, data = mf_test)[, des$keep, drop = FALSE] + }, error = function(e) NULL) + if (is.null(x_test) || nrow(x_test) != nrow(data_test) || anyNA(x_test)) { + return(NULL) + } + x_test_c <- sweep(x_test, 2, des$x_means) + + # Quantities of the projected reference draws: + wobs <- des$wobs + sum_w <- des$sum_w + mu_ref <- p_ref$mu_offs + ybar <- colSums(mu_ref * wobs) / sum_w + sum_wr2 <- colSums(sweep(mu_ref, 2, ybar)^2 * wobs) + sum_wvar <- colSums(p_ref$var * wobs) + c_mat <- crossprod(des$x_c, mu_ref * wobs) + nprj <- ncol(mu_ref) + + # Columns of the submodel of each size (the submodels are nested): + labels_size <- lapply(seq_along(ranking), function(j) { + unique(unlist(lapply(ranking[seq_len(j)], labels_gauss))) + }) + if (!all(unlist(labels_size) %in% names(des$cols_by_term))) { + return(NULL) + } + + out_size <- function(cols, r_chol, y_mat) { + if (length(cols) == 0) { + quad <- rep(0, nprj) + mu_test <- matrix(ybar, nrow = nrow(x_test), ncol = nprj, byrow = TRUE) + } else { + beta <- backsolve(r_chol, y_mat) + quad <- colSums(y_mat^2) + colSums(des$lambda_x[cols] * beta^2) + mu_test <- sweep(x_test_c[, cols, drop = FALSE] %*% beta, 2, ybar, "+") + } + rownames(mu_test) <- rownames(x_test) + dis <- sqrt((sum_wvar + sum_wr2 - quad) / sum_w) + ce <- weighted.mean((-ybar^2 * sum_w - quad) / sum_w, p_ref$wdraws_prj) + sub_summary <- weighted_summary_means( + y_wobs_test = data.frame(y = y_test, y_oscale = y_oscale_test, + wobs = wobs_test), + data_aux_test = data_test, + family = refmodel$family, + wdraws = p_ref$wdraws_prj, + mu = mu_test, + dis = dis, + cl_ref = p_ref$cl, + wdraws_ref = p_ref$wdraws_orig + ) + return(list(sub_summary = sub_summary, ce = ce, dis_sub = dis)) + } + + # Grow the Cholesky factor along the ranking (block Cholesky): + out_by_size <- vector("list", length(nterms)) + cols <- integer(0) + r_chol <- NULL + y_mat <- NULL + for (size_j in seq(0, max(nterms))) { + if (size_j > 0) { + cols_new <- setdiff( + unlist(des$cols_by_term[labels_size[[size_j]]], use.names = FALSE), + cols + ) + if (length(cols_new) > 0) { + a_nn <- des$a_full[cols_new, cols_new, drop = FALSE] + rhs_n <- c_mat[cols_new, , drop = FALSE] + if (length(cols) > 0) { + r_cn <- forwardsolve(t(r_chol), + des$a_full[cols, cols_new, drop = FALSE]) + a_nn <- a_nn - crossprod(r_cn) + rhs_n <- rhs_n - crossprod(r_cn, y_mat) + } + r_n <- chol(a_nn) + y_n <- forwardsolve(t(r_n), rhs_n) + if (length(cols) == 0) { + r_chol <- r_n + y_mat <- y_n + } else { + r_chol <- rbind( + cbind(r_chol, r_cn), + cbind(matrix(0, length(cols_new), length(cols)), r_n) + ) + y_mat <- rbind(y_mat, y_n) + } + cols <- c(cols, cols_new) + } + } + for (k in which(nterms == size_j)) { + out_by_size[[k]] <- out_size(cols, r_chol, y_mat) + } + } + return(out_by_size) +} diff --git a/R/projfun.R b/R/projfun.R index aa79dbf31..61625eb47 100644 --- a/R/projfun.R +++ b/R/projfun.R @@ -106,7 +106,20 @@ perf_eval <- function(search_path, txt_clust_draws(p_ref[["clust_used"]], p_ref[["nprjdraws"]]), " (`refit_prj = ", refit_prj, "`) ...", verbose = verbose && !return_submodls) - out_by_size <- lapply(nterms, function(size_j) { + out_by_size <- NULL + if (refit_prj && !is.null(reweighting_args) && !return_submodls && + !return_preds && getOption("projpred.search_gauss_crossprod", TRUE)) { + # In the LOO CV with validated search, compute the evaluation of Gaussian + # submodels from cross-products (`NULL` if this case does not apply): + out_by_size <- perf_eval_gauss( + search_path = search_path, nterms = nterms, refmodel = refmodel, + p_ref = p_ref, refmodel_fulldata = refmodel_fulldata, + indices_test = indices_test, newdata_test = newdata_test, + offset_test = offset_test, wobs_test = wobs_test, y_test = y_test, + y_oscale_test = y_oscale_test, divmin_args = list(...) + ) + } + if (is.null(out_by_size)) out_by_size <- lapply(nterms, function(size_j) { # Fetch the init_submodl() output (of class `submodl`) for the submodel at # position `size_j + 1` of the predictor ranking: submodl <- fetch_submodl(size_j, ...) diff --git a/R/projpred-package.R b/R/projpred-package.R index fb1c0508e..a88d9534b 100644 --- a/R/projpred-package.R +++ b/R/projpred-package.R @@ -266,6 +266,12 @@ #' [#207](https://github.com/stan-dev/projpred/issues/207) and GitHub pull #' request [#353](https://github.com/stan-dev/projpred/pull/353)). Currently, #' this is an experimental feature. +#' * `projpred.search_gauss_crossprod`: A single logical value indicating +#' whether to compute the forward search and the evaluation in the LOO CV with +#' validated search for Gaussian submodels with the identity link (including +#' the latent projection) from cross-products of the model matrix (default +#' `TRUE`). With `FALSE`, \pkg{projpred} projects onto each candidate submodel +#' draw by draw. #' #' @name internal-global-options #' @keywords internal diff --git a/R/search.R b/R/search.R index 468898a63..2d9cce0e4 100644 --- a/R/search.R +++ b/R/search.R @@ -1,12 +1,30 @@ +# `return_outdmins = FALSE` tells the search that the caller does not use the +# submodel fits along the ranking (element `outdmins`), so +# search_forward_gauss() does not compute them. The draw-wise search computes +# them anyway. search_forward <- function(p_ref, refmodel, nterms_max, verbose = TRUE, search_terms, est_runtime = TRUE, - search_terms_was_null, ...) { + search_terms_was_null, return_outdmins = TRUE, + ...) { nterms_max_with_icpt <- nterms_max + 1L iq <- ceiling(quantile(seq_len(nterms_max_with_icpt), 1:10 / 10)) if (is.null(search_terms)) { stop("Did not expect `search_terms` to be `NULL`. Please report this.") } + # Gaussian submodels: rank the terms from cross-products (see + # search_forward_gauss()); `NULL` if this case does not apply: + if (getOption("projpred.search_gauss_crossprod", TRUE)) { + out_gauss <- search_forward_gauss( + p_ref = p_ref, refmodel = refmodel, nterms_max = nterms_max, + verbose = verbose, search_terms = search_terms, + return_outdmins = return_outdmins, ... + ) + if (!is.null(out_gauss)) { + return(out_gauss) + } + } + chosen <- character() outdmins <- c() @@ -142,6 +160,193 @@ search_forward <- function(p_ref, refmodel, nterms_max, verbose = TRUE, return(nlist(predictor_ranking = setdiff(chosen, "1"), outdmins)) } +# Forward search for Gaussian submodels with the identity link, from the +# cross-products of the model matrix of all search terms. Let J be the +# model-matrix columns of a candidate. For draw s, the projection solves the +# weighted ridge problem of glm_ridge(): +# (G[J, J] + regul * diag(sd[J]^2)) beta_s = C[J, s]. +# Here, G = Xc' W Xc and C = Xc' W M. Xc is the model matrix, centered at the +# weighted column means. W holds the observation weights, M the fitted values +# of the projected reference draws, and sd the weighted column standard +# deviations from .weighted_sd(). glm_ridge() centers and scales the columns +# in the same way and does not penalize the intercept. The part of ce_gauss() +# that differs between the candidates is +# -beta_s' (G[J, J] + 2 * regul * diag(sd[J]^2)) beta_s / sum(W). +# Thus, the cost per candidate does not depend on the number of observations. +# This function gives only the ranking. The submodel fits along the ranking come +# from proj_to_submodl(), as in search_forward(). The function returns `NULL` +# if one of its conditions does not hold; then search_forward() continues with +# its own search. +search_forward_gauss <- function(p_ref, refmodel, nterms_max, verbose, + search_terms, return_outdmins = TRUE, ...) { + # The arguments for the divergence minimizer, as in proj_to_submodl(): + dot_args <- list(...) + divmin_args <- dot_args$search_control %||% + dot_args[setdiff(names(dot_args), "search_control")] + regul <- gauss_crossprod_regul(refmodel, divmin_args) + if (is.null(regul)) { + return(NULL) + } + # Custom search terms are possible, because the candidates come from + # select_possible_terms_size(), as in search_forward(): + des <- gauss_crossprod_design(refmodel, search_terms, regul) + if (is.null(des)) { + return(NULL) + } + cols_by_term <- des$cols_by_term + c_mat <- crossprod(des$x_c, p_ref$mu * des$wobs) + lambda_x <- des$lambda_x + a_full <- des$a_full + wdraws <- p_ref$wdraws_prj + norm_label <- norm_label_gauss + + # The term labels of a candidate string. The same strings come back in every + # step, so the labels are computed once: + label_cache <- new.env(hash = TRUE, parent = emptyenv()) + labels_of <- function(cand) { + lbls <- label_cache[[cand]] + if (is.null(lbls)) { + lbls <- character(0) + if (!identical(cand, "1")) { + lbls <- unique(norm_label(trimws(strsplit(cand, "+", fixed = TRUE)[[1]]))) + } + assign(cand, lbls, envir = label_cache) + } + return(lbls) + } + + # The criterion is -sum_s wdraws_s * beta_s' (G[J, J] + 2 * Lambda) beta_s + # (smaller is better, as for the cross-entropy). With + # G[J, J] + Lambda = R' R and y_s = R^-T C[J, s], it equals + # -sum_s wdraws_s * (||y_s||^2 + beta_s' Lambda beta_s). In one step, all + # candidates share the columns of the chosen terms, so the Cholesky factor + # `r_chosen` and `y_chosen` of these columns are computed once per step. Each + # candidate extends them by its new columns `cols_new` (block Cholesky). + crit_cand <- function(cols_new, cols_chosen, r_chosen, y_chosen) { + quad <- 0 + r_cn <- NULL + r_n <- NULL + y_n <- NULL + beta_n <- NULL + if (length(cols_new) > 0) { + a_nn <- a_full[cols_new, cols_new, drop = FALSE] + rhs_n <- c_mat[cols_new, , drop = FALSE] + if (length(cols_chosen) > 0) { + r_cn <- forwardsolve(t(r_chosen), + a_full[cols_chosen, cols_new, drop = FALSE]) + a_nn <- a_nn - crossprod(r_cn) + rhs_n <- rhs_n - crossprod(r_cn, y_chosen) + } + r_n <- chol(a_nn) + y_n <- forwardsolve(t(r_n), rhs_n) + beta_n <- backsolve(r_n, y_n) + quad <- colSums(y_n^2) + colSums(lambda_x[cols_new] * beta_n^2) + } + if (length(cols_chosen) > 0) { + y_c <- y_chosen + if (!is.null(beta_n)) { + y_c <- y_c - r_cn %*% beta_n + } + beta_c <- backsolve(r_chosen, y_c) + quad <- quad + colSums(y_chosen^2) + + colSums(lambda_x[cols_chosen] * beta_c^2) + } + return(list(crit = -sum(wdraws * quad), r_cn = r_cn, r_n = r_n, y_n = y_n)) + } + + # crit_cand() for many candidates that each add one column (`cols_single`) + # at once, with matrix operations instead of one call per candidate. For + # candidate j, beta_c = b0 - u_j beta_n[j, ]' with b0 = R^-1 y_chosen and + # u_j = R^-1 r_cj, so the penalty term of the chosen columns is + # sum(lambda b0^2) - 2 beta_n[j, ] p[j, ] + beta_n[j, ]^2 l[j], with + # p = (u * lambda)' b0 and l = colSums(lambda u^2). + crit_single <- function(cols_single, cols_chosen, r_chosen, y_chosen) { + a_jj <- a_full[cbind(cols_single, cols_single)] + rhs <- c_mat[cols_single, , drop = FALSE] + if (length(cols_chosen) > 0) { + r_cj <- forwardsolve(t(r_chosen), + a_full[cols_chosen, cols_single, drop = FALSE]) + a_jj <- a_jj - colSums(r_cj^2) + rhs <- rhs - crossprod(r_cj, y_chosen) + } + y_n <- rhs / sqrt(a_jj) + beta_n <- y_n / sqrt(a_jj) + quad <- y_n^2 + lambda_x[cols_single] * beta_n^2 + if (length(cols_chosen) > 0) { + lambda_c <- lambda_x[cols_chosen] + b0 <- backsolve(r_chosen, y_chosen) + u <- backsolve(r_chosen, r_cj) + q0 <- colSums(y_chosen^2) + colSums(lambda_c * b0^2) + quad <- quad + + matrix(q0, nrow = length(cols_single), ncol = length(q0), byrow = TRUE) - + 2 * beta_n * crossprod(u * lambda_c, b0) + + colSums(lambda_c * u^2) * beta_n^2 + } + return(-drop(quad %*% wdraws)) + } + + chosen <- character() + labels_chosen <- character() + cols_chosen <- integer(0) + r_chosen <- NULL + y_chosen <- NULL + for (size in seq_len(nterms_max + 1L)) { + cands <- select_possible_terms_size(chosen, search_terms, size = size) + if (is.null(cands)) { + next + } + labels_new <- lapply(cands, function(cand) { + setdiff(labels_of(cand), labels_chosen) + }) + if (!all(unlist(labels_new) %in% names(cols_by_term))) { + return(NULL) + } + cols_cands <- lapply(labels_new, function(lbls) { + unlist(cols_by_term[lbls], use.names = FALSE) + }) + crits <- rep(NA_real_, length(cands)) + is_single <- lengths(cols_cands) == 1L + if (any(is_single)) { + crits[is_single] <- crit_single(unlist(cols_cands[is_single]), + cols_chosen, r_chosen, y_chosen) + } + for (k in which(!is_single)) { + crits[k] <- crit_cand(cols_cands[[k]], cols_chosen, r_chosen, + y_chosen)$crit + } + imin <- which.min(crits) + chosen <- c(chosen, cands[imin]) + cols_new <- cols_cands[[imin]] + if (length(cols_new) > 0) { + # The factors of the chosen candidate, for the next step: + best <- crit_cand(cols_new, cols_chosen, r_chosen, y_chosen) + if (length(cols_chosen) == 0) { + r_chosen <- best$r_n + y_chosen <- best$y_n + } else { + r_chosen <- rbind( + cbind(r_chosen, best$r_cn), + cbind(matrix(0, length(cols_new), length(cols_chosen)), best$r_n) + ) + y_chosen <- rbind(y_chosen, best$y_n) + } + cols_chosen <- c(cols_chosen, cols_new) + labels_chosen <- c(labels_chosen, labels_new[[imin]]) + } + } + + # The submodel fits along the ranking (as stored by search_forward()), if the + # caller uses them: + outdmins <- NULL + if (return_outdmins) { + outdmins <- lapply(seq_along(chosen), function(j) { + proj_to_submodl(predictor_terms = chosen[seq_len(j)], p_ref = p_ref, + refmodel = refmodel, verbose = verbose, ...)$outdmin + }) + } + return(nlist(predictor_ranking = setdiff(chosen, "1"), outdmins)) +} + #' Force search terms #' #' A helper function to construct the input for argument `search_terms` of diff --git a/R/varsel.R b/R/varsel.R index 6b71cb7e4..b3a990121 100644 --- a/R/varsel.R +++ b/R/varsel.R @@ -535,7 +535,7 @@ varsel.refmodel <- function( .select <- function(refmodel, ndraws, nclusters, reweighting_args = NULL, method, nterms_max, penalty, verbose, verbose_line_length = 5, verbose_txt_add = "", - search_control, ...) { + search_control, return_outdmins = TRUE, ...) { if (is.null(reweighting_args)) { p_sel <- get_refdist(refmodel, ndraws = ndraws, nclusters = nclusters) } else { @@ -564,7 +564,8 @@ varsel.refmodel <- function( } else if (method == "forward") { search_path <- search_forward( p_ref = p_sel, refmodel = refmodel, nterms_max = nterms_max, - verbose = verbose, search_control = search_control, ... + verbose = verbose, search_control = search_control, + return_outdmins = return_outdmins, ... ) } verb_out(rep("-", verbose_line_length), verbose = verbose) diff --git a/man/internal-global-options.Rd b/man/internal-global-options.Rd index e3b9d2560..8d74685eb 100644 --- a/man/internal-global-options.Rd +++ b/man/internal-global-options.Rd @@ -40,6 +40,12 @@ projections (see GitHub issue \href{https://github.com/stan-dev/projpred/issues/207}{#207} and GitHub pull request \href{https://github.com/stan-dev/projpred/pull/353}{#353}). Currently, this is an experimental feature. +\item \code{projpred.search_gauss_crossprod}: A single logical value indicating +whether to compute the forward search and the evaluation in the LOO CV with +validated search for Gaussian submodels with the identity link (including +the latent projection) from cross-products of the model matrix (default +\code{TRUE}). With \code{FALSE}, \pkg{projpred} projects onto each candidate submodel +draw by draw. } } \keyword{internal} diff --git a/tests/testthat/test_search_gauss.R b/tests/testthat/test_search_gauss.R new file mode 100644 index 000000000..98aad939c --- /dev/null +++ b/tests/testthat/test_search_gauss.R @@ -0,0 +1,213 @@ +context("forward search from cross-products for Gaussian submodels") + +# A Gaussian reference model with simulated draws of the coefficients (no model +# fit is necessary). `fml` gives the predictor terms. +refm_search_gauss <- function(fml, n = 80, ndraws = 200, seed = 6283) { + set.seed(seed) + dat <- data.frame(x1 = rnorm(n), x2 = rnorm(n), x3 = rnorm(n), + x4 = rnorm(n), + f = factor(sample(c("a", "b", "c"), n, replace = TRUE))) + dat$y <- 0.8 * dat$x1 - 0.4 * dat$x2 + 0.3 * dat$x1 * dat$x2 + + c(a = 0, b = 0.5, c = -0.5)[as.character(dat$f)] + rnorm(n) + mm <- model.matrix(fml, dat) + beta_hat <- drop(solve(crossprod(mm) + diag(1e-6, ncol(mm)), + crossprod(mm, dat$y))) + beta <- beta_hat + matrix(rnorm(ncol(mm) * ndraws, sd = 0.05), ncol(mm)) + fit <- list(beta = beta, data = dat, fml = fml) + ref_predfun <- function(fit, newdata = NULL) { + nd <- if (is.null(newdata)) fit$data else newdata + model.matrix(delete.response(terms(fit$fml)), nd) %*% fit$beta + } + suppressWarnings(init_refmodel(fit, data = dat, formula = fml, + family = gaussian(), ref_predfun = ref_predfun, + dis = rep(1, ndraws))) +} + +fml_search_gauss <- y ~ x1 + x2 + x3 + x4 + f + x1:x2 + poly(x3, 2) + +test_that(paste( + "the cross-product search gives the same results as the search with", + "draw-wise projections" +), { + refm <- refm_search_gauss(fml_search_gauss) + for (nclusters_crr in c(1, 5)) { + vs_crossprod <- suppressWarnings(varsel( + refm, nclusters = nclusters_crr, nclusters_pred = 7, nterms_max = 6, + seed = 3317, verbose = FALSE + )) + opts_old <- options(projpred.search_gauss_crossprod = FALSE) + on.exit(options(opts_old), add = TRUE) + vs_drawwise <- suppressWarnings(varsel( + refm, nclusters = nclusters_crr, nclusters_pred = 7, nterms_max = 6, + seed = 3317, verbose = FALSE + )) + options(opts_old) + expect_identical(ranking(vs_crossprod)$fulldata, + ranking(vs_drawwise)$fulldata, info = nclusters_crr) + expect_equal(performances(vs_crossprod, deltas = FALSE), + performances(vs_drawwise, deltas = FALSE), + info = nclusters_crr) + expect_equal(vs_crossprod$search_path$outdmins, + vs_drawwise$search_path$outdmins, info = nclusters_crr) + } +}) + +test_that(paste( + "the LOO CV with validated search gives the same results with and without", + "the cross-product computations" +), { + refm <- refm_search_gauss(fml_search_gauss, n = 40) + args_cv <- list(refm, cv_method = "LOO", validate_search = TRUE, nloo = 10, + nclusters = 3, ndraws_pred = 20, nterms_max = 4, seed = 5821, + verbose = FALSE) + cv_crossprod <- suppressWarnings(do.call(cv_varsel, args_cv)) + opts_old <- options(projpred.search_gauss_crossprod = FALSE) + on.exit(options(opts_old), add = TRUE) + cv_drawwise <- suppressWarnings(do.call(cv_varsel, args_cv)) + options(opts_old) + expect_identical(ranking(cv_crossprod), ranking(cv_drawwise)) + expect_equal(cv_crossprod$summaries$sub, cv_drawwise$summaries$sub, + tolerance = 1e-10) + expect_equal(performances(cv_crossprod, deltas = FALSE), + performances(cv_drawwise, deltas = FALSE), tolerance = 1e-10) +}) + +test_that("perf_eval_gauss() applies in the LOO CV with validated search", { + refm <- refm_search_gauss(fml_search_gauss, n = 40) + n_used <- 0L + trace("perf_eval_gauss", exit = function() { + if (!is.null(returnValue())) n_used <<- n_used + 1L + }, print = FALSE, where = asNamespace("projpred")) + on.exit(untrace("perf_eval_gauss", where = asNamespace("projpred"))) + cv <- suppressWarnings(cv_varsel( + refm, cv_method = "LOO", validate_search = TRUE, nloo = 5, nclusters = 3, + ndraws_pred = 20, nterms_max = 3, seed = 5821, verbose = FALSE + )) + expect_identical(n_used, 5L) +}) + +test_that("search_forward_gauss() applies to Gaussian GLM-type submodels", { + refm <- refm_search_gauss(fml_search_gauss) + p_ref <- get_refdist(refm, nclusters = 3) + search_terms <- union("1", split_formula(refm$formula, + data = refm$fetch_data())) + out <- search_forward_gauss(p_ref = p_ref, refmodel = refm, nterms_max = 4, + verbose = FALSE, search_terms = search_terms) + expect_false(is.null(out)) + expect_length(out$predictor_ranking, 4) + expect_length(out$outdmins, 5) +}) + +test_that(paste( + "search_forward_gauss() with `return_outdmins = FALSE` gives the same", + "ranking and no submodel fits" +), { + refm <- refm_search_gauss(fml_search_gauss) + p_ref <- get_refdist(refm, nclusters = 3) + search_terms <- union("1", split_formula(refm$formula, + data = refm$fetch_data())) + out <- search_forward_gauss(p_ref = p_ref, refmodel = refm, nterms_max = 4, + verbose = FALSE, search_terms = search_terms) + out_nofits <- search_forward_gauss(p_ref = p_ref, refmodel = refm, + nterms_max = 4, verbose = FALSE, + search_terms = search_terms, + return_outdmins = FALSE) + expect_identical(out_nofits$predictor_ranking, out$predictor_ranking) + expect_null(out_nofits$outdmins) +}) + +test_that(paste( + "the search of each left-out observation in the LOO CV with validated", + "search does not compute the submodel fits along the ranking" +), { + refm <- refm_search_gauss(y ~ x1 + x2 + x3 + x4 + f + x1:x2, n = 40) + calls <- new.env() + calls$n_nofits <- 0L + ns <- asNamespace("projpred") + suppressMessages(trace( + "search_forward_gauss", where = ns, print = FALSE, + tracer = bquote(if (!return_outdmins) { + assign("n_nofits", get("n_nofits", envir = .(calls)) + 1L, + envir = .(calls)) + }) + )) + on.exit(suppressMessages(untrace("search_forward_gauss", where = ns)), + add = TRUE) + cv <- suppressWarnings(cv_varsel( + refm, cv_method = "LOO", validate_search = TRUE, nloo = 5, nclusters = 3, + ndraws_pred = 20, nterms_max = 3, seed = 5821, verbose = FALSE + )) + # One search per left-out observation; the search on the full data keeps + # its submodel fits: + expect_identical(calls$n_nofits, 5L) +}) + +test_that(paste( + "the cross-product search and evaluation apply to the latent projection", + "and give the same results as the draw-wise projections" +), { + set.seed(7391) + n <- 50 + dat <- data.frame(x1 = rnorm(n), x2 = rnorm(n), x3 = rnorm(n), + x4 = rnorm(n)) + eta <- 0.9 * dat$x1 - 0.6 * dat$x2 + dat$y <- rbinom(n, 1, plogis(eta)) + ndraws <- 100 + mm <- model.matrix(~ x1 + x2 + x3 + x4, dat) + beta <- c(0, 0.9, -0.6, 0, 0) + + matrix(rnorm(ncol(mm) * ndraws, sd = 0.1), 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 + x3 + x4, nd) %*% fit$beta + } + refm <- suppressWarnings(init_refmodel( + fit, data = dat, formula = y ~ x1 + x2 + x3 + x4, family = binomial(), + ref_predfun = ref_predfun, dis = rep(1, ndraws), latent = TRUE + )) + expect_true(refm$family$for_latent) + calls <- new.env() + calls$n <- 0L + suppressMessages(trace( + "search_forward_gauss", where = asNamespace("projpred"), print = FALSE, + exit = bquote(if (!is.null(returnValue(NULL))) { + assign("n", get("n", envir = .(calls)) + 1L, envir = .(calls)) + }) + )) + on.exit(suppressMessages(untrace("search_forward_gauss", + where = asNamespace("projpred"))), + add = TRUE) + args_vs <- list(refm, nclusters = 3, nclusters_pred = 5, nterms_max = 3, + seed = 4410, verbose = FALSE) + args_cv <- c(args_vs, list(cv_method = "LOO", validate_search = TRUE, + nloo = 5)) + vs_crossprod <- suppressWarnings(do.call(varsel, args_vs)) + cv_crossprod <- suppressWarnings(do.call(cv_varsel, args_cv)) + opts_old <- options(projpred.search_gauss_crossprod = FALSE) + on.exit(options(opts_old), add = TRUE) + vs_drawwise <- suppressWarnings(do.call(varsel, args_vs)) + cv_drawwise <- suppressWarnings(do.call(cv_varsel, args_cv)) + options(opts_old) + expect_identical(ranking(vs_crossprod)$fulldata, + ranking(vs_drawwise)$fulldata) + expect_equal(performances(vs_crossprod, deltas = FALSE), + performances(vs_drawwise, deltas = FALSE)) + expect_identical(ranking(cv_crossprod), ranking(cv_drawwise)) + expect_equal(performances(cv_crossprod, deltas = FALSE), + performances(cv_drawwise, deltas = FALSE), tolerance = 1e-10) + # The cross-product search ran: once in varsel(), once for the full data in + # cv_varsel() and once per left-out observation: + expect_identical(calls$n, 7L) +}) + +test_that(paste( + "search_forward_gauss() returns `NULL` for a factor in an interaction term" +), { + refm <- refm_search_gauss(y ~ x1 + x2 + f + x1:f) + p_ref <- get_refdist(refm, nclusters = 3) + search_terms <- union("1", split_formula(refm$formula, + data = refm$fetch_data())) + out <- search_forward_gauss(p_ref = p_ref, refmodel = refm, nterms_max = 3, + verbose = FALSE, search_terms = search_terms) + expect_null(out) +})