From 95e282523e3a8f6c1c5bc6d3fe3050327a15da6c Mon Sep 17 00:00:00 2001 From: fweber144 Date: Thu, 26 Sep 2024 23:16:25 +0200 Subject: [PATCH 01/40] Adapt the existing tests to work with the new implementation of subsampled LOO-CV. --- tests/testthat/test_varsel.R | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test_varsel.R b/tests/testthat/test_varsel.R index a4c580683..6db89b676 100644 --- a/tests/testthat/test_varsel.R +++ b/tests/testthat/test_varsel.R @@ -1384,6 +1384,10 @@ test_that("invalid `nloo` fails", { skip_if_not(run_cvvs) tstsetups_nonkfold <- grep("\\.kfold", names(cvvss), value = TRUE, invert = TRUE) + valsearch_arg <- lapply(args_cvvs[tstsetups_nonkfold], "[[", + "validate_search") + tstsetups_nonkfold <- tstsetups_nonkfold[!sapply(valsearch_arg, isFALSE)] + skip_if(length(tstsetups_nonkfold) == 0) for (tstsetup in head(tstsetups_nonkfold, 1)) { args_cvvs_i <- args_cvvs[[tstsetup]] # Use suppressWarnings() because test_that() somehow redirects stderr() and @@ -1410,6 +1414,9 @@ test_that(paste( "\\.glm\\.gauss\\..*\\.default_cvmeth\\.default_search_trms", names(cvvss), value = TRUE ) + valsearch_arg <- lapply(args_cvvs[tstsetups], "[[", "validate_search") + tstsetups <- tstsetups[!sapply(valsearch_arg, isFALSE)] + skip_if(length(tstsetups) == 0) for (tstsetup in tstsetups) { args_cvvs_i <- args_cvvs[[tstsetup]] # Use suppressWarnings() because test_that() somehow redirects stderr() and @@ -1428,7 +1435,8 @@ test_that("setting `nloo` smaller than the number of observations works", { skip_if_not(run_cvvs) nloo_tst <- nobsv %/% 5L # Output elements of `vsel` objects that may be influenced by `nloo`: - vsel_nms_nloo <- c("summaries", "summaries_fast","predictor_ranking_cv", "nloo", "loo_inds", "ce") + vsel_nms_nloo <- c("summaries", "summaries_fast","predictor_ranking_cv", + "nloo", "loo_inds", "ce") # In general, element `ce` is affected as well (because the PRNG state when # doing the clustering for the performance evaluation is different when `nloo` # is smaller than the number of observations compared to when `nloo` is equal @@ -1440,6 +1448,9 @@ test_that("setting `nloo` smaller than the number of observations works", { "\\.glm\\.gauss\\..*\\.default_cvmeth\\.default_search_trms", names(cvvss), value = TRUE ) + valsearch_arg <- lapply(args_cvvs[tstsetups], "[[", "validate_search") + tstsetups <- tstsetups[!sapply(valsearch_arg, isFALSE)] + skip_if(length(tstsetups) == 0) for (tstsetup in tstsetups) { args_cvvs_i <- args_cvvs[[tstsetup]] tstsetup_ref <- args_cvvs_i$tstsetup_ref @@ -1461,7 +1472,7 @@ test_that("setting `nloo` smaller than the number of observations works", { prd_trms_len_expected = args_cvvs_i$nterms_max, method_expected = meth_exp_crr, cv_method_expected = "LOO", - valsearch_expected = args_cvvs_i$validate_search, + valsearch_expected = TRUE, nloo_expected = nloo_tst, search_terms_expected = args_cvvs_i$search_terms, search_trms_empty_size = From c101badab002f68b5eededdfec4f91d4b21ff70a Mon Sep 17 00:00:00 2001 From: fweber144 Date: Mon, 4 Nov 2024 22:16:54 +0100 Subject: [PATCH 02/40] fix the early check for all-`NA`s in `get_stat()` --- R/summary_funs.R | 23 ++++++++++++++++++++--- 1 file changed, 20 insertions(+), 3 deletions(-) diff --git a/R/summary_funs.R b/R/summary_funs.R index de02ccf85..42c09dd7c 100644 --- a/R/summary_funs.R +++ b/R/summary_funs.R @@ -278,11 +278,28 @@ get_stat <- function(summaries, summaries_baseline = NULL, y_wobs_test, stat, alpha = 0.1, ...) { mu <- summaries$mu lppd <- summaries$lppd - n_full <- length(lppd) - n_loo <- if (is.null(loo_inds)) n_full else length(loo_inds) - if (all(is.na(lppd)) || all(is.na(y_wobs_test$y_prop %||% y_wobs_test$y))) { + all_na_baseline <- NULL + if (stat %in% c("elpd", "mlpd", "gmpd")) { + if (!is.null(summaries_baseline)) { + all_na_baseline <- all(is.na(summaries_baseline$lppd)) + } + all_na <- all(is.na(lppd)) + } else { + hasNA_y <- is.na(y_wobs_test$y_prop %||% y_wobs_test$y) + if (!is.null(summaries_baseline)) { + all_na_baseline <- all(is.na(summaries_baseline$mu) | hasNA_y) + } + all_na <- all(is.na(mu) | hasNA_y) + } + if (!is.null(all_na_baseline) && + getOption("projpred.additional_checks", FALSE)) { + stopifnot(all_na == all_na_baseline) + } + if (all_na) { return(list(value = NA, se = NA, lq = NA, uq = NA)) } + n_full <- length(lppd) + n_loo <- if (is.null(loo_inds)) n_full else length(loo_inds) alpha_half <- alpha / 2 one_minus_alpha_half <- 1 - alpha_half From 06cae7d601f25372322a20e5e031c0d3bf6c55c7 Mon Sep 17 00:00:00 2001 From: fweber144 Date: Mon, 4 Nov 2024 22:31:22 +0100 Subject: [PATCH 03/40] Revert "fix the early check for all-`NA`s in `get_stat()`" This reverts commit 0bd3830f4c87548335ea33819e52193c301d7940. Reason for the revert: The "fix" from that commit would have to take `summaries_fast` into account as well. For simplicity, we will avoid the early check for all-`NA`s in `get_stat()` in its entirety. --- R/summary_funs.R | 23 +++-------------------- 1 file changed, 3 insertions(+), 20 deletions(-) diff --git a/R/summary_funs.R b/R/summary_funs.R index 42c09dd7c..de02ccf85 100644 --- a/R/summary_funs.R +++ b/R/summary_funs.R @@ -278,28 +278,11 @@ get_stat <- function(summaries, summaries_baseline = NULL, y_wobs_test, stat, alpha = 0.1, ...) { mu <- summaries$mu lppd <- summaries$lppd - all_na_baseline <- NULL - if (stat %in% c("elpd", "mlpd", "gmpd")) { - if (!is.null(summaries_baseline)) { - all_na_baseline <- all(is.na(summaries_baseline$lppd)) - } - all_na <- all(is.na(lppd)) - } else { - hasNA_y <- is.na(y_wobs_test$y_prop %||% y_wobs_test$y) - if (!is.null(summaries_baseline)) { - all_na_baseline <- all(is.na(summaries_baseline$mu) | hasNA_y) - } - all_na <- all(is.na(mu) | hasNA_y) - } - if (!is.null(all_na_baseline) && - getOption("projpred.additional_checks", FALSE)) { - stopifnot(all_na == all_na_baseline) - } - if (all_na) { - return(list(value = NA, se = NA, lq = NA, uq = NA)) - } n_full <- length(lppd) n_loo <- if (is.null(loo_inds)) n_full else length(loo_inds) + if (all(is.na(lppd)) || all(is.na(y_wobs_test$y_prop %||% y_wobs_test$y))) { + return(list(value = NA, se = NA, lq = NA, uq = NA)) + } alpha_half <- alpha / 2 one_minus_alpha_half <- 1 - alpha_half From b65d3417a46ac1ccc02227e19caed09949641bc5 Mon Sep 17 00:00:00 2001 From: fweber144 Date: Mon, 4 Nov 2024 22:40:45 +0100 Subject: [PATCH 04/40] avoid the early check for all-`NA`s in `get_stat()` --- R/misc.R | 10 +- R/summary_funs.R | 33 +++-- tests/testthat/helpers/creators.R | 7 ++ tests/testthat/setup.R | 1 + tests/testthat/test_misc.R | 192 +++++++++++++++++++++++++++++- 5 files changed, 224 insertions(+), 19 deletions(-) create mode 100644 tests/testthat/helpers/creators.R diff --git a/R/misc.R b/R/misc.R index 2b31d5d0c..e09f95b91 100644 --- a/R/misc.R +++ b/R/misc.R @@ -68,11 +68,10 @@ ilinkfun_raw <- function(x, link_nm) { pred <- x[, 2] wobs <- x[, 3] - # Make it explicit that `x` should not be used anymore (due to the possibility - # of `NA`s, but also due to the re-ordering): + # Make it explicit that `x` should not be used anymore: rm(x) - ord <- order(pred, decreasing = TRUE, na.last = NA) + ord <- order(pred, decreasing = TRUE, na.last = FALSE) n <- length(ord) resp <- resp[ord] @@ -80,8 +79,9 @@ ilinkfun_raw <- function(x, link_nm) { wobs <- wobs[ord] w0 <- w1 <- wobs - # CAUTION: The following check also ensures that `resp` does not have `NA`s: - stopifnot(all(resp %in% c(0, 1))) + stopifnot(all(na.omit(resp) %in% c(0, 1))) + w0[is.na(resp)] <- NA # ensure that `NA`s in `resp` propagate to the output + w1[is.na(resp)] <- NA # ensure that `NA`s in `resp` propagate to the output w0[resp == 1] <- 0 # for calculating the false positive rate (fpr) w1[resp == 0] <- 0 # for calculating the true positive rate (tpr) cum_w0 <- cumsum(w0) diff --git a/R/summary_funs.R b/R/summary_funs.R index de02ccf85..1faa42450 100644 --- a/R/summary_funs.R +++ b/R/summary_funs.R @@ -280,9 +280,6 @@ get_stat <- function(summaries, summaries_baseline = NULL, lppd <- summaries$lppd n_full <- length(lppd) n_loo <- if (is.null(loo_inds)) n_full else length(loo_inds) - if (all(is.na(lppd)) || all(is.na(y_wobs_test$y_prop %||% y_wobs_test$y))) { - return(list(value = NA, se = NA, lq = NA, uq = NA)) - } alpha_half <- alpha / 2 one_minus_alpha_half <- 1 - alpha_half @@ -302,8 +299,8 @@ get_stat <- function(summaries, summaries_baseline = NULL, value_se <- sqrt(srs_diffe$v_y_hat + srs_diffe$hat_v_y) } else { # full LOO estimator - value <- sum(lppd - lppd_baseline, na.rm = TRUE) - value_se <-sd(lppd - lppd_baseline, na.rm = TRUE) * sqrt(n_full) + value <- sum(lppd - lppd_baseline) + value_se <- sd(lppd - lppd_baseline) * sqrt(n_full) } if (stat %in% c("mlpd", "gmpd")) { value <- value / n_full @@ -517,18 +514,28 @@ get_stat <- function(summaries, summaries_baseline = NULL, }, ... ) - value_se <- sd(diffvalue.bootstrap, na.rm = TRUE) - lq_uq <- quantile(diffvalue.bootstrap, - probs = c(alpha_half, one_minus_alpha_half), - names = FALSE, na.rm = TRUE) + value_se <- sd(diffvalue.bootstrap) + if (any(is.na(diffvalue.bootstrap))) { + # quantile() is not able to deal with `NA`s + lq_uq <- rep(NA_real_, 2) + } else { + lq_uq <- quantile(diffvalue.bootstrap, + probs = c(alpha_half, one_minus_alpha_half), + names = FALSE) + } } else { auc_data <- cbind(y, mu, wobs) value <- .auc(auc_data) value.bootstrap <- bootstrap(auc_data, .auc, ...) - value_se <- sd(value.bootstrap, na.rm = TRUE) - lq_uq <- quantile(value.bootstrap, - probs = c(alpha_half, one_minus_alpha_half), - names = FALSE, na.rm = TRUE) + value_se <- sd(value.bootstrap) + if (any(is.na(value.bootstrap))) { + # quantile() is not able to deal with `NA`s + lq_uq <- rep(NA_real_, 2) + } else { + lq_uq <- quantile(value.bootstrap, + probs = c(alpha_half, one_minus_alpha_half), + names = FALSE) + } } } } diff --git a/tests/testthat/helpers/creators.R b/tests/testthat/helpers/creators.R new file mode 100644 index 000000000..59f1bae3f --- /dev/null +++ b/tests/testthat/helpers/creators.R @@ -0,0 +1,7 @@ +imperfect_alternation <- function(unq_vals, n_tail = 6L, ...) { + x <- rep_len(unq_vals, ...) + tail_idxs_x <- tail(seq_along(x), n = n_tail) + x[tail_idxs_x] <- x[c(tail_idxs_x[-1], + tail_idxs_x[1])] + return(x) +} diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R index c0181c176..733e6618e 100755 --- a/tests/testthat/setup.R +++ b/tests/testthat/setup.R @@ -163,6 +163,7 @@ source(testthat::test_path("helpers", "getters.R"), local = TRUE) source(testthat::test_path("helpers", "formul_handlers.R"), local = TRUE) source(testthat::test_path("helpers", "predictor_handlers.R"), local = TRUE) source(testthat::test_path("helpers", "dummies.R"), local = TRUE) +source(testthat::test_path("helpers", "creators.R"), local = TRUE) # Note: The following `mod_nms` refer to *generalized* (linear/additive, # multilevel) models. This is due to history (when these tests were written, diff --git a/tests/testthat/test_misc.R b/tests/testthat/test_misc.R index d9f0f7768..982ac5e42 100644 --- a/tests/testthat/test_misc.R +++ b/tests/testthat/test_misc.R @@ -168,7 +168,7 @@ test_that("rstanarm: special formulas work", { } }) -# pseudo_data() ----------------------------------------------------------- +# Other internal functions ------------------------------------------------ test_that(paste( "`pseudo_data(f = 0, [...], family = extend_family(gaussian()), [...])` is", @@ -188,3 +188,193 @@ test_that(paste( expect_false(isTRUE(all.equal(psdat$z, mu_crr))) expect_equal(psdat$wobs, wobs_crr) }) + +test_that(".srs_diff_est_w() propagates input `NA`s to its output", { + nloo_tst <- nobsv %/% 5L + loo_inds_tst <- ceiling(seq(1L, nobsv, length.out = nloo_tst)) + skip_if_not(length(unique(loo_inds_tst)) == nloo_tst) + srs_diff_res_NAs <- list( + m = nloo_tst, + N = nobsv, + y_hat = NA_real_, + v_y_hat = NA_real_, + hat_v_y = NA_real_ + ) + expect_identical( + .srs_diff_est_w(y_approx = rep(NA, nobsv), + y = rep(NA, nloo_tst), + y_idx = loo_inds_tst, + wobs = rep(NA, nobsv)), + srs_diff_res_NAs, + info = "all inputs (except `y_idx`) are all-`NA`s" + ) + expect_identical( + .srs_diff_est_w(y_approx = rep(-0.8, nobsv), + y = rep(NA, nloo_tst), + y_idx = loo_inds_tst, + wobs = rep_len(c(2, 3), length.out = nobsv)), + srs_diff_res_NAs, + info = "`y_approx` without `NA`s, `y` only `NA`s, `wobs` without `NA`s" + ) + expect_identical( + .srs_diff_est_w(y_approx = rep(-0.8, nobsv), + y = c(rep(-0.7, nloo_tst - 1L), NA), + y_idx = loo_inds_tst, + wobs = rep_len(c(2, 3), length.out = nobsv)), + srs_diff_res_NAs, + info = paste0("`y_approx` without `NA`s, `y` with a single `NA`, ", + "`wobs` without `NA`s") + ) + expect_identical( + .srs_diff_est_w(y_approx = c(NA, rep(-0.8, nobsv - 1L)), + y = c(rep(NA, nloo_tst - 1L), -0.7), + y_idx = loo_inds_tst, + wobs = rep_len(c(2, 3), length.out = nobsv)), + srs_diff_res_NAs, + info = paste0("`y_approx` with a single `NA`, `y` with a single `NA`, ", + "`wobs` without `NA`s") + ) +}) + +test_that("weighted.sd() with `na.rm = FALSE` propagates input `NA`s to its output", { + expect_identical( + weighted.sd(x = rep(NA, nobsv), + w = rep(NA, nobsv)), + NA_real_, + info = "all inputs are all-`NA`s" + ) + expect_identical( + weighted.sd(x = rep(-0.8, nobsv), + w = c(rep_len(c(2, 3), length.out = nobsv - 1L), NA)), + NA_real_, + info = "`x` without `NA`s, `w` with a single `NA`" + ) + expect_identical( + weighted.sd(x = c(rep(-0.8, nobsv - 1L), NA), + w = rep_len(c(2, 3), length.out = nobsv)), + NA_real_, + info = "`x` with a single `NA`, `w` without `NA`s" + ) + expect_identical( + weighted.sd(x = rep(-0.8, nobsv), + w = rep(NA, nobsv)), + NA_real_, + info = "`x` without `NA`s, `w` all-`NA`s" + ) + expect_identical( + weighted.sd(x = rep(NA, nobsv), + w = rep_len(c(2, 3), length.out = nobsv)), + NA_real_, + info = "`x` all-`NA`s, `w` without `NA`s" + ) +}) + +test_that("auc() works", { + nobsv_auc <- 19L + expect_equal( + auc(cbind(rep_len(c(0, 1), length.out = nobsv_auc), + imperfect_alternation(c(0.3, 0.8), length.out = nobsv_auc), + rep(1, nobsv_auc))), + 0.6833333333333333333333, + info = "`wobs` column only `1`s" + ) + expect_equal( + auc(cbind(rep_len(c(0, 1), length.out = nobsv_auc), + imperfect_alternation(c(0.3, 0.8), length.out = nobsv_auc), + imperfect_alternation(c(1, 2), n_tail = 11L, length.out = nobsv_auc))), + 0.717948717948718062587, + info = "`wobs` column with imperfect alternation of `1`s and `2`s" + ) +}) + +test_that("auc() propagates input `NA`s to its output", { + nobsv_auc <- 19L + expect_identical( + auc(cbind(rep(NA, nobsv_auc), + rep(NA, nobsv_auc), + rep(NA, nobsv_auc))), + NA_real_, + info = "all inputs are all-`NA`s" + ) + expect_identical( + auc(cbind(rep(1, nobsv_auc), + rep(NA, nobsv_auc), + rep(1, nobsv_auc))), + NA_real_, + info = paste0("`resp` column without `NA`s, ", + "`pred` column only `NA`s, ", + "`wobs` column without `NA`s") + ) + expect_identical( + auc(cbind(rep(1, nobsv_auc), + rep(NA, nobsv_auc), + c(rep(1, nobsv_auc - 1L), NA))), + NA_real_, + info = paste0("`resp` column without `NA`s, ", + "`pred` column only `NA`s, ", + "`wobs` column with a single `NA`") + ) + expect_identical( + auc(cbind(c(rep_len(c(0, 1), length.out = nobsv_auc - 1L), NA), + rep(0.7, nobsv_auc), + rep(1, nobsv_auc))), + NA_real_, + info = paste0("`resp` column with a single `NA`, ", + "`pred` column without `NA`s, ", + "`wobs` column without `NA`s") + ) + expect_identical( + auc(cbind(rep(NA, nobsv_auc), + rep(0.7, nobsv_auc), + rep(1, nobsv_auc))), + NA_real_, + info = paste0("`resp` column only `NA`s, ", + "`pred` column without `NA`s, ", + "`wobs` column without `NA`s") + ) + expect_identical( + auc(cbind(c(rep_len(c(0, 1), length.out = nobsv_auc - 1L), NA), + c(rep(0.7, nobsv_auc - 1L), NA), + rep(1, nobsv_auc))), + NA_real_, + info = paste0("`resp` column with a single `NA`, ", + "`pred` column with a single `NA`, ", + "`wobs` column without `NA`s") + ) + expect_identical( + auc(cbind(c(rep_len(c(0, 1), length.out = nobsv_auc - 1L), NA), + c(NA, rep(0.7, nobsv_auc - 1L)), + rep(1, nobsv_auc))), + NA_real_, + info = paste0("`resp` column with a single `NA`, ", + "`pred` column with a single `NA` (at different position), ", + "`wobs` column without `NA`s") + ) + expect_identical( + auc(cbind(c(rep_len(c(0, 1), length.out = nobsv_auc - 1L), NA), + rep(NA, nobsv_auc), + rep(1, nobsv_auc))), + NA_real_, + info = paste0("`resp` column with a single `NA`, ", + "`pred` column only `NA`s, ", + "`wobs` column without `NA`s") + ) + expect_identical( + auc(cbind(rep(NA, nobsv_auc), + rep(0.7, nobsv_auc), + c(rep(1, nobsv_auc - 1L), NA))), + NA_real_, + info = paste0("`resp` column only `NA`s, ", + "`pred` column without `NA`s, ", + "`wobs` column with a single `NA`") + ) + expect_identical( + auc(cbind(rep_len(c(0, 1), length.out = nobsv_auc), + imperfect_alternation(c(0.3, 0.8), length.out = nobsv_auc), + rep(NA, nobsv_auc))), + NA_real_, + info = paste0("`resp` column without `NA`s, ", + "`pred` column without `NA`s, ", + "`wobs` column only `NA`s") + ) +}) From bd6673419c57cfe6a20de7f08ee95d7fb01f9fc6 Mon Sep 17 00:00:00 2001 From: fweber144 Date: Sun, 17 Nov 2024 14:56:50 +0100 Subject: [PATCH 05/40] replace `n_full <- sum(!is.na(mu))` with `n_full <- length(mu)` because this ensures `all(wobs == 1)` after lines ``` wobs <- rep(wobs, y_wobs_test$wobs) wobs <- n_full * wobs / sum(wobs) ``` (but note that `NA`s in `mu` should lead to `NA` output anyway). --- R/summary_funs.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/summary_funs.R b/R/summary_funs.R index 1faa42450..234b83ce6 100644 --- a/R/summary_funs.R +++ b/R/summary_funs.R @@ -438,7 +438,7 @@ get_stat <- function(summaries, summaries_baseline = NULL, mu_fast <- rep(summaries_fast$mu, y_wobs_test$wobs) # CAUTION: If `y` is allowed to have `NA`s here, then the following # definition of `n_full` needs to be adapted: - n_full <- sum(!is.na(mu)) + n_full <- length(mu) if (!is.null(loo_inds)) { stopifnot(all(y_wobs_test$wobs > 0)) loo_inds <- unlist(lapply(loo_inds, function(loo_idx) { From af7bd97c064642714e66ee0c0e67d21accdab43b Mon Sep 17 00:00:00 2001 From: fweber144 Date: Sun, 17 Nov 2024 16:15:05 +0100 Subject: [PATCH 06/40] tests: `NA`s in summaries should be `NA_real_`s now --- tests/testthat/helpers/testers.R | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/tests/testthat/helpers/testers.R b/tests/testthat/helpers/testers.R index 1998abadc..629db3ac2 100644 --- a/tests/testthat/helpers/testers.R +++ b/tests/testthat/helpers/testers.R @@ -2736,11 +2736,7 @@ smmry_ref_tester <- function( is_lat_kfold <- latent_expected && !resp_oscale_expected && identical(cv_method_expected, "kfold") - if (is_lat_kfold) { - expect_true(is.vector(smmry_ref, "logical"), info = info_str) - } else { - expect_true(is.vector(smmry_ref, "numeric"), info = info_str) - } + expect_true(is.vector(smmry_ref, "numeric"), info = info_str) smmry_nms <- character() stats_mean_name <- stats_expected smmry_nms <- c(smmry_nms, From 8564b4a6f333015596f58c672d718b0a6caf6d98 Mon Sep 17 00:00:00 2001 From: fweber144 Date: Sun, 17 Nov 2024 16:54:11 +0100 Subject: [PATCH 07/40] in this branch, `weighted.sd()` and `auc()` have been renamed --- tests/testthat/test_misc.R | 98 +++++++++++++++++++------------------- 1 file changed, 49 insertions(+), 49 deletions(-) diff --git a/tests/testthat/test_misc.R b/tests/testthat/test_misc.R index 982ac5e42..93284f123 100644 --- a/tests/testthat/test_misc.R +++ b/tests/testthat/test_misc.R @@ -236,142 +236,142 @@ test_that(".srs_diff_est_w() propagates input `NA`s to its output", { ) }) -test_that("weighted.sd() with `na.rm = FALSE` propagates input `NA`s to its output", { +test_that(".weighted_sd() with `na.rm = FALSE` propagates input `NA`s to its output", { expect_identical( - weighted.sd(x = rep(NA, nobsv), - w = rep(NA, nobsv)), + .weighted_sd(x = rep(NA, nobsv), + w = rep(NA, nobsv)), NA_real_, info = "all inputs are all-`NA`s" ) expect_identical( - weighted.sd(x = rep(-0.8, nobsv), - w = c(rep_len(c(2, 3), length.out = nobsv - 1L), NA)), + .weighted_sd(x = rep(-0.8, nobsv), + w = c(rep_len(c(2, 3), length.out = nobsv - 1L), NA)), NA_real_, info = "`x` without `NA`s, `w` with a single `NA`" ) expect_identical( - weighted.sd(x = c(rep(-0.8, nobsv - 1L), NA), - w = rep_len(c(2, 3), length.out = nobsv)), + .weighted_sd(x = c(rep(-0.8, nobsv - 1L), NA), + w = rep_len(c(2, 3), length.out = nobsv)), NA_real_, info = "`x` with a single `NA`, `w` without `NA`s" ) expect_identical( - weighted.sd(x = rep(-0.8, nobsv), - w = rep(NA, nobsv)), + .weighted_sd(x = rep(-0.8, nobsv), + w = rep(NA, nobsv)), NA_real_, info = "`x` without `NA`s, `w` all-`NA`s" ) expect_identical( - weighted.sd(x = rep(NA, nobsv), - w = rep_len(c(2, 3), length.out = nobsv)), + .weighted_sd(x = rep(NA, nobsv), + w = rep_len(c(2, 3), length.out = nobsv)), NA_real_, info = "`x` all-`NA`s, `w` without `NA`s" ) }) -test_that("auc() works", { +test_that(".auc() works", { nobsv_auc <- 19L expect_equal( - auc(cbind(rep_len(c(0, 1), length.out = nobsv_auc), - imperfect_alternation(c(0.3, 0.8), length.out = nobsv_auc), - rep(1, nobsv_auc))), + .auc(cbind(rep_len(c(0, 1), length.out = nobsv_auc), + imperfect_alternation(c(0.3, 0.8), length.out = nobsv_auc), + rep(1, nobsv_auc))), 0.6833333333333333333333, info = "`wobs` column only `1`s" ) expect_equal( - auc(cbind(rep_len(c(0, 1), length.out = nobsv_auc), - imperfect_alternation(c(0.3, 0.8), length.out = nobsv_auc), - imperfect_alternation(c(1, 2), n_tail = 11L, length.out = nobsv_auc))), + .auc(cbind(rep_len(c(0, 1), length.out = nobsv_auc), + imperfect_alternation(c(0.3, 0.8), length.out = nobsv_auc), + imperfect_alternation(c(1, 2), n_tail = 11L, length.out = nobsv_auc))), 0.717948717948718062587, info = "`wobs` column with imperfect alternation of `1`s and `2`s" ) }) -test_that("auc() propagates input `NA`s to its output", { +test_that(".auc() propagates input `NA`s to its output", { nobsv_auc <- 19L expect_identical( - auc(cbind(rep(NA, nobsv_auc), - rep(NA, nobsv_auc), - rep(NA, nobsv_auc))), + .auc(cbind(rep(NA, nobsv_auc), + rep(NA, nobsv_auc), + rep(NA, nobsv_auc))), NA_real_, info = "all inputs are all-`NA`s" ) expect_identical( - auc(cbind(rep(1, nobsv_auc), - rep(NA, nobsv_auc), - rep(1, nobsv_auc))), + .auc(cbind(rep(1, nobsv_auc), + rep(NA, nobsv_auc), + rep(1, nobsv_auc))), NA_real_, info = paste0("`resp` column without `NA`s, ", "`pred` column only `NA`s, ", "`wobs` column without `NA`s") ) expect_identical( - auc(cbind(rep(1, nobsv_auc), - rep(NA, nobsv_auc), - c(rep(1, nobsv_auc - 1L), NA))), + .auc(cbind(rep(1, nobsv_auc), + rep(NA, nobsv_auc), + c(rep(1, nobsv_auc - 1L), NA))), NA_real_, info = paste0("`resp` column without `NA`s, ", "`pred` column only `NA`s, ", "`wobs` column with a single `NA`") ) expect_identical( - auc(cbind(c(rep_len(c(0, 1), length.out = nobsv_auc - 1L), NA), - rep(0.7, nobsv_auc), - rep(1, nobsv_auc))), + .auc(cbind(c(rep_len(c(0, 1), length.out = nobsv_auc - 1L), NA), + rep(0.7, nobsv_auc), + rep(1, nobsv_auc))), NA_real_, info = paste0("`resp` column with a single `NA`, ", "`pred` column without `NA`s, ", "`wobs` column without `NA`s") ) expect_identical( - auc(cbind(rep(NA, nobsv_auc), - rep(0.7, nobsv_auc), - rep(1, nobsv_auc))), + .auc(cbind(rep(NA, nobsv_auc), + rep(0.7, nobsv_auc), + rep(1, nobsv_auc))), NA_real_, info = paste0("`resp` column only `NA`s, ", "`pred` column without `NA`s, ", "`wobs` column without `NA`s") ) expect_identical( - auc(cbind(c(rep_len(c(0, 1), length.out = nobsv_auc - 1L), NA), - c(rep(0.7, nobsv_auc - 1L), NA), - rep(1, nobsv_auc))), + .auc(cbind(c(rep_len(c(0, 1), length.out = nobsv_auc - 1L), NA), + c(rep(0.7, nobsv_auc - 1L), NA), + rep(1, nobsv_auc))), NA_real_, info = paste0("`resp` column with a single `NA`, ", "`pred` column with a single `NA`, ", "`wobs` column without `NA`s") ) expect_identical( - auc(cbind(c(rep_len(c(0, 1), length.out = nobsv_auc - 1L), NA), - c(NA, rep(0.7, nobsv_auc - 1L)), - rep(1, nobsv_auc))), + .auc(cbind(c(rep_len(c(0, 1), length.out = nobsv_auc - 1L), NA), + c(NA, rep(0.7, nobsv_auc - 1L)), + rep(1, nobsv_auc))), NA_real_, info = paste0("`resp` column with a single `NA`, ", "`pred` column with a single `NA` (at different position), ", "`wobs` column without `NA`s") ) expect_identical( - auc(cbind(c(rep_len(c(0, 1), length.out = nobsv_auc - 1L), NA), - rep(NA, nobsv_auc), - rep(1, nobsv_auc))), + .auc(cbind(c(rep_len(c(0, 1), length.out = nobsv_auc - 1L), NA), + rep(NA, nobsv_auc), + rep(1, nobsv_auc))), NA_real_, info = paste0("`resp` column with a single `NA`, ", "`pred` column only `NA`s, ", "`wobs` column without `NA`s") ) expect_identical( - auc(cbind(rep(NA, nobsv_auc), - rep(0.7, nobsv_auc), - c(rep(1, nobsv_auc - 1L), NA))), + .auc(cbind(rep(NA, nobsv_auc), + rep(0.7, nobsv_auc), + c(rep(1, nobsv_auc - 1L), NA))), NA_real_, info = paste0("`resp` column only `NA`s, ", "`pred` column without `NA`s, ", "`wobs` column with a single `NA`") ) expect_identical( - auc(cbind(rep_len(c(0, 1), length.out = nobsv_auc), - imperfect_alternation(c(0.3, 0.8), length.out = nobsv_auc), - rep(NA, nobsv_auc))), + .auc(cbind(rep_len(c(0, 1), length.out = nobsv_auc), + imperfect_alternation(c(0.3, 0.8), length.out = nobsv_auc), + rep(NA, nobsv_auc))), NA_real_, info = paste0("`resp` column without `NA`s, ", "`pred` column without `NA`s, ", From ebfc665e545431f4462c85d3ba7428f6b601203f Mon Sep 17 00:00:00 2001 From: fweber144 Date: Sun, 22 Dec 2024 20:09:07 +0100 Subject: [PATCH 08/40] divide by `(n_full - 1)` instead of `n_full` where necessary, see --- R/summary_funs.R | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/R/summary_funs.R b/R/summary_funs.R index 234b83ce6..2990173d5 100644 --- a/R/summary_funs.R +++ b/R/summary_funs.R @@ -348,10 +348,10 @@ get_stat <- function(summaries, summaries_baseline = NULL, ((mu_baseline - y)^2 - mse_b))[loo_inds], y_idx = loo_inds, wobs = wobs) - cov_mse_e_b <- srs_diffe$y_hat / n_full^2 + cov_mse_e_b <- srs_diffe$y_hat / (n_full * (n_full - 1)) } else { cov_mse_e_b <- mean(wobs * ((mu - y)^2 - mse_e) * - ((mu_baseline - y)^2 - mse_b)) / n_full + ((mu_baseline - y)^2 - mse_b)) / (n_full - 1) } value_se <- sqrt(value_se^2 - 2 * cov_mse_e_b + var_mse_b) } @@ -389,15 +389,15 @@ get_stat <- function(summaries, summaries_baseline = NULL, y_idx = loo_inds, wobs = wobs) } - cov_mse_e_y <- srs_diffe$y_hat / n_full^2 + cov_mse_e_y <- srs_diffe$y_hat / (n_full * (n_full - 1)) } else { if (is.null(summaries_baseline)) { cov_mse_e_y <- mean(wobs * ((mu - y)^2 - mse_e) * - ((mean(y) - y)^2 - mse_y)) / n_full + ((mean(y) - y)^2 - mse_y)) / (n_full - 1) } else { cov_mse_e_y <- mean(wobs * ((mu - y)^2 - mse_e - ((mu_baseline - y)^2 - mse_b)) * - ((mean(y) - y)^2 - mse_y)) / n_full + ((mean(y) - y)^2 - mse_y)) / (n_full - 1) } } # part of delta se comes automatically via mse From e841bcd1e74281a5e5439c88fe3b33ff95d97506 Mon Sep 17 00:00:00 2001 From: fweber144 Date: Sun, 22 Dec 2024 21:02:42 +0100 Subject: [PATCH 09/40] add test `".srs_diff_est_w() works as expected"` (copied from 'loo', see ) --- tests/testthat/test_misc.R | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/tests/testthat/test_misc.R b/tests/testthat/test_misc.R index 93284f123..ecbfbc2c8 100644 --- a/tests/testthat/test_misc.R +++ b/tests/testthat/test_misc.R @@ -189,6 +189,43 @@ test_that(paste( expect_equal(psdat$wobs, wobs_crr) }) +# Test copied from 'loo' package and changed function name, see +# +test_that(".srs_diff_est_w() works as expected", { + set.seed(1234) + N <- 1000 + y_true <- 1:N + sigma_hat_true <- sqrt(N * sum((y_true - mean(y_true))^2) / length(y_true)) + y_approx <- rnorm(N, y_true, 0.1) + m <- 100 + sigma_hat <- y_hat <- se_y_hat <- numeric(10000) + for(i in 1:10000){ + y_idx <- sample(1:N, size = m) + y <- y_true[y_idx] + res <- .srs_diff_est_w(y_approx, y, y_idx) + y_hat[i] <- res$y_hat + se_y_hat[i] <- sqrt(res$v_y_hat) + sigma_hat[i] <- sqrt(res$hat_v_y) + } + expect_equal(mean(y_hat), sum(y_true), tol = 0.1) + + in_ki <- y_hat + 2 * se_y_hat > sum(y_true) & y_hat - 2*se_y_hat < sum(y_true) + expect_equal(mean(in_ki), 0.95, tol = 0.01) + + # Should be unbiased + expect_equal(mean(sigma_hat), sigma_hat_true, tol = 0.1) + + m <- N + y_idx <- sample(1:N, size = m) + y <- y_true[y_idx] + res <- .srs_diff_est_w(y_approx, y, y_idx) + expect_equal(res$y_hat, 500500, tol = 0.0001) + expect_equal(res$v_y_hat, 0, tol = 0.0001) + expect_equal(sqrt(res$hat_v_y), sigma_hat_true, tol = 0.1) +}) + +# TODO: Add test for .srs_diff_est_w() with `wobs` not full of ones. + test_that(".srs_diff_est_w() propagates input `NA`s to its output", { nloo_tst <- nobsv %/% 5L loo_inds_tst <- ceiling(seq(1L, nobsv, length.out = nloo_tst)) From ae82ffcbb4565d6b51287b7fbf54aab3206b8021 Mon Sep 17 00:00:00 2001 From: fweber144 Date: Sat, 4 Jan 2025 21:00:29 +0100 Subject: [PATCH 10/40] explain "fast LOO" and "full LOO" in the docs, see --- R/cv_varsel.R | 14 ++++++++------ man/cv_varsel.Rd | 13 +++++++------ 2 files changed, 15 insertions(+), 12 deletions(-) diff --git a/R/cv_varsel.R b/R/cv_varsel.R index 40889d153..2f77e1504 100644 --- a/R/cv_varsel.R +++ b/R/cv_varsel.R @@ -22,12 +22,14 @@ #' performed. See also section "Note" below. #' @param nloo Only relevant if `cv_method = "LOO"` and `validate_search = #' TRUE`. If `nloo > 0` is smaller than the number of all observations, full -#' LOO-CV is approximated by combining the fast LOO result for the selected -#' models and `nloo` leave-one-out searches using the difference estimator -#' with simple random sampling (SRS) without replacement (WOR) (Magnusson et -#' al., 2020). Smaller values lead to faster computation, but higher -#' uncertainty in the evaluation part. If `NULL`, all observations are used -#' (as by default). +#' LOO-CV (i.e., PSIS-LOO CV with `validate_search = TRUE` and with `nloo = +#' n` where `n` denotes the number of all observations) is approximated by +#' combining the fast (i.e., `validate_search = FALSE`) LOO result for the +#' selected models and `nloo` leave-one-out searches using the difference +#' estimator with simple random sampling (SRS) without replacement (WOR) +#' (Magnusson et al., 2020). Smaller `nloo` values lead to faster computation, +#' but higher uncertainty in the evaluation part. If `NULL`, all observations +#' are used (as by default). #' @param K Only relevant if `cv_method = "kfold"` and if `cvfits` is `NULL` #' (which is the case for reference model objects created by #' [get_refmodel.stanreg()] or [brms::get_refmodel.brmsfit()]). Number of diff --git a/man/cv_varsel.Rd b/man/cv_varsel.Rd index 6a1c43f61..9e6bc0efd 100644 --- a/man/cv_varsel.Rd +++ b/man/cv_varsel.Rd @@ -68,12 +68,13 @@ contrast to a standard LOO-CV). In the \code{"kfold"} case, a \eqn{K}-fold-CV is performed. See also section "Note" below.} \item{nloo}{Only relevant if \code{cv_method = "LOO"} and \code{validate_search = TRUE}. If \code{nloo > 0} is smaller than the number of all observations, full -LOO-CV is approximated by combining the fast LOO result for the selected -models and \code{nloo} leave-one-out searches using the difference estimator -with simple random sampling (SRS) without replacement (WOR) (Magnusson et -al., 2020). Smaller values lead to faster computation, but higher -uncertainty in the evaluation part. If \code{NULL}, all observations are used -(as by default).} +LOO-CV (i.e., PSIS-LOO CV with \code{validate_search = TRUE} and with \code{nloo = n} where \code{n} denotes the number of all observations) is approximated by +combining the fast (i.e., \code{validate_search = FALSE}) LOO result for the +selected models and \code{nloo} leave-one-out searches using the difference +estimator with simple random sampling (SRS) without replacement (WOR) +(Magnusson et al., 2020). Smaller \code{nloo} values lead to faster computation, +but higher uncertainty in the evaluation part. If \code{NULL}, all observations +are used (as by default).} \item{K}{Only relevant if \code{cv_method = "kfold"} and if \code{cvfits} is \code{NULL} (which is the case for reference model objects created by From 006ff44cdcaaa7b9295b99160def8162326c825f Mon Sep 17 00:00:00 2001 From: fweber144 Date: Sun, 5 Jan 2025 21:44:42 +0100 Subject: [PATCH 11/40] make computation of `var_e_i` in `.srs_diff_est_w()` numerically more stable, see --- R/misc.R | 3 +++ R/summary_funs.R | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/R/misc.R b/R/misc.R index e09f95b91..20f7d75ad 100644 --- a/R/misc.R +++ b/R/misc.R @@ -15,6 +15,9 @@ nms_y_wobs_test <- function(wobs_nm = "wobs") { } .weighted_sd <- function(x, w, na.rm = FALSE) { + if (length(w) == 1 && length(x) > 1) { + w <- rep_len(w, length.out = length(x)) + } if (na.rm) { ind <- !is.na(w) & !is.na(x) n <- sum(ind) diff --git a/R/summary_funs.R b/R/summary_funs.R index fa30a0ecf..b275082bf 100644 --- a/R/summary_funs.R +++ b/R/summary_funs.R @@ -631,7 +631,7 @@ get_nfeat_baseline <- function(object, baseline, stat, ...) { # eq (7) est_list$y_hat <- t_pi_tilde + t_e # eq (8) - var_e_i <- m / (m - 1) * (mean(wobs_m * e_i^2) - mean(wobs_m * e_i)^2) + var_e_i <- .weighted_sd(e_i, w = wobs_m)^2 est_list$v_y_hat <- N^2 * (1 - m / N) * var_e_i / m # eq (9) first row second `+` should be `-` # Supplementary material eq (6) has this correct From 09487e3bd1867ce16c2e0bb4db368921eec32b3d Mon Sep 17 00:00:00 2001 From: fweber144 Date: Sun, 5 Jan 2025 21:54:27 +0100 Subject: [PATCH 12/40] fix a comment, see --- R/summary_funs.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/summary_funs.R b/R/summary_funs.R index b275082bf..287bdb150 100644 --- a/R/summary_funs.R +++ b/R/summary_funs.R @@ -328,7 +328,7 @@ get_stat <- function(summaries, summaries_baseline = NULL, if (!is.null(summaries_baseline)) { mu_baseline <- summaries_baseline$mu } - # Use normal approximation for mse and delta method for rmse and R2 + # Use "exact" standard error for mse and delta method for rmse and R2 if (n_loo < n_full) { # subsampling difference estimator (Magnusson et al., 2020) srs_diffe <- .srs_diff_est_w(y_approx = (summaries_fast$mu - y)^2, From aecaa4ee04f7ced5469278f27c1d3e3312f2151d Mon Sep 17 00:00:00 2001 From: fweber144 Date: Mon, 6 Jan 2025 15:28:30 +0100 Subject: [PATCH 13/40] minor cleaning --- R/summary_funs.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/summary_funs.R b/R/summary_funs.R index 287bdb150..cd6e128bf 100644 --- a/R/summary_funs.R +++ b/R/summary_funs.R @@ -550,7 +550,7 @@ get_stat <- function(summaries, summaries_baseline = NULL, } } - if (stat %in% c("mse","rmse") && is.null(summaries_baseline)) { + if (stat %in% c("mse", "rmse") && is.null(summaries_baseline)) { # Compute mean and variance in log scale by matching the variance of a # log-normal approximation # https://en.wikipedia.org/wiki/Log-normal_distribution#Arithmetic_moments From 35280fbbd48ed31b263fbec5cb84ae4f8e7486be Mon Sep 17 00:00:00 2001 From: fweber144 Date: Mon, 6 Jan 2025 15:40:09 +0100 Subject: [PATCH 14/40] fix a comment, see --- R/summary_funs.R | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/R/summary_funs.R b/R/summary_funs.R index cd6e128bf..2091305d4 100644 --- a/R/summary_funs.R +++ b/R/summary_funs.R @@ -346,7 +346,10 @@ get_stat <- function(summaries, summaries_baseline = NULL, # store for later calculations mse_e <- value if (!is.null(summaries_baseline)) { - # delta=TRUE, variance of difference of two normally distributed + # delta=TRUE, variance of difference of two normally distributed random + # variables (log-normally in case of MSE and RMSE, although the central + # limit theorem would ensure convergence -- probably slower, though -- to + # a normal distribution even for MSE and RMSE) mse_b <- mean(wobs * (mu_baseline - y)^2) var_mse_b <- .weighted_sd((mu_baseline - y)^2, wobs)^2 / n_full if (n_loo < n_full) { From 6dcc223b57416906f3ea0c9e92657ee21f14724e Mon Sep 17 00:00:00 2001 From: fweber144 Date: Mon, 6 Jan 2025 20:39:54 +0100 Subject: [PATCH 15/40] replace `mean(y)` with `mean(wobs * y)`, see --- R/summary_funs.R | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/R/summary_funs.R b/R/summary_funs.R index 2091305d4..c3ae668ad 100644 --- a/R/summary_funs.R +++ b/R/summary_funs.R @@ -377,28 +377,28 @@ get_stat <- function(summaries, summaries_baseline = NULL, value_se <- sqrt(value_se^2 / mse_e / 4) } else if (stat == "R2") { # simple transformation of mse - mse_y <- mean(wobs * (mean(y) - y)^2) + mse_y <- mean(wobs * (mean(wobs * y) - y)^2) value <- 1 - mse_e / mse_y - ifelse(is.null(summaries_baseline), 0, 1 - mse_b / mse_y) # the first-order Taylor approximation of the variance - var_mse_y <- .weighted_sd((mean(y) - y)^2, wobs)^2 / n_full + var_mse_y <- .weighted_sd((mean(wobs * y) - y)^2, wobs)^2 / n_full if (n_loo < n_full) { mse_e_fast <- mean(wobs * (summaries_fast$mu - y)^2) if (is.null(summaries_baseline)) { srs_diffe <- .srs_diff_est_w(y_approx = ((summaries_fast$mu - y)^2 - mse_e_fast) * - ((mean(y) - y)^2 - mse_y), + ((mean(wobs * y) - y)^2 - mse_y), y = (((mu - y)^2 - mse_e) * - ((mean(y) - y)^2 - mse_y))[loo_inds], + ((mean(wobs * y) - y)^2 - mse_y))[loo_inds], y_idx = loo_inds, wobs = wobs) } else { srs_diffe <- .srs_diff_est_w(y_approx = ((summaries_fast$mu - y)^2 - mse_e_fast - ((mu_baseline - y)^2 - mse_b)) * - ((mean(y) - y)^2 - mse_y), + ((mean(wobs * y) - y)^2 - mse_y), y = (((mu - y)^2 - mse_e - ((mu_baseline - y)^2 - mse_b)) * - ((mean(y) - y)^2 - mse_y))[loo_inds], + ((mean(wobs * y) - y)^2 - mse_y))[loo_inds], y_idx = loo_inds, wobs = wobs) } @@ -406,11 +406,11 @@ get_stat <- function(summaries, summaries_baseline = NULL, } else { if (is.null(summaries_baseline)) { cov_mse_e_y <- mean(wobs * ((mu - y)^2 - mse_e) * - ((mean(y) - y)^2 - mse_y)) / (n_full - 1) + ((mean(wobs * y) - y)^2 - mse_y)) / (n_full - 1) } else { cov_mse_e_y <- mean(wobs * ((mu - y)^2 - mse_e - ((mu_baseline - y)^2 - mse_b)) * - ((mean(y) - y)^2 - mse_y)) / (n_full - 1) + ((mean(wobs * y) - y)^2 - mse_y)) / (n_full - 1) } } # part of delta se comes automatically via mse From 1ee2bbd7387a4924c9bc8bdc477c933248146f5e Mon Sep 17 00:00:00 2001 From: fweber144 Date: Mon, 6 Jan 2025 20:42:10 +0100 Subject: [PATCH 16/40] avoid redundant computations by introducing object `y_mean_w` --- R/summary_funs.R | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/R/summary_funs.R b/R/summary_funs.R index c3ae668ad..6a2a2eed2 100644 --- a/R/summary_funs.R +++ b/R/summary_funs.R @@ -376,29 +376,30 @@ get_stat <- function(summaries, summaries_baseline = NULL, # the first-order Taylor approximation of the variance value_se <- sqrt(value_se^2 / mse_e / 4) } else if (stat == "R2") { + y_mean_w <- mean(wobs * y) # simple transformation of mse - mse_y <- mean(wobs * (mean(wobs * y) - y)^2) + mse_y <- mean(wobs * (y_mean_w - y)^2) value <- 1 - mse_e / mse_y - ifelse(is.null(summaries_baseline), 0, 1 - mse_b / mse_y) # the first-order Taylor approximation of the variance - var_mse_y <- .weighted_sd((mean(wobs * y) - y)^2, wobs)^2 / n_full + var_mse_y <- .weighted_sd((y_mean_w - y)^2, wobs)^2 / n_full if (n_loo < n_full) { mse_e_fast <- mean(wobs * (summaries_fast$mu - y)^2) if (is.null(summaries_baseline)) { srs_diffe <- .srs_diff_est_w(y_approx = ((summaries_fast$mu - y)^2 - mse_e_fast) * - ((mean(wobs * y) - y)^2 - mse_y), + ((y_mean_w - y)^2 - mse_y), y = (((mu - y)^2 - mse_e) * - ((mean(wobs * y) - y)^2 - mse_y))[loo_inds], + ((y_mean_w - y)^2 - mse_y))[loo_inds], y_idx = loo_inds, wobs = wobs) } else { srs_diffe <- .srs_diff_est_w(y_approx = ((summaries_fast$mu - y)^2 - mse_e_fast - ((mu_baseline - y)^2 - mse_b)) * - ((mean(wobs * y) - y)^2 - mse_y), + ((y_mean_w - y)^2 - mse_y), y = (((mu - y)^2 - mse_e - ((mu_baseline - y)^2 - mse_b)) * - ((mean(wobs * y) - y)^2 - mse_y))[loo_inds], + ((y_mean_w - y)^2 - mse_y))[loo_inds], y_idx = loo_inds, wobs = wobs) } @@ -406,11 +407,11 @@ get_stat <- function(summaries, summaries_baseline = NULL, } else { if (is.null(summaries_baseline)) { cov_mse_e_y <- mean(wobs * ((mu - y)^2 - mse_e) * - ((mean(wobs * y) - y)^2 - mse_y)) / (n_full - 1) + ((y_mean_w - y)^2 - mse_y)) / (n_full - 1) } else { cov_mse_e_y <- mean(wobs * ((mu - y)^2 - mse_e - ((mu_baseline - y)^2 - mse_b)) * - ((mean(wobs * y) - y)^2 - mse_y)) / (n_full - 1) + ((y_mean_w - y)^2 - mse_y)) / (n_full - 1) } } # part of delta se comes automatically via mse From df4bad0f4c627736fb532ecf767d071a0d1af834 Mon Sep 17 00:00:00 2001 From: fweber144 Date: Mon, 6 Jan 2025 20:59:03 +0100 Subject: [PATCH 17/40] use `log1p()` for numerical stability (in `get_stat()`), see --- R/summary_funs.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/summary_funs.R b/R/summary_funs.R index 6a2a2eed2..db5f500ba 100644 --- a/R/summary_funs.R +++ b/R/summary_funs.R @@ -559,7 +559,7 @@ get_stat <- function(summaries, summaries_baseline = NULL, # log-normal approximation # https://en.wikipedia.org/wiki/Log-normal_distribution#Arithmetic_moments mul <- log(value^2 / sqrt(value_se^2 + value^2)) - varl <- log(1 + value_se^2 / value^2) + varl <- log1p(value_se^2 / value^2) lq <- qnorm(alpha_half, mean = mul, sd = sqrt(varl)) uq <- qnorm(one_minus_alpha_half, mean = mul, sd = sqrt(varl)) # Go back to linear scale From 971d5b68481d56b8649349d178f42a15f3e8c996 Mon Sep 17 00:00:00 2001 From: fweber144 Date: Sun, 12 Jan 2025 10:39:01 +0100 Subject: [PATCH 18/40] drop quotes from "exact" in a comment, see --- R/summary_funs.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/summary_funs.R b/R/summary_funs.R index db5f500ba..3b01b3855 100644 --- a/R/summary_funs.R +++ b/R/summary_funs.R @@ -328,7 +328,7 @@ get_stat <- function(summaries, summaries_baseline = NULL, if (!is.null(summaries_baseline)) { mu_baseline <- summaries_baseline$mu } - # Use "exact" standard error for mse and delta method for rmse and R2 + # Use exact standard error for mse and delta method for rmse and R2 if (n_loo < n_full) { # subsampling difference estimator (Magnusson et al., 2020) srs_diffe <- .srs_diff_est_w(y_approx = (summaries_fast$mu - y)^2, From a2fed535070a7f66de38ba546b3c24fc9d41c4ab Mon Sep 17 00:00:00 2001 From: fweber144 Date: Sun, 12 Jan 2025 21:49:09 +0100 Subject: [PATCH 19/40] fix docs for RMSE and R2 --- R/methods.R | 10 +++++++--- man/plot.vsel.Rd | 10 +++++++--- man/summary.vsel.Rd | 10 +++++++--- 3 files changed, 21 insertions(+), 9 deletions(-) diff --git a/R/methods.R b/R/methods.R index 17ae387cb..9fea3697f 100644 --- a/R/methods.R +++ b/R/methods.R @@ -1173,8 +1173,11 @@ plot.vsel <- function( #' * `"mse"`: mean squared error (only available in the situations mentioned #' in section "Details" below). #' * `"rmse"`: root mean squared error (only available in the situations -#' mentioned in section "Details" below). For the corresponding standard error -#' and lower and upper confidence interval bounds, the delta method is used. +#' mentioned in section "Details" below). For the corresponding standard +#' error, the delta method is used. +#' * `"R2"`: R-squared, i.e., coefficient of determination (only available in +#' the situations mentioned in section "Details" below). For the corresponding +#' standard error, the delta method is used. #' * `"acc"` (or its alias, `"pctcorr"`): classification accuracy (only #' available in the situations mentioned in section "Details" below). By #' "classification accuracy", we mean the proportion of correctly classified @@ -1222,7 +1225,8 @@ plot.vsel <- function( #' and `seed` (see [set.seed()], but defaulting to `NA` so that [set.seed()] #' is not called within that function at all). #' -#' @details The `stats` options `"mse"` and `"rmse"` are only available for: +#' @details The `stats` options `"mse"`, `"rmse"`, and `"R2"` are only available +#' for: #' * the traditional projection, #' * the latent projection with `resp_oscale = FALSE`, #' * the latent projection with `resp_oscale = TRUE` in combination with diff --git a/man/plot.vsel.Rd b/man/plot.vsel.Rd index 3ee685bdc..ce96ee898 100644 --- a/man/plot.vsel.Rd +++ b/man/plot.vsel.Rd @@ -61,8 +61,11 @@ of the MLPD. \item \code{"mse"}: mean squared error (only available in the situations mentioned in section "Details" below). \item \code{"rmse"}: root mean squared error (only available in the situations -mentioned in section "Details" below). For the corresponding standard error -and lower and upper confidence interval bounds, the delta method is used. +mentioned in section "Details" below). For the corresponding standard +error, the delta method is used. +\item \code{"R2"}: R-squared, i.e., coefficient of determination (only available in +the situations mentioned in section "Details" below). For the corresponding +standard error, the delta method is used. \item \code{"acc"} (or its alias, \code{"pctcorr"}): classification accuracy (only available in the situations mentioned in section "Details" below). By "classification accuracy", we mean the proportion of correctly classified @@ -191,7 +194,8 @@ available; inferred from \code{\link[=cv_proportions]{cv_proportions()}}). For a see \code{\link[=summary.vsel]{summary.vsel()}} and \code{\link[=performances]{performances()}}. } \details{ -The \code{stats} options \code{"mse"} and \code{"rmse"} are only available for: +The \code{stats} options \code{"mse"}, \code{"rmse"}, and \code{"R2"} are only available +for: \itemize{ \item the traditional projection, \item the latent projection with \code{resp_oscale = FALSE}, diff --git a/man/summary.vsel.Rd b/man/summary.vsel.Rd index 5aaa8191a..61ec97e1d 100644 --- a/man/summary.vsel.Rd +++ b/man/summary.vsel.Rd @@ -51,8 +51,11 @@ of the MLPD. \item \code{"mse"}: mean squared error (only available in the situations mentioned in section "Details" below). \item \code{"rmse"}: root mean squared error (only available in the situations -mentioned in section "Details" below). For the corresponding standard error -and lower and upper confidence interval bounds, the delta method is used. +mentioned in section "Details" below). For the corresponding standard +error, the delta method is used. +\item \code{"R2"}: R-squared, i.e., coefficient of determination (only available in +the situations mentioned in section "Details" below). For the corresponding +standard error, the delta method is used. \item \code{"acc"} (or its alias, \code{"pctcorr"}): classification accuracy (only available in the situations mentioned in section "Details" below). By "classification accuracy", we mean the proportion of correctly classified @@ -125,7 +128,8 @@ results printed at the bottom of the output created by this \code{\link[=summary method, see \code{\link[=performances]{performances()}}. } \details{ -The \code{stats} options \code{"mse"} and \code{"rmse"} are only available for: +The \code{stats} options \code{"mse"}, \code{"rmse"}, and \code{"R2"} are only available +for: \itemize{ \item the traditional projection, \item the latent projection with \code{resp_oscale = FALSE}, From 21317229fdf49b2dd113986c1137b64780f73070 Mon Sep 17 00:00:00 2001 From: fweber144 Date: Sun, 12 Jan 2025 21:37:15 +0100 Subject: [PATCH 20/40] add `stat = "R2"` to the tests; also take into account that `stat = "rmse"` does not need a PRNG seed anymore (and should produce a point estimate between the CI bounds) --- tests/testthat/helpers/testers.R | 20 ++++++++++++-------- tests/testthat/setup.R | 13 ++++++------- tests/testthat/test_augdat.R | 8 ++++---- tests/testthat/test_methods_vsel.R | 10 +++++----- 4 files changed, 27 insertions(+), 24 deletions(-) diff --git a/tests/testthat/helpers/testers.R b/tests/testthat/helpers/testers.R index 629db3ac2..00b9fdb96 100644 --- a/tests/testthat/helpers/testers.R +++ b/tests/testthat/helpers/testers.R @@ -2670,8 +2670,9 @@ smmry_sub_tester <- function( if ("lower" %in% type_expected && !is_lat_kfold) { lower_nm <- paste(stats_expected, "lower", sep = ".") for (stat_idx in seq_along(stats_expected)) { - if (!stats_expected[stat_idx] %in% c("rmse", "auc")) { - # RMSE and AUC are excluded here because of PR #347. + if (!stats_expected[stat_idx] %in% c("auc")) { + # AUC is excluded here because of PR #347 (originally, RMSE was excluded + # as well, but PR #496 switched to the delta method for RMSE). expect_true(all(smmry_sub[, stats_mean_name[stat_idx]] >= smmry_sub[, lower_nm[stat_idx]]), info = info_str) @@ -2681,8 +2682,9 @@ smmry_sub_tester <- function( if ("upper" %in% type_expected && !is_lat_kfold) { upper_nm <- paste(stats_expected, "upper", sep = ".") for (stat_idx in seq_along(stats_expected)) { - if (!stats_expected[stat_idx] %in% c("rmse", "auc")) { - # RMSE and AUC are excluded here because of PR #347. + if (!stats_expected[stat_idx] %in% c("auc")) { + # AUC is excluded here because of PR #347 (originally, RMSE was excluded + # as well, but PR #496 switched to the delta method for RMSE). expect_true(all(smmry_sub[, stats_mean_name[stat_idx]] <= smmry_sub[, upper_nm[stat_idx]]), info = info_str) @@ -2758,8 +2760,9 @@ smmry_ref_tester <- function( if ("lower" %in% type_expected && !is_lat_kfold && !from_datafit) { lower_nm <- paste(stats_expected, "lower", sep = ".") for (stat_idx in seq_along(stats_expected)) { - if (!stats_expected[stat_idx] %in% c("rmse", "auc")) { - # RMSE and AUC are excluded here because of PR #347. + if (!stats_expected[stat_idx] %in% c("auc")) { + # AUC is excluded here because of PR #347 (originally, RMSE was excluded + # as well, but PR #496 switched to the delta method for RMSE). expect_true(all(smmry_ref[stats_mean_name[stat_idx]] >= smmry_ref[lower_nm[stat_idx]]), info = info_str) @@ -2769,8 +2772,9 @@ smmry_ref_tester <- function( if ("upper" %in% type_expected && !is_lat_kfold && !from_datafit) { upper_nm <- paste(stats_expected, "upper", sep = ".") for (stat_idx in seq_along(stats_expected)) { - if (!stats_expected[stat_idx] %in% c("rmse", "auc")) { - # RMSE and AUC are excluded here because of PR #347. + if (!stats_expected[stat_idx] %in% c("auc")) { + # AUC is excluded here because of PR #347 (originally, RMSE was excluded + # as well, but PR #496 switched to the delta method for RMSE). expect_true(all(smmry_ref[stats_mean_name[stat_idx]] <= smmry_ref[upper_nm[stat_idx]]), info = info_str) diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R index 733e6618e..623a33ba5 100755 --- a/tests/testthat/setup.R +++ b/tests/testthat/setup.R @@ -947,11 +947,10 @@ vsel_funs <- nlist("summary.vsel", "plot.vsel", "suggest_size.vsel") # projection (or the latent projection with `resp_oscale = FALSE` or the latent # projection with `resp_oscale = TRUE`, but the latter only in combination with # `$family$cats` being `NULL`): -stats_common <- c("elpd", "mlpd", "gmpd", "mse", "rmse") -# Performance statistics for the binomial() family only, when using the -# traditional projection (or the latent projection with `resp_oscale = TRUE`, -# but the latter only in combination with `$family$cats` being -# `NULL`): +stats_common <- c("elpd", "mlpd", "gmpd", "mse", "rmse", "R2") +# Performance statistics for the binomial() family when using the traditional +# projection (or the latent projection with `resp_oscale = TRUE` and +# `$family$cats` being `NULL`): stats_binom <- c(stats_common, "acc", "auc") # For creating test setups: stats_tst <- list( @@ -1700,7 +1699,7 @@ if (run_vs) { identical, 0L)) smmrys_vs <- lapply(args_smmry_vs, function(args_smmry_vs_i) { - if (any(c("rmse", "auc") %in% args_smmry_vs_i$stats)) { + if (any(c("auc") %in% args_smmry_vs_i$stats)) { smmry_seed <- list(seed = seed3_tst) } else { smmry_seed <- list() @@ -1723,7 +1722,7 @@ if (run_cvvs) { identical, 0L)) smmrys_cvvs <- lapply(args_smmry_cvvs, function(args_smmry_cvvs_i) { - if (any(c("rmse", "auc") %in% args_smmry_cvvs_i$stats)) { + if (any(c("auc") %in% args_smmry_cvvs_i$stats)) { smmry_seed <- list(seed = seed3_tst) } else { smmry_seed <- list() diff --git a/tests/testthat/test_augdat.R b/tests/testthat/test_augdat.R index a33831db1..95c46bf41 100644 --- a/tests/testthat/test_augdat.R +++ b/tests/testthat/test_augdat.R @@ -418,10 +418,10 @@ test_that(paste( # Exclude statistics which are not supported for the augmented-data # projection: smmry_vs_trad$perf_sub <- smmry_vs_trad$perf_sub[ - , -grep("mse|auc", names(smmry_vs_trad$perf_sub)), drop = FALSE + , -grep("mse|R2|auc", names(smmry_vs_trad$perf_sub)), drop = FALSE ] smmry_vs_trad$perf_ref <- smmry_vs_trad$perf_ref[ - -grep("mse|auc", names(smmry_vs_trad$perf_ref)) + -grep("mse|R2|auc", names(smmry_vs_trad$perf_ref)) ] expect_equal(smmry_vs$perf_sub, smmry_vs_trad$perf_sub, tolerance = 1e-6, info = tstsetup) @@ -592,10 +592,10 @@ test_that(paste( # Exclude statistics which are not supported for the augmented-data # projection: smmry_cvvs_trad$perf_sub <- smmry_cvvs_trad$perf_sub[ - , -grep("mse|auc", names(smmry_cvvs_trad$perf_sub)), drop = FALSE + , -grep("mse|R2|auc", names(smmry_cvvs_trad$perf_sub)), drop = FALSE ] smmry_cvvs_trad$perf_ref <- smmry_cvvs_trad$perf_ref[ - -grep("mse|auc", names(smmry_cvvs_trad$perf_ref)) + -grep("mse|R2|auc", names(smmry_cvvs_trad$perf_ref)) ] is_kfold <- identical( args_cvvs[[args_smmry_cvvs[[tstsetup]]$tstsetup_vsel]]$cv_method, diff --git a/tests/testthat/test_methods_vsel.R b/tests/testthat/test_methods_vsel.R index bef80b020..d5c0523dc 100644 --- a/tests/testthat/test_methods_vsel.R +++ b/tests/testthat/test_methods_vsel.R @@ -155,7 +155,7 @@ test_that("performances.vsel() is a shortcut", { skip_if_not(run_cvvs) for (tstsetup in names(smmrys_vs)) { args_smmry_i <- args_smmry_vs[[tstsetup]] - if (any(c("rmse", "auc") %in% args_smmry_i$stats)) { + if (any(c("auc") %in% args_smmry_i$stats)) { smmry_seed <- list(seed = seed3_tst) } else { smmry_seed <- list() @@ -169,7 +169,7 @@ test_that("performances.vsel() is a shortcut", { } for (tstsetup in names(smmrys_cvvs)) { args_smmry_i <- args_smmry_cvvs[[tstsetup]] - if (any(c("rmse", "auc") %in% args_smmry_i$stats)) { + if (any(c("auc") %in% args_smmry_i$stats)) { smmry_seed <- list(seed = seed3_tst) } else { smmry_seed <- list() @@ -239,7 +239,7 @@ test_that(paste( skip_if_not(run_vs) for (tstsetup in head(names(smmrys_vs), 1)) { args_smmry_vs_i <- args_smmry_vs[[tstsetup]] - if (any(c("rmse", "auc") %in% args_smmry_vs_i$stats)) { + if (any(c("auc") %in% args_smmry_vs_i$stats)) { smmry_seed <- list(seed = seed3_tst) } else { smmry_seed <- list() @@ -264,7 +264,7 @@ test_that(paste( skip_if_not(run_cvvs) for (tstsetup in head(names(smmrys_cvvs), 1)) { args_smmry_cvvs_i <- args_smmry_cvvs[[tstsetup]] - if (any(c("rmse", "auc") %in% args_smmry_cvvs_i$stats)) { + if (any(c("auc") %in% args_smmry_cvvs_i$stats)) { smmry_seed <- list(seed = seed3_tst) } else { smmry_seed <- list() @@ -527,7 +527,7 @@ test_that("`stat` works", { "common_stats")) stat_vec <- stats_tst[[stat_crr_nm]]$stats for (stat_crr in stat_vec) { - if (stat_crr %in% c("rmse", "auc")) { + if (stat_crr %in% c("auc")) { suggsize_seed <- seed3_tst } else { suggsize_seed <- NULL From 93c3f8ecdd92c7045a7002ef414df1d96a1fc8f6 Mon Sep 17 00:00:00 2001 From: fweber144 Date: Thu, 16 Jan 2025 15:39:42 +0100 Subject: [PATCH 21/40] set a negative squared standard error which is numerically equal to zero to a value of zero; this was necessary due to a failing `stat = "R2"` test for `datafit`s where such negative squared standard errors numerically equal to zero occurred --- R/summary_funs.R | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/R/summary_funs.R b/R/summary_funs.R index 3b01b3855..83591ac3d 100644 --- a/R/summary_funs.R +++ b/R/summary_funs.R @@ -419,9 +419,17 @@ get_stat <- function(summaries, summaries_baseline = NULL, # delta=TRUE mse_e <- mse_e - mse_b } - value_se <- sqrt((value_se^2 - - 2 * mse_e / mse_y * cov_mse_e_y + - (mse_e / mse_y)^2 * var_mse_y) / mse_y^2) + value_se_sq <- (value_se^2 - + 2 * mse_e / mse_y * cov_mse_e_y + + (mse_e / mse_y)^2 * var_mse_y) / mse_y^2 + if (!is.na(value_se_sq) && sign(value_se_sq) == -1) { + if (abs(value_se_sq) < sqrt(.Machine$double.eps)) { + value_se_sq <- 0 + } else { + stop("Negative (and numerically non-zero) `value_se_sq`.") + } + } + value_se <- sqrt(value_se_sq) } } else if (stat %in% c("acc", "pctcorr", "auc")) { y <- y_wobs_test$y From cc71813c131528258ca035f19fac72b07134e833 Mon Sep 17 00:00:00 2001 From: fweber144 Date: Sun, 12 Jan 2025 20:42:54 +0100 Subject: [PATCH 22/40] fix first-order Taylor approximation of the variance (delta method) for RMSE in the case `!is.null(summaries_baseline)`, see --- R/summary_funs.R | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/R/summary_funs.R b/R/summary_funs.R index 83591ac3d..205023039 100644 --- a/R/summary_funs.R +++ b/R/summary_funs.R @@ -366,7 +366,9 @@ get_stat <- function(summaries, summaries_baseline = NULL, cov_mse_e_b <- mean(wobs * ((mu - y)^2 - mse_e) * ((mu_baseline - y)^2 - mse_b)) / (n_full - 1) } - value_se <- sqrt(value_se^2 - 2 * cov_mse_e_b + var_mse_b) + if (stat != "rmse") { + value_se <- sqrt(value_se^2 - 2 * cov_mse_e_b + var_mse_b) + } } if (stat == "mse") { value <- mse_e - ifelse(is.null(summaries_baseline), 0, mse_b) @@ -374,7 +376,13 @@ get_stat <- function(summaries, summaries_baseline = NULL, # simple transformation of mse value <- sqrt(mse_e) - ifelse(is.null(summaries_baseline), 0, sqrt(mse_b)) # the first-order Taylor approximation of the variance - value_se <- sqrt(value_se^2 / mse_e / 4) + if (is.null(summaries_baseline)) { + value_se <- sqrt(value_se^2 / mse_e / 4) + } else { + value_se <- sqrt((value_se^2 / mse_e - + 2 * cov_mse_e_b / sqrt(mse_e * mse_b) + + var_mse_b / mse_b) / 4) + } } else if (stat == "R2") { y_mean_w <- mean(wobs * y) # simple transformation of mse From 26b36ab444548e693d41a76936fb89c660661986 Mon Sep 17 00:00:00 2001 From: fweber144 Date: Tue, 21 Jan 2025 21:35:15 +0100 Subject: [PATCH 23/40] change wording in a comment, see --- R/summary_funs.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/summary_funs.R b/R/summary_funs.R index 205023039..9224a021a 100644 --- a/R/summary_funs.R +++ b/R/summary_funs.R @@ -346,8 +346,8 @@ get_stat <- function(summaries, summaries_baseline = NULL, # store for later calculations mse_e <- value if (!is.null(summaries_baseline)) { - # delta=TRUE, variance of difference of two normally distributed random - # variables (log-normally in case of MSE and RMSE, although the central + # delta=TRUE, variance of difference of two normally distributed + # quantities (log-normally in case of MSE and RMSE, although the central # limit theorem would ensure convergence -- probably slower, though -- to # a normal distribution even for MSE and RMSE) mse_b <- mean(wobs * (mu_baseline - y)^2) From 4a3149975eda107ae157171d3c571ae8da1b7366 Mon Sep 17 00:00:00 2001 From: fweber144 Date: Fri, 31 Jan 2025 20:10:54 +0100 Subject: [PATCH 24/40] subsampled LOO and AUC: switch from warning to error, see --- R/summary_funs.R | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/R/summary_funs.R b/R/summary_funs.R index 9224a021a..9b7f8dfd6 100644 --- a/R/summary_funs.R +++ b/R/summary_funs.R @@ -93,11 +93,13 @@ weighted_summary_means <- function(y_wobs_test, family, wdraws, mu, dis, cl_ref, summaries_ref <- varsel$summaries$ref summaries_sub <- varsel$summaries$sub summaries_fast_sub <- varsel$summaries_fast$sub - if (!is.null(summaries_fast_sub)) { - if (any(stats %in% c("auc"))) { - warning("Subsampling LOO with AUC not implemented. Using fast LOO for ", - "submodel AUC.") - } + if (!is.null(summaries_fast_sub) && any(stats %in% c("auc"))) { + stop("Subsampled LOO-CV with AUC not implemented. Alternatives using ", + "`validate_search = TRUE` are full (i.e., non-subsampled) LOO-CV and ", + "K-fold CV. Otherwise, results from `validate_search = FALSE` (which ", + "often already exist at this point of the workflow) can be used, ", + "with the downside that the search part is not cross-validated in ", + "that case.") } if (!varsel$refmodel$family$for_latent && !resp_oscale) { From eda0cc1fb4720cadd83fd59ef11a5544033a5445 Mon Sep 17 00:00:00 2001 From: fweber144 Date: Fri, 31 Jan 2025 20:54:40 +0100 Subject: [PATCH 25/40] docs: mention that the `stats` option `"auc"` is not supported in case of subsampled LOO-CV, see --- R/cv_varsel.R | 18 ++++++++++-------- R/methods.R | 7 ++++++- man/cv_varsel.Rd | 16 +++++++++------- man/plot.vsel.Rd | 7 ++++++- man/summary.vsel.Rd | 7 ++++++- 5 files changed, 37 insertions(+), 18 deletions(-) diff --git a/R/cv_varsel.R b/R/cv_varsel.R index 2f77e1504..ddd84d179 100644 --- a/R/cv_varsel.R +++ b/R/cv_varsel.R @@ -22,14 +22,16 @@ #' performed. See also section "Note" below. #' @param nloo Only relevant if `cv_method = "LOO"` and `validate_search = #' TRUE`. If `nloo > 0` is smaller than the number of all observations, full -#' LOO-CV (i.e., PSIS-LOO CV with `validate_search = TRUE` and with `nloo = -#' n` where `n` denotes the number of all observations) is approximated by -#' combining the fast (i.e., `validate_search = FALSE`) LOO result for the -#' selected models and `nloo` leave-one-out searches using the difference -#' estimator with simple random sampling (SRS) without replacement (WOR) -#' (Magnusson et al., 2020). Smaller `nloo` values lead to faster computation, -#' but higher uncertainty in the evaluation part. If `NULL`, all observations -#' are used (as by default). +#' LOO-CV (i.e., PSIS-LOO CV with `validate_search = TRUE` and with `nloo = n` +#' where `n` denotes the number of all observations) is approximated by +#' subsampled LOO-CV, i.e., by combining the fast (i.e., `validate_search = +#' FALSE`) LOO result for the selected models and `nloo` leave-one-out +#' searches using the difference estimator with simple random sampling (SRS) +#' without replacement (WOR) (Magnusson et al., 2020). Smaller `nloo` values +#' lead to faster computation, but higher uncertainty in the evaluation part. +#' If `NULL`, all observations are used (as by default). Note that performance +#' statistic `"auc"` (see argument `stats` of [summary.vsel()] and +#' [plot.vsel()]) is not supported in case of subsampled LOO-CV. #' @param K Only relevant if `cv_method = "kfold"` and if `cvfits` is `NULL` #' (which is the case for reference model objects created by #' [get_refmodel.stanreg()] or [brms::get_refmodel.brmsfit()]). Number of diff --git a/R/methods.R b/R/methods.R index 9fea3697f..fc6ceea25 100644 --- a/R/methods.R +++ b/R/methods.R @@ -1186,7 +1186,9 @@ plot.vsel <- function( #' ("classification") for an observation. #' * `"auc"`: area under the ROC curve (only available in the situations #' mentioned in section "Details" below). For the corresponding standard error -#' and lower and upper confidence interval bounds, bootstrapping is used. +#' and lower and upper confidence interval bounds, bootstrapping is used. Not +#' supported in case of subsampled LOO-CV (see argument `nloo` of +#' [cv_varsel()]). #' @param type One or more items from `"mean"`, `"se"`, `"lower"`, `"upper"`, #' `"diff"`, and `"diff.se"` indicating which of these to compute for each #' item from `stats` (mean, standard error, lower and upper confidence @@ -1248,6 +1250,9 @@ plot.vsel <- function( #' latent projection with `resp_oscale = TRUE` in combination with #' `$family$cats` being `NULL`. #' +#' Note that the `stats` option `"auc"` is not supported in case of subsampled +#' LOO-CV (see argument `nloo` of [cv_varsel()]). +#' #' @return An object of class `vselsummary`. The elements of this object are not #' meant to be accessed directly but instead via helper functions #' ([print.vselsummary()] and [performances.vselsummary()]). diff --git a/man/cv_varsel.Rd b/man/cv_varsel.Rd index 9e6bc0efd..b60089d6e 100644 --- a/man/cv_varsel.Rd +++ b/man/cv_varsel.Rd @@ -68,13 +68,15 @@ contrast to a standard LOO-CV). In the \code{"kfold"} case, a \eqn{K}-fold-CV is performed. See also section "Note" below.} \item{nloo}{Only relevant if \code{cv_method = "LOO"} and \code{validate_search = TRUE}. If \code{nloo > 0} is smaller than the number of all observations, full -LOO-CV (i.e., PSIS-LOO CV with \code{validate_search = TRUE} and with \code{nloo = n} where \code{n} denotes the number of all observations) is approximated by -combining the fast (i.e., \code{validate_search = FALSE}) LOO result for the -selected models and \code{nloo} leave-one-out searches using the difference -estimator with simple random sampling (SRS) without replacement (WOR) -(Magnusson et al., 2020). Smaller \code{nloo} values lead to faster computation, -but higher uncertainty in the evaluation part. If \code{NULL}, all observations -are used (as by default).} +LOO-CV (i.e., PSIS-LOO CV with \code{validate_search = TRUE} and with \code{nloo = n} +where \code{n} denotes the number of all observations) is approximated by +subsampled LOO-CV, i.e., by combining the fast (i.e., \code{validate_search = FALSE}) LOO result for the selected models and \code{nloo} leave-one-out +searches using the difference estimator with simple random sampling (SRS) +without replacement (WOR) (Magnusson et al., 2020). Smaller \code{nloo} values +lead to faster computation, but higher uncertainty in the evaluation part. +If \code{NULL}, all observations are used (as by default). Note that performance +statistic \code{"auc"} (see argument \code{stats} of \code{\link[=summary.vsel]{summary.vsel()}} and +\code{\link[=plot.vsel]{plot.vsel()}}) is not supported in case of subsampled LOO-CV.} \item{K}{Only relevant if \code{cv_method = "kfold"} and if \code{cvfits} is \code{NULL} (which is the case for reference model objects created by diff --git a/man/plot.vsel.Rd b/man/plot.vsel.Rd index ce96ee898..6edc4fbd2 100644 --- a/man/plot.vsel.Rd +++ b/man/plot.vsel.Rd @@ -74,7 +74,9 @@ probability (the probabilities are model-based) is taken as the prediction ("classification") for an observation. \item \code{"auc"}: area under the ROC curve (only available in the situations mentioned in section "Details" below). For the corresponding standard error -and lower and upper confidence interval bounds, bootstrapping is used. +and lower and upper confidence interval bounds, bootstrapping is used. Not +supported in case of subsampled LOO-CV (see argument \code{nloo} of +\code{\link[=cv_varsel]{cv_varsel()}}). }} \item{deltas}{If \code{TRUE}, the submodel statistics are estimated relatively to @@ -222,6 +224,9 @@ The \code{stats} option \code{"auc"} is only available for: latent projection with \code{resp_oscale = TRUE} in combination with \verb{$family$cats} being \code{NULL}. } + +Note that the \code{stats} option \code{"auc"} is not supported in case of subsampled +LOO-CV (see argument \code{nloo} of \code{\link[=cv_varsel]{cv_varsel()}}). } \section{Horizontal lines}{ As long as the reference model's performance is computable, it is always diff --git a/man/summary.vsel.Rd b/man/summary.vsel.Rd index 61ec97e1d..5db51d0a9 100644 --- a/man/summary.vsel.Rd +++ b/man/summary.vsel.Rd @@ -64,7 +64,9 @@ probability (the probabilities are model-based) is taken as the prediction ("classification") for an observation. \item \code{"auc"}: area under the ROC curve (only available in the situations mentioned in section "Details" below). For the corresponding standard error -and lower and upper confidence interval bounds, bootstrapping is used. +and lower and upper confidence interval bounds, bootstrapping is used. Not +supported in case of subsampled LOO-CV (see argument \code{nloo} of +\code{\link[=cv_varsel]{cv_varsel()}}). }} \item{type}{One or more items from \code{"mean"}, \code{"se"}, \code{"lower"}, \code{"upper"}, @@ -156,6 +158,9 @@ The \code{stats} option \code{"auc"} is only available for: latent projection with \code{resp_oscale = TRUE} in combination with \verb{$family$cats} being \code{NULL}. } + +Note that the \code{stats} option \code{"auc"} is not supported in case of subsampled +LOO-CV (see argument \code{nloo} of \code{\link[=cv_varsel]{cv_varsel()}}). } \examples{ \dontshow{if (requireNamespace("rstanarm", quietly = TRUE)) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} From 95ee8cd3c41a988d4495a02acd0d2d35b7f07bcd Mon Sep 17 00:00:00 2001 From: fweber144 Date: Fri, 31 Jan 2025 21:21:37 +0100 Subject: [PATCH 26/40] docs: mention log-normal approximation for MSE and RMSE, see --- R/methods.R | 46 ++++++++++++++++++++++----------------------- man/plot.vsel.Rd | 13 ++++++++----- man/suggest_size.Rd | 26 ++++++++++++------------- man/summary.vsel.Rd | 19 ++++++++++--------- 4 files changed, 54 insertions(+), 50 deletions(-) diff --git a/R/methods.R b/R/methods.R index fc6ceea25..393d56701 100644 --- a/R/methods.R +++ b/R/methods.R @@ -1171,10 +1171,14 @@ plot.vsel <- function( #' confidence interval bounds are the exponentiated confidence interval bounds #' of the MLPD. #' * `"mse"`: mean squared error (only available in the situations mentioned -#' in section "Details" below). +#' in section "Details" below). For the corresponding confidence interval, a +#' log-normal approximation is used if `deltas` is `FALSE` and a normal +#' approximation is used if `deltas` is `TRUE`. #' * `"rmse"`: root mean squared error (only available in the situations #' mentioned in section "Details" below). For the corresponding standard -#' error, the delta method is used. +#' error, the delta method is used. For the corresponding confidence interval, +#' a log-normal approximation is used if `deltas` is `FALSE` and a normal +#' approximation is used if `deltas` is `TRUE`. #' * `"R2"`: R-squared, i.e., coefficient of determination (only available in #' the situations mentioned in section "Details" below). For the corresponding #' standard error, the delta method is used. @@ -1196,10 +1200,8 @@ plot.vsel <- function( #' reference model, and standard error of this difference, respectively; note #' that for the GMPD, `"diff"`, and `"diff.se"` actually refer to the ratio #' vs. the reference model, not the difference). The confidence interval -#' bounds belong to normal-approximation (or bootstrap or exponentiated -#' normal-approximation; see argument `stats`) confidence intervals with -#' (nominal) coverage `1 - alpha`. Items `"diff"` and `"diff.se"` are only -#' supported if `deltas` is `FALSE`. +#' bounds belong to confidence intervals with (nominal) coverage `1 - alpha`. +#' Items `"diff"` and `"diff.se"` are only supported if `deltas` is `FALSE`. #' @param deltas If `TRUE`, the submodel statistics are estimated relatively to #' the baseline model (see argument `baseline`). For the GMPD, the term #' "relatively" refers to the ratio vs. the baseline model (i.e., the submodel @@ -1207,9 +1209,8 @@ plot.vsel <- function( #' "relatively" refers to the difference from the baseline model (i.e., the #' submodel statistic minus the baseline model statistic). #' @param alpha A number determining the (nominal) coverage `1 - alpha` of the -#' normal-approximation (or bootstrap or exponentiated normal-approximation; -#' see argument `stats`) confidence intervals. For example, in case of the -#' normal approximation, `alpha = 2 * pnorm(-1)` corresponds to a confidence +#' confidence intervals. For example, in case of a normal-approximation +#' confidence interval, `alpha = 2 * pnorm(-1)` corresponds to a confidence #' interval stretching by one standard error on either side of the point #' estimate. #' @param baseline For [summary.vsel()]: Only relevant if `deltas` is `TRUE`. @@ -1551,7 +1552,8 @@ print.vsel <- function(x, digits = getOption("projpred.digits", 2), ...) { #' @param object An object of class `vsel` (returned by [varsel()] or #' [cv_varsel()]). #' @param stat Performance statistic (i.e., utility or loss) used for the -#' decision. See argument `stats` of [summary.vsel()] for possible choices. +#' decision. See argument `stats` of [summary.vsel()] and [plot.vsel()] for +#' possible choices. #' @param pct A number giving the proportion (*not* percents) of the *relative* #' null model utility one is willing to sacrifice. See section "Details" below #' for more information. @@ -1575,13 +1577,11 @@ print.vsel <- function(x, digits = getOption("projpred.digits", 2), ...) { #' @details In general (beware of special cases below), the suggested model #' size is the smallest model size \eqn{j \in \{0, 1, ..., #' \texttt{nterms\_max}\}}{{j = 0, 1, ..., nterms_max}} for which either the -#' lower or upper bound (depending on argument `type`) of the -#' normal-approximation (or bootstrap or exponentiated normal-approximation; -#' see argument `stat`) confidence interval (with nominal coverage `1 - -#' alpha`; see argument `alpha` of [summary.vsel()]) for \eqn{U_j - -#' U_{\mathrm{base}}}{U_j - U_base} (with \eqn{U_j} denoting the \eqn{j}-th -#' submodel's true utility and \eqn{U_{\mathrm{base}}}{U_base} denoting the -#' baseline model's true utility) +#' lower or upper bound (depending on argument `type`) of the confidence +#' interval (with nominal coverage `1 - alpha`; see argument `alpha` of +#' [summary.vsel()]) for \eqn{U_j - U_{\mathrm{base}}}{U_j - U_base} (with +#' \eqn{U_j} denoting the \eqn{j}-th submodel's true utility and +#' \eqn{U_{\mathrm{base}}}{U_base} denoting the baseline model's true utility) #' falls above (or is equal to) \deqn{\texttt{pct} \cdot (u_0 - #' u_{\mathrm{base}})}{pct * (u_0 - u_base)} where \eqn{u_0} denotes the null #' model's estimated utility and \eqn{u_{\mathrm{base}}}{u_base} the baseline @@ -1631,15 +1631,15 @@ print.vsel <- function(x, digits = getOption("projpred.digits", 2), ...) { #' `!is.na(thres_elpd)` with `stat %in% c("elpd", "mlpd", "gmpd")`), `alpha = #' 2 * pnorm(-1)`, `pct = 0`, and `type = "upper"` means that we select the #' smallest model size for which the upper bound of the `1 - 2 * pnorm(-1)` -#' (approximately 68.3%) confidence interval for \eqn{U_j - +#' (approximately 68.3 %) confidence interval for \eqn{U_j - #' U_{\mathrm{base}}}{U_j - U_base} #' (\eqn{\frac{U^\ast_j}{U^\ast_{\mathrm{base}}}}{U^*_j / U^*_base} in case of #' the GMPD) exceeds (or is equal to) zero (one in case of the GMPD), that is -#' (if `stat` is a performance statistic for which the normal approximation is -#' used, not the bootstrap and not the exponentiated normal approximation), -#' for which the submodel's utility estimate is at most one standard error -#' smaller than the baseline model's utility estimate (with that standard -#' error referring to the utility *difference*). +#' (if `stat` is a performance statistic for which a normal-approximation +#' confidence interval is used, see argument `stats` of [summary.vsel()] and +#' [plot.vsel()]), for which the submodel's utility estimate is at most one +#' standard error smaller than the baseline model's utility estimate (with +#' that standard error referring to the utility *difference*). #' #' Apart from the two [summary.vsel()] arguments mentioned above (`alpha` and #' `baseline`), `resp_oscale` is another important [summary.vsel()] argument diff --git a/man/plot.vsel.Rd b/man/plot.vsel.Rd index 6edc4fbd2..69d318d96 100644 --- a/man/plot.vsel.Rd +++ b/man/plot.vsel.Rd @@ -59,10 +59,14 @@ interval type is "exponentiated normal approximation" because the confidence interval bounds are the exponentiated confidence interval bounds of the MLPD. \item \code{"mse"}: mean squared error (only available in the situations mentioned -in section "Details" below). +in section "Details" below). For the corresponding confidence interval, a +log-normal approximation is used if \code{deltas} is \code{FALSE} and a normal +approximation is used if \code{deltas} is \code{TRUE}. \item \code{"rmse"}: root mean squared error (only available in the situations mentioned in section "Details" below). For the corresponding standard -error, the delta method is used. +error, the delta method is used. For the corresponding confidence interval, +a log-normal approximation is used if \code{deltas} is \code{FALSE} and a normal +approximation is used if \code{deltas} is \code{TRUE}. \item \code{"R2"}: R-squared, i.e., coefficient of determination (only available in the situations mentioned in section "Details" below). For the corresponding standard error, the delta method is used. @@ -87,9 +91,8 @@ statistic divided by the baseline model statistic). For all other \code{stats}, submodel statistic minus the baseline model statistic).} \item{alpha}{A number determining the (nominal) coverage \code{1 - alpha} of the -normal-approximation (or bootstrap or exponentiated normal-approximation; -see argument \code{stats}) confidence intervals. For example, in case of the -normal approximation, \code{alpha = 2 * pnorm(-1)} corresponds to a confidence +confidence intervals. For example, in case of a normal-approximation +confidence interval, \code{alpha = 2 * pnorm(-1)} corresponds to a confidence interval stretching by one standard error on either side of the point estimate.} diff --git a/man/suggest_size.Rd b/man/suggest_size.Rd index 3353a2358..460b5801e 100644 --- a/man/suggest_size.Rd +++ b/man/suggest_size.Rd @@ -27,7 +27,8 @@ See section "Details" below for some important arguments which may be passed here.} \item{stat}{Performance statistic (i.e., utility or loss) used for the -decision. See argument \code{stats} of \code{\link[=summary.vsel]{summary.vsel()}} for possible choices.} +decision. See argument \code{stats} of \code{\link[=summary.vsel]{summary.vsel()}} and \code{\link[=plot.vsel]{plot.vsel()}} for +possible choices.} \item{pct}{A number giving the proportion (\emph{not} percents) of the \emph{relative} null model utility one is willing to sacrifice. See section "Details" below @@ -67,12 +68,11 @@ decision based on what is most appropriate for the problem at hand. In general (beware of special cases below), the suggested model size is the smallest model size \eqn{j \in \{0, 1, ..., \texttt{nterms\_max}\}}{{j = 0, 1, ..., nterms_max}} for which either the -lower or upper bound (depending on argument \code{type}) of the -normal-approximation (or bootstrap or exponentiated normal-approximation; -see argument \code{stat}) confidence interval (with nominal coverage \code{1 - alpha}; see argument \code{alpha} of \code{\link[=summary.vsel]{summary.vsel()}}) for \eqn{U_j - - U_{\mathrm{base}}}{U_j - U_base} (with \eqn{U_j} denoting the \eqn{j}-th -submodel's true utility and \eqn{U_{\mathrm{base}}}{U_base} denoting the -baseline model's true utility) +lower or upper bound (depending on argument \code{type}) of the confidence +interval (with nominal coverage \code{1 - alpha}; see argument \code{alpha} of +\code{\link[=summary.vsel]{summary.vsel()}}) for \eqn{U_j - U_{\mathrm{base}}}{U_j - U_base} (with +\eqn{U_j} denoting the \eqn{j}-th submodel's true utility and +\eqn{U_{\mathrm{base}}}{U_base} denoting the baseline model's true utility) falls above (or is equal to) \deqn{\texttt{pct} \cdot (u_0 - u_{\mathrm{base}})}{pct * (u_0 - u_base)} where \eqn{u_0} denotes the null model's estimated utility and \eqn{u_{\mathrm{base}}}{u_base} the baseline @@ -120,15 +120,15 @@ above \emph{or} \eqn{\frac{u^\ast_j}{u^\ast_{\mathrm{base}}} > For example (disregarding the special extensions in case of \code{!is.na(thres_elpd)} with \code{stat \%in\% c("elpd", "mlpd", "gmpd")}), \code{alpha = 2 * pnorm(-1)}, \code{pct = 0}, and \code{type = "upper"} means that we select the smallest model size for which the upper bound of the \code{1 - 2 * pnorm(-1)} -(approximately 68.3\%) confidence interval for \eqn{U_j - +(approximately 68.3 \%) confidence interval for \eqn{U_j - U_{\mathrm{base}}}{U_j - U_base} (\eqn{\frac{U^\ast_j}{U^\ast_{\mathrm{base}}}}{U^*_j / U^*_base} in case of the GMPD) exceeds (or is equal to) zero (one in case of the GMPD), that is -(if \code{stat} is a performance statistic for which the normal approximation is -used, not the bootstrap and not the exponentiated normal approximation), -for which the submodel's utility estimate is at most one standard error -smaller than the baseline model's utility estimate (with that standard -error referring to the utility \emph{difference}). +(if \code{stat} is a performance statistic for which a normal-approximation +confidence interval is used, see argument \code{stats} of \code{\link[=summary.vsel]{summary.vsel()}} and +\code{\link[=plot.vsel]{plot.vsel()}}), for which the submodel's utility estimate is at most one +standard error smaller than the baseline model's utility estimate (with +that standard error referring to the utility \emph{difference}). Apart from the two \code{\link[=summary.vsel]{summary.vsel()}} arguments mentioned above (\code{alpha} and \code{baseline}), \code{resp_oscale} is another important \code{\link[=summary.vsel]{summary.vsel()}} argument diff --git a/man/summary.vsel.Rd b/man/summary.vsel.Rd index 5db51d0a9..5534f4142 100644 --- a/man/summary.vsel.Rd +++ b/man/summary.vsel.Rd @@ -49,10 +49,14 @@ interval type is "exponentiated normal approximation" because the confidence interval bounds are the exponentiated confidence interval bounds of the MLPD. \item \code{"mse"}: mean squared error (only available in the situations mentioned -in section "Details" below). +in section "Details" below). For the corresponding confidence interval, a +log-normal approximation is used if \code{deltas} is \code{FALSE} and a normal +approximation is used if \code{deltas} is \code{TRUE}. \item \code{"rmse"}: root mean squared error (only available in the situations mentioned in section "Details" below). For the corresponding standard -error, the delta method is used. +error, the delta method is used. For the corresponding confidence interval, +a log-normal approximation is used if \code{deltas} is \code{FALSE} and a normal +approximation is used if \code{deltas} is \code{TRUE}. \item \code{"R2"}: R-squared, i.e., coefficient of determination (only available in the situations mentioned in section "Details" below). For the corresponding standard error, the delta method is used. @@ -76,10 +80,8 @@ interval bounds, mean difference to the corresponding statistic of the reference model, and standard error of this difference, respectively; note that for the GMPD, \code{"diff"}, and \code{"diff.se"} actually refer to the ratio vs. the reference model, not the difference). The confidence interval -bounds belong to normal-approximation (or bootstrap or exponentiated -normal-approximation; see argument \code{stats}) confidence intervals with -(nominal) coverage \code{1 - alpha}. Items \code{"diff"} and \code{"diff.se"} are only -supported if \code{deltas} is \code{FALSE}.} +bounds belong to confidence intervals with (nominal) coverage \code{1 - alpha}. +Items \code{"diff"} and \code{"diff.se"} are only supported if \code{deltas} is \code{FALSE}.} \item{deltas}{If \code{TRUE}, the submodel statistics are estimated relatively to the baseline model (see argument \code{baseline}). For the GMPD, the term @@ -89,9 +91,8 @@ statistic divided by the baseline model statistic). For all other \code{stats}, submodel statistic minus the baseline model statistic).} \item{alpha}{A number determining the (nominal) coverage \code{1 - alpha} of the -normal-approximation (or bootstrap or exponentiated normal-approximation; -see argument \code{stats}) confidence intervals. For example, in case of the -normal approximation, \code{alpha = 2 * pnorm(-1)} corresponds to a confidence +confidence intervals. For example, in case of a normal-approximation +confidence interval, \code{alpha = 2 * pnorm(-1)} corresponds to a confidence interval stretching by one standard error on either side of the point estimate.} From 6c6a596154d7963c709a08b9d3b57052b6ce9e83 Mon Sep 17 00:00:00 2001 From: fweber144 Date: Fri, 31 Jan 2025 21:30:44 +0100 Subject: [PATCH 27/40] `plot.vsel()`: mention log-normal approximation for MSE and RMSE, see --- R/methods.R | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/R/methods.R b/R/methods.R index 393d56701..d4f76a406 100644 --- a/R/methods.R +++ b/R/methods.R @@ -1069,7 +1069,10 @@ plot.vsel <- function( ci_type <- "bootstrap " } else if (all(stats %in% c("gmpd"))) { ci_type <- "exponentiated normal-approximation " - } else if (all(!stats %in% c("auc", "gmpd"))) { + } else if (all(stats %in% c("mse", "rmse")) && !deltas) { + ci_type <- "log-normal-approximation " + } else if (all(!stats %in% c("auc", "gmpd", "mse", "rmse")) || + (all(!stats %in% c("auc", "gmpd")) && deltas)) { ci_type <- "normal-approximation " } else { ci_type <- "" From 25a3ba7a0b8a562c0d03bf9218d81446354630bf Mon Sep 17 00:00:00 2001 From: fweber144 Date: Sat, 1 Feb 2025 06:45:24 +0100 Subject: [PATCH 28/40] docs: `baseline = "best"` is not supported in case of subsampled LOO-CV --- R/cv_varsel.R | 4 +++- R/methods.R | 3 ++- man/cv_varsel.Rd | 4 +++- man/plot.vsel.Rd | 2 +- man/summary.vsel.Rd | 2 +- 5 files changed, 10 insertions(+), 5 deletions(-) diff --git a/R/cv_varsel.R b/R/cv_varsel.R index ddd84d179..05e4e203f 100644 --- a/R/cv_varsel.R +++ b/R/cv_varsel.R @@ -31,7 +31,9 @@ #' lead to faster computation, but higher uncertainty in the evaluation part. #' If `NULL`, all observations are used (as by default). Note that performance #' statistic `"auc"` (see argument `stats` of [summary.vsel()] and -#' [plot.vsel()]) is not supported in case of subsampled LOO-CV. +#' [plot.vsel()]) is not supported in case of subsampled LOO-CV. Furthermore, +#' option `"best"` for argument `baseline` of [summary.vsel()] and +#' [plot.vsel()] is not supported in case of subsampled LOO-CV. #' @param K Only relevant if `cv_method = "kfold"` and if `cvfits` is `NULL` #' (which is the case for reference model objects created by #' [get_refmodel.stanreg()] or [brms::get_refmodel.brmsfit()]). Number of diff --git a/R/methods.R b/R/methods.R index d4f76a406..2ea886961 100644 --- a/R/methods.R +++ b/R/methods.R @@ -1219,7 +1219,8 @@ plot.vsel <- function( #' @param baseline For [summary.vsel()]: Only relevant if `deltas` is `TRUE`. #' For [plot.vsel()]: Always relevant. Either `"ref"` or `"best"`, indicating #' whether the baseline is the reference model or the best submodel found (in -#' terms of `stats[1]`), respectively. +#' terms of `stats[1]`), respectively. In case of subsampled LOO-CV, `baseline +#' = "best"` is not supported. #' @param resp_oscale Only relevant for the latent projection. A single logical #' value indicating whether to calculate the performance statistics on the #' original response scale (`TRUE`) or on latent scale (`FALSE`). diff --git a/man/cv_varsel.Rd b/man/cv_varsel.Rd index b60089d6e..dae0d10cc 100644 --- a/man/cv_varsel.Rd +++ b/man/cv_varsel.Rd @@ -76,7 +76,9 @@ without replacement (WOR) (Magnusson et al., 2020). Smaller \code{nloo} values lead to faster computation, but higher uncertainty in the evaluation part. If \code{NULL}, all observations are used (as by default). Note that performance statistic \code{"auc"} (see argument \code{stats} of \code{\link[=summary.vsel]{summary.vsel()}} and -\code{\link[=plot.vsel]{plot.vsel()}}) is not supported in case of subsampled LOO-CV.} +\code{\link[=plot.vsel]{plot.vsel()}}) is not supported in case of subsampled LOO-CV. Furthermore, +option \code{"best"} for argument \code{baseline} of \code{\link[=summary.vsel]{summary.vsel()}} and +\code{\link[=plot.vsel]{plot.vsel()}} is not supported in case of subsampled LOO-CV.} \item{K}{Only relevant if \code{cv_method = "kfold"} and if \code{cvfits} is \code{NULL} (which is the case for reference model objects created by diff --git a/man/plot.vsel.Rd b/man/plot.vsel.Rd index 69d318d96..89aeb43c8 100644 --- a/man/plot.vsel.Rd +++ b/man/plot.vsel.Rd @@ -99,7 +99,7 @@ estimate.} \item{baseline}{For \code{\link[=summary.vsel]{summary.vsel()}}: Only relevant if \code{deltas} is \code{TRUE}. For \code{\link[=plot.vsel]{plot.vsel()}}: Always relevant. Either \code{"ref"} or \code{"best"}, indicating whether the baseline is the reference model or the best submodel found (in -terms of \code{stats[1]}), respectively.} +terms of \code{stats[1]}), respectively. In case of subsampled LOO-CV, \code{baseline = "best"} is not supported.} \item{thres_elpd}{Only relevant if \code{any(stats \%in\% c("elpd", "mlpd", "gmpd"))}. The threshold for the ELPD difference (taking the submodel's ELPD minus the baseline model's ELPD) above which the submodel's ELPD is diff --git a/man/summary.vsel.Rd b/man/summary.vsel.Rd index 5534f4142..336d91eb6 100644 --- a/man/summary.vsel.Rd +++ b/man/summary.vsel.Rd @@ -99,7 +99,7 @@ estimate.} \item{baseline}{For \code{\link[=summary.vsel]{summary.vsel()}}: Only relevant if \code{deltas} is \code{TRUE}. For \code{\link[=plot.vsel]{plot.vsel()}}: Always relevant. Either \code{"ref"} or \code{"best"}, indicating whether the baseline is the reference model or the best submodel found (in -terms of \code{stats[1]}), respectively.} +terms of \code{stats[1]}), respectively. In case of subsampled LOO-CV, \code{baseline = "best"} is not supported.} \item{resp_oscale}{Only relevant for the latent projection. A single logical value indicating whether to calculate the performance statistics on the From dd0390a7a9625443257caec3a0b8a22ee81bf0a6 Mon Sep 17 00:00:00 2001 From: fweber144 Date: Sat, 1 Feb 2025 07:03:44 +0100 Subject: [PATCH 29/40] vignettes: mention subsampled LOO --- vignettes/projpred.Rmd | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/vignettes/projpred.Rmd b/vignettes/projpred.Rmd index 39c7aded2..b770434c4 100755 --- a/vignettes/projpred.Rmd +++ b/vignettes/projpred.Rmd @@ -123,7 +123,7 @@ The search part determines the predictor ranking (also known as solution path), The evaluation part determines the predictive performance of the increasingly complex submodels along the predictor ranking. There are two functions for running the combination of search and evaluation: `varsel()` and `cv_varsel()`. -In contrast to `varsel()`, `cv_varsel()` performs a cross-validation (CV). With `cv_method = "LOO"` (the default), `cv_varsel()` runs a Pareto-smoothed importance sampling leave-one-out CV [PSIS-LOO CV, see @vehtari_practical_2017; @vehtari_pareto_2022]. With `cv_method = "kfold"`, `cv_varsel()` runs a $K$-fold CV. The extent of the CV depends on `cv_varsel()`'s argument `validate_search`: If `validate_search = TRUE` (the default), the search part is run with the training data of each CV fold separately and the evaluation part is run with the corresponding test data of each CV fold. +In contrast to `varsel()`, `cv_varsel()` performs a cross-validation (CV). With `cv_method = "LOO"` (the default), `cv_varsel()` runs a Pareto-smoothed importance sampling leave-one-out CV [PSIS-LOO CV, see @vehtari_practical_2017; @vehtari_pareto_2022]. With `cv_method = "kfold"`, `cv_varsel()` runs a $K$-fold CV. The extent of the CV mainly depends on `cv_varsel()`'s argument `validate_search`: If `validate_search = TRUE` (the default), the search part is run with the training data of each CV fold separately and the evaluation part is run with the corresponding test data of each CV fold. If `validate_search = FALSE`, the search is excluded from the CV so that only a single full-data search is run. Because of its most thorough protection against overfitting^[Currently, neither `varsel()` nor `cv_varsel()` (not even `cv_varsel()` with `validate_search = TRUE`) guard against overfitting in the selection of the submodel *size*. This is why we added "approximately" to "valid post-selection inference" in section ["Introduction"](#intro). Typically, however, the overfitting induced by the size selection should be comparatively small [@piironen_comparison_2017].], `cv_varsel()` with `validate_search = TRUE` is recommended over `varsel()` and `cv_varsel()` with `validate_search = FALSE`. Nonetheless, a preliminary and comparatively fast run of `varsel()` or `cv_varsel()` with `validate_search = FALSE` can give a rough idea of the performance of the submodels and can be used for finding a suitable value for argument `nterms_max` in subsequent runs (argument `nterms_max` imposes a limit on the submodel size up to which the search is continued and is thus able to reduce the runtime considerably). @@ -190,7 +190,7 @@ Here, we skip this for the sake of brevity and instead head over to the final `c ### Final `cv_varsel()` run For this final `cv_varsel()` run (with `validate_search = TRUE`, as recommended), we use a $K$-fold CV with a small number of folds (`K = 2`) to make this vignette build faster. -In practice, we recommend using either the default of `cv_method = "LOO"` or a larger value for `K` if this is possible in terms of computation time. +In practice, we recommend using either the default of `cv_method = "LOO"` (possibly subsampled, see argument `nloo` of `cv_varsel()`) or a larger value for `K` if this is possible in terms of computation time. Here, we also perform the $K$ reference model refits outside of `cv_varsel()`. Although not strictly necessary here, this is helpful in practice because often, `cv_varsel()` needs to be re-run multiple times in order to try out different argument settings. @@ -554,9 +554,11 @@ Some speed-up possibilities are: In case of `cv_method = "LOO"`^[In case of `cv_method = "kfold"`, the runtime advantage of `validate_search = FALSE` compared to `validate_search = TRUE` is not as large as in case of `cv_method = "LOO"`, but even for `cv_method = "kfold"`, such a runtime advantage still exists.], `cv_varsel()` with `validate_search = FALSE` has comparable runtime to `varsel()`, but accounts for some overfitting, namely that induced by `varsel()`'s in-sample predictions during the predictive performance evaluation. However, as explained in section ["Variable selection"](#variableselection) (see also section ["Overfitting"](#overfitting)), `cv_varsel()` with `validate_search = FALSE` is more prone to overfitting than `cv_varsel()` with `validate_search = TRUE`. +1. Using `cv_varsel()` with subsampled PSIS-LOO CV, see argument `nloo` of `cv_varsel()`. + 1. Using `cv_varsel()` with $K$-fold CV instead of PSIS-LOO CV. - Whether this provides a speed improvement mainly depends on the number of observations and the complexity of the reference model. - Note that PSIS-LOO CV is often more accurate than $K$-fold CV if argument `K` is (much) smaller than the number of observations. + Whether this provides a speed improvement mainly depends on the number of observations, argument `nloo` of `cv_varsel()`, and the complexity of the reference model. + Note that PSIS-LOO CV is often more accurate than $K$-fold CV if argument `K` is (much) smaller than argument `nloo` of `cv_varsel()`. 1. Using a "custom" reference model object with a dimension reduction technique for the predictor data (e.g., by computing principal components from the original predictors, using these principal components as predictors when fitting the reference model, and then performing the variable selection in terms of the *original* predictor terms). Examples are given in @piironen_projective_2020 and @pavone_using_2022. From 79c7b25ecf7356c427ece6cbfa895eef25723094 Mon Sep 17 00:00:00 2001 From: fweber144 Date: Sun, 2 Feb 2025 21:08:26 +0100 Subject: [PATCH 30/40] update `NEWS.md` --- NEWS.md | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 400383f95..38bc390f8 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,9 +4,16 @@ If you read this from a place other than ). This threshold depends on the Monte Carlo sample size and is often close to the former fixed threshold of 0.7 (a short introduction may also be found in the [LOO glossary](https://mc-stan.org/loo/reference/loo-glossary.html)). Correspondingly, the former "secondary" threshold of 0.5 is not used anymore either. +* Use the updated threshold for high Pareto-$\hat{k}$ values presented by Vehtari et al. (2024, "Pareto smoothed importance sampling", *Journal of Machine Learning Research*, 25(72):1-58, ). This threshold depends on the Monte Carlo sample size and is often close to the former fixed threshold of 0.7 (a short introduction may also be found in the [LOO glossary](https://mc-stan.org/loo/reference/loo-glossary.html)). Correspondingly, the former "secondary" threshold of 0.5 is not used anymore either. (GitHub: #490, #498) # projpred 2.8.0 From 6bdade65a0eabed4c34dd6a510cd0ecab6db7be6 Mon Sep 17 00:00:00 2001 From: fweber144 Date: Sun, 2 Feb 2025 21:10:27 +0100 Subject: [PATCH 31/40] fixup! docs: mention log-normal approximation for MSE and RMSE, see --- R/methods.R | 12 ++++++++---- man/plot.vsel.Rd | 12 ++++++++---- man/summary.vsel.Rd | 12 ++++++++---- 3 files changed, 24 insertions(+), 12 deletions(-) diff --git a/R/methods.R b/R/methods.R index 2ea886961..5dc9a2442 100644 --- a/R/methods.R +++ b/R/methods.R @@ -1163,9 +1163,11 @@ plot.vsel <- function( #' * `"elpd"`: expected log (pointwise) predictive density (for a new #' dataset) (ELPD). Estimated by the sum of the observation-specific log #' predictive density values (with each of these predictive density values -#' being a---possibly weighted---average across the parameter draws). +#' being a---possibly weighted---average across the parameter draws). For the +#' corresponding confidence interval, a normal approximation is used. #' * `"mlpd"`: mean log predictive density (MLPD), that is, the ELPD divided -#' by the number of observations. +#' by the number of observations. For the corresponding confidence interval, a +#' normal approximation is used. #' * `"gmpd"`: geometric mean predictive density (GMPD), that is, [exp()] of #' the MLPD. The GMPD is especially helpful for discrete response families #' (because there, the GMPD is bounded by zero and one). For the corresponding @@ -1184,13 +1186,15 @@ plot.vsel <- function( #' approximation is used if `deltas` is `TRUE`. #' * `"R2"`: R-squared, i.e., coefficient of determination (only available in #' the situations mentioned in section "Details" below). For the corresponding -#' standard error, the delta method is used. +#' standard error, the delta method is used. For the corresponding confidence +#' interval, a normal approximation is used. #' * `"acc"` (or its alias, `"pctcorr"`): classification accuracy (only #' available in the situations mentioned in section "Details" below). By #' "classification accuracy", we mean the proportion of correctly classified #' observations. For this, the response category ("class") with highest #' probability (the probabilities are model-based) is taken as the prediction -#' ("classification") for an observation. +#' ("classification") for an observation. For the corresponding confidence +#' interval, a normal approximation is used. #' * `"auc"`: area under the ROC curve (only available in the situations #' mentioned in section "Details" below). For the corresponding standard error #' and lower and upper confidence interval bounds, bootstrapping is used. Not diff --git a/man/plot.vsel.Rd b/man/plot.vsel.Rd index 89aeb43c8..ee30aab2d 100644 --- a/man/plot.vsel.Rd +++ b/man/plot.vsel.Rd @@ -48,9 +48,11 @@ set). Available statistics are: \item \code{"elpd"}: expected log (pointwise) predictive density (for a new dataset) (ELPD). Estimated by the sum of the observation-specific log predictive density values (with each of these predictive density values -being a---possibly weighted---average across the parameter draws). +being a---possibly weighted---average across the parameter draws). For the +corresponding confidence interval, a normal approximation is used. \item \code{"mlpd"}: mean log predictive density (MLPD), that is, the ELPD divided -by the number of observations. +by the number of observations. For the corresponding confidence interval, a +normal approximation is used. \item \code{"gmpd"}: geometric mean predictive density (GMPD), that is, \code{\link[=exp]{exp()}} of the MLPD. The GMPD is especially helpful for discrete response families (because there, the GMPD is bounded by zero and one). For the corresponding @@ -69,13 +71,15 @@ a log-normal approximation is used if \code{deltas} is \code{FALSE} and a normal approximation is used if \code{deltas} is \code{TRUE}. \item \code{"R2"}: R-squared, i.e., coefficient of determination (only available in the situations mentioned in section "Details" below). For the corresponding -standard error, the delta method is used. +standard error, the delta method is used. For the corresponding confidence +interval, a normal approximation is used. \item \code{"acc"} (or its alias, \code{"pctcorr"}): classification accuracy (only available in the situations mentioned in section "Details" below). By "classification accuracy", we mean the proportion of correctly classified observations. For this, the response category ("class") with highest probability (the probabilities are model-based) is taken as the prediction -("classification") for an observation. +("classification") for an observation. For the corresponding confidence +interval, a normal approximation is used. \item \code{"auc"}: area under the ROC curve (only available in the situations mentioned in section "Details" below). For the corresponding standard error and lower and upper confidence interval bounds, bootstrapping is used. Not diff --git a/man/summary.vsel.Rd b/man/summary.vsel.Rd index 336d91eb6..cc131e62d 100644 --- a/man/summary.vsel.Rd +++ b/man/summary.vsel.Rd @@ -38,9 +38,11 @@ set). Available statistics are: \item \code{"elpd"}: expected log (pointwise) predictive density (for a new dataset) (ELPD). Estimated by the sum of the observation-specific log predictive density values (with each of these predictive density values -being a---possibly weighted---average across the parameter draws). +being a---possibly weighted---average across the parameter draws). For the +corresponding confidence interval, a normal approximation is used. \item \code{"mlpd"}: mean log predictive density (MLPD), that is, the ELPD divided -by the number of observations. +by the number of observations. For the corresponding confidence interval, a +normal approximation is used. \item \code{"gmpd"}: geometric mean predictive density (GMPD), that is, \code{\link[=exp]{exp()}} of the MLPD. The GMPD is especially helpful for discrete response families (because there, the GMPD is bounded by zero and one). For the corresponding @@ -59,13 +61,15 @@ a log-normal approximation is used if \code{deltas} is \code{FALSE} and a normal approximation is used if \code{deltas} is \code{TRUE}. \item \code{"R2"}: R-squared, i.e., coefficient of determination (only available in the situations mentioned in section "Details" below). For the corresponding -standard error, the delta method is used. +standard error, the delta method is used. For the corresponding confidence +interval, a normal approximation is used. \item \code{"acc"} (or its alias, \code{"pctcorr"}): classification accuracy (only available in the situations mentioned in section "Details" below). By "classification accuracy", we mean the proportion of correctly classified observations. For this, the response category ("class") with highest probability (the probabilities are model-based) is taken as the prediction -("classification") for an observation. +("classification") for an observation. For the corresponding confidence +interval, a normal approximation is used. \item \code{"auc"}: area under the ROC curve (only available in the situations mentioned in section "Details" below). For the corresponding standard error and lower and upper confidence interval bounds, bootstrapping is used. Not From d06f43d4a20525860661fc0039ab08a5d0de7bfc Mon Sep 17 00:00:00 2001 From: fweber144 Date: Sun, 2 Feb 2025 21:23:24 +0100 Subject: [PATCH 32/40] fixup! subsampled LOO and AUC: switch from warning to error, see --- R/summary_funs.R | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/R/summary_funs.R b/R/summary_funs.R index 9b7f8dfd6..8903e373c 100644 --- a/R/summary_funs.R +++ b/R/summary_funs.R @@ -528,8 +528,9 @@ get_stat <- function(summaries, summaries_baseline = NULL, } } else if (stat == "auc") { if (n_loo < n_full) { - # subsampling LOO with AUC not implemented. Using fast LOO result. - mu <- mu_fast + # Note: Previously, subsampled LOO with AUC caused the fast LOO results + # to be used automatically (via `mu <- mu_fast`), see PR #496. + stop("Subsampled LOO-CV with AUC not implemented.") } if (!is.null(mu_baseline)) { auc_data <- cbind(y, mu, wobs) From a33e614502a1ffe6ca51c2cf2ad00648d81a476b Mon Sep 17 00:00:00 2001 From: fweber144 Date: Tue, 4 Feb 2025 20:58:17 +0100 Subject: [PATCH 33/40] set `nloo` to `NULL` in case of K-fold CV --- R/cv_varsel.R | 3 ++- man/cv_varsel.Rd | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/R/cv_varsel.R b/R/cv_varsel.R index 05e4e203f..6f6f6d117 100644 --- a/R/cv_varsel.R +++ b/R/cv_varsel.R @@ -272,7 +272,7 @@ cv_varsel.refmodel <- function( nterms_max = NULL, penalty = NULL, verbose = getOption("projpred.verbose", interactive()), - nloo = object$nobs, + nloo = if (cv_method == "LOO") object$nobs else NULL, K = if (!inherits(object, "datafit")) 5 else 10, cvfits = object$cvfits, search_control = NULL, @@ -560,6 +560,7 @@ parse_args_cv_varsel <- function(refmodel, cv_method, nloo, K, cvfits, stop("For K-fold-CV, `validate_search = FALSE` may not be combined with ", "`refit_prj = FALSE`.") } + nloo <- NULL } else { stopifnot(!is.null(refmodel[["nobs"]])) nloo <- min(nloo, refmodel[["nobs"]]) diff --git a/man/cv_varsel.Rd b/man/cv_varsel.Rd index dae0d10cc..eafb9ab63 100644 --- a/man/cv_varsel.Rd +++ b/man/cv_varsel.Rd @@ -33,7 +33,7 @@ cv_varsel(object, ...) nterms_max = NULL, penalty = NULL, verbose = getOption("projpred.verbose", interactive()), - nloo = object$nobs, + nloo = if (cv_method == "LOO") object$nobs else NULL, K = if (!inherits(object, "datafit")) 5 else 10, cvfits = object$cvfits, search_control = NULL, From 944899feda30f0ef490d440a6d1eed005cb7d453 Mon Sep 17 00:00:00 2001 From: fweber144 Date: Mon, 3 Feb 2025 20:11:57 +0100 Subject: [PATCH 34/40] modify the tests so that subsampled LOO is tested more thoroughly --- tests/testthat/helpers/testers.R | 11 +- tests/testthat/setup.R | 35 +++-- tests/testthat/test_datafit.R | 5 + tests/testthat/test_varsel.R | 211 ++++++------------------------- 4 files changed, 80 insertions(+), 182 deletions(-) diff --git a/tests/testthat/helpers/testers.R b/tests/testthat/helpers/testers.R index 00b9fdb96..4d2422fe9 100644 --- a/tests/testthat/helpers/testers.R +++ b/tests/testthat/helpers/testers.R @@ -1968,7 +1968,11 @@ vsel_tester <- function( nclusters_pred_tst }, seed_expected = seed_tst, - nloo_expected = if (with_cv) refmod_expected$nobs else NULL, + nloo_expected = if (with_cv && !identical(cv_method_expected, "kfold")) { + refmod_expected$nobs + } else { + NULL + }, K_expected = NULL, penalty_expected = NULL, search_terms_expected = NULL, @@ -2001,7 +2005,6 @@ vsel_tester <- function( # size (see issue #307): prd_trms_len_expected <- prd_trms_len_expected - 1L } - nloo_expected_orig <- nloo_expected # Test the general structure of the object: expect_s3_class(vs, "vsel") @@ -2250,7 +2253,7 @@ vsel_tester <- function( expect_named(vs$summaries, c("sub", "ref"), info = info_str) expect_type(vs$summaries$sub, "list") expect_length(vs$summaries$sub, prd_trms_len_expected + 1) - if (with_cv) { + if (with_cv && identical(cv_method_expected, "LOO")) { if (is.null(nloo_expected) || nloo_expected > nobsv) { nloo_expected <- nobsv } @@ -2374,7 +2377,7 @@ vsel_tester <- function( expect_identical(vs$cv_method, cv_method_expected, info = info_str) # nloo - expect_identical(vs$nloo, nloo_expected_orig, info = info_str) + expect_identical(vs$nloo, nloo_expected, info = info_str) # K if (!is.null(K_expected)) { diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R index 623a33ba5..5f562e3aa 100755 --- a/tests/testthat/setup.R +++ b/tests/testthat/setup.R @@ -937,6 +937,11 @@ cvmeth_tst <- list( kfold = list(cv_method = "kfold") ) +nloo_tst <- list( + default_nloo = list(), + subsmpl = list(nloo = as.integer(nobsv %/% 10)) +) + resp_oscale_tst <- list( default_r_oscale = list(), r_oscale_F = list(resp_oscale = FALSE) @@ -1233,6 +1238,14 @@ if (run_cvvs) { !run_valsearch_aug_lat))) { cvmeth_i <- c(cvmeth_i, list(validate_search = FALSE)) } + if (identical(cvmeth_i$cv_method, "kfold")) { + nloo_tst <- nloo_tst["default_nloo"] + } else if (!((prj_crr == "trad" && mod_crr == "glm" && + fam_crr == "gauss") || + (prj_crr %in% c("augdat", "latent") && mod_crr == "glm" && + fam_crr == "cumul"))) { + nloo_tst <- nloo_tst["subsmpl"] + } if (run_more && mod_crr == "glm" && fam_crr == "gauss" && grepl("\\.stdformul\\.", tstsetup_ref)) { # Here, we also test non-NULL `search_terms`: @@ -1247,14 +1260,20 @@ if (run_cvvs) { nterms_max_tst <- count_terms_chosen(search_trms_i$search_terms) - 1L } - return(c( - nlist(tstsetup_ref), only_nonargs(args_ref[[tstsetup_ref]]), - list( - nclusters = nclusters_tst, nclusters_pred = nclusters_pred_tst, - nterms_max = nterms_max_tst, verbose = FALSE, seed = seed_tst - ), - meth_i, cvmeth_i, search_trms_i - )) + lapply(nloo_tst, function(nloo_i) { + if (!is.null(nloo_i$nloo) && nloo_i$nloo < nobsv && + identical(cvmeth_i$validate_search, FALSE)) { + cvmeth_i$validate_search <- NULL + } + return(c( + nlist(tstsetup_ref), only_nonargs(args_ref[[tstsetup_ref]]), + list( + nclusters = nclusters_tst, nclusters_pred = nclusters_pred_tst, + nterms_max = nterms_max_tst, verbose = FALSE, seed = seed_tst + ), + meth_i, cvmeth_i, nloo_i, search_trms_i + )) + }) }) }) }) diff --git a/tests/testthat/test_datafit.R b/tests/testthat/test_datafit.R index 0859ae027..68f119851 100644 --- a/tests/testthat/test_datafit.R +++ b/tests/testthat/test_datafit.R @@ -114,12 +114,17 @@ if (run_cvvs) { args_cvvs_datafit <- lapply(args_cvvs_datafit, function(args_cvvs_i) { args_cvvs_i$cv_method <- NULL args_cvvs_i$K <- NULL + args_cvvs_i$nloo <- NULL args_cvvs_i$validate_search <- TRUE return(c(args_cvvs_i, list(cv_method = "kfold"))) }) + names(args_cvvs_datafit) <- gsub("(\\.default_cvmeth\\..*)\\.default_nloo", + "\\1.subsmpl", names(args_cvvs_datafit)) names(args_cvvs_datafit) <- gsub("default_cvmeth", "kfold", names(args_cvvs_datafit)) args_cvvs_datafit <- args_cvvs_datafit[unique(names(args_cvvs_datafit))] + names(args_cvvs_datafit) <- gsub("\\.subsmpl", "\\.default_nloo", + names(args_cvvs_datafit)) # For `datafit`s, we always have 1 cluster by default, so omit related # arguments: args_cvvs_datafit <- lapply(args_cvvs_datafit, function(args_cvvs_i) { diff --git a/tests/testthat/test_varsel.R b/tests/testthat/test_varsel.R index 6db89b676..dc35ba87b 100644 --- a/tests/testthat/test_varsel.R +++ b/tests/testthat/test_varsel.R @@ -1195,6 +1195,7 @@ test_that(paste( prd_trms_len_expected = args_cvvs[[tstsetup]]$nterms_max, method_expected = meth_exp_crr, cv_method_expected = args_cvvs[[tstsetup]]$cv_method, + nloo_expected = args_cvvs[[tstsetup]]$nloo, valsearch_expected = args_cvvs[[tstsetup]]$validate_search, search_terms_expected = args_cvvs[[tstsetup]]$search_terms, search_trms_empty_size = @@ -1321,6 +1322,7 @@ test_that("`refit_prj` works", { method_expected = meth_exp_crr, refit_prj_expected = FALSE, cv_method_expected = args_cvvs_i$cv_method, + nloo_expected = args_cvvs_i$nloo, valsearch_expected = args_cvvs_i$validate_search, search_terms_expected = args_cvvs_i$search_terms, search_trms_empty_size = @@ -1396,7 +1398,7 @@ test_that("invalid `nloo` fails", { suppressWarnings(do.call(cv_varsel, c( list(object = refmods[[args_cvvs_i$tstsetup_ref]], nloo = -1), - excl_nonargs(args_cvvs_i) + excl_nonargs(args_cvvs_i, nms_excl_add = "nloo") ))), "^nloo must be at least 1$", info = tstsetup @@ -1411,7 +1413,7 @@ test_that(paste( skip_if_not(run_cvvs) nloo_tst <- nobsv + 1L tstsetups <- grep( - "\\.glm\\.gauss\\..*\\.default_cvmeth\\.default_search_trms", + "\\.glm\\.gauss\\..*\\.default_cvmeth\\.default_search_trms\\.default_nloo", names(cvvss), value = TRUE ) valsearch_arg <- lapply(args_cvvs[tstsetups], "[[", "validate_search") @@ -1426,26 +1428,23 @@ test_that(paste( nloo = nloo_tst), excl_nonargs(args_cvvs_i) ))) - cvvs_nloo[["nloo"]] <- nobsv expect_equal(cvvs_nloo, cvvss[[tstsetup]], info = tstsetup) } }) test_that("setting `nloo` smaller than the number of observations works", { skip_if_not(run_cvvs) - nloo_tst <- nobsv %/% 5L # Output elements of `vsel` objects that may be influenced by `nloo`: - vsel_nms_nloo <- c("summaries", "summaries_fast","predictor_ranking_cv", + vsel_nms_nloo <- c("summaries", "summaries_fast", "predictor_ranking_cv", "nloo", "loo_inds", "ce") # In general, element `ce` is affected as well (because the PRNG state when # doing the clustering for the performance evaluation is different when `nloo` # is smaller than the number of observations compared to when `nloo` is equal - # to the number of observations), but the changes in `ce` may be so small that - # they are not detected by all.equal(): + # to the number of observations): vsel_nms_nloo_opt <- c("ce") # The setups that should be tested: tstsetups <- grep( - "\\.glm\\.gauss\\..*\\.default_cvmeth\\.default_search_trms", + "\\.default_cvmeth\\.default_search_trms\\.default_nloo", names(cvvss), value = TRUE ) valsearch_arg <- lapply(args_cvvs[tstsetups], "[[", "validate_search") @@ -1453,34 +1452,9 @@ test_that("setting `nloo` smaller than the number of observations works", { skip_if(length(tstsetups) == 0) for (tstsetup in tstsetups) { args_cvvs_i <- args_cvvs[[tstsetup]] - tstsetup_ref <- args_cvvs_i$tstsetup_ref - mod_crr <- args_cvvs_i$mod_nm - fam_crr <- args_cvvs_i$fam_nm - prj_crr <- args_cvvs_i$prj_nm - meth_exp_crr <- args_cvvs_i$method %||% "forward" - # Use suppressWarnings() because test_that() somehow redirects stderr() and - # so throws warnings that projpred wants to capture internally: - cvvs_nloo <- suppressWarnings(do.call(cv_varsel, c( - list(object = refmods[[args_cvvs_i$tstsetup_ref]], - nloo = nloo_tst), - excl_nonargs(args_cvvs_i) - ))) - vsel_tester( - cvvs_nloo, - with_cv = TRUE, - refmod_expected = refmods[[tstsetup_ref]], - prd_trms_len_expected = args_cvvs_i$nterms_max, - method_expected = meth_exp_crr, - cv_method_expected = "LOO", - valsearch_expected = TRUE, - nloo_expected = nloo_tst, - search_terms_expected = args_cvvs_i$search_terms, - search_trms_empty_size = - length(args_cvvs_i$search_terms) && - all(grepl("\\+", args_cvvs_i$search_terms)), - search_control_expected = args_cvvs_i[c("avoid.increase")], - info_str = tstsetup - ) + tstsetup_nloo <- sub("\\.default_nloo", "\\.subsmpl", tstsetup) + stopifnot(tstsetup_nloo %in% names(cvvss)) + cvvs_nloo <- cvvss[[tstsetup_nloo]] # Expected equality for most elements with a few exceptions: vsel_nms_nloo_crr <- vsel_nms_nloo if (isFALSE(args_cvvs_i$validate_search)) { @@ -1511,7 +1485,8 @@ test_that("`validate_search` works", { tstsetups <- names(cvvss) if (!run_valsearch_always) { has_valsearch_true <- sapply(tstsetups, function(tstsetup_cvvs) { - !isFALSE(args_cvvs[[tstsetup_cvvs]]$validate_search) + !isFALSE(args_cvvs[[tstsetup_cvvs]]$validate_search) && + (args_cvvs[[tstsetup_cvvs]]$nloo %||% nobsv) == nobsv }) tstsetups <- tstsetups[has_valsearch_true] } @@ -1610,7 +1585,9 @@ test_that("`validate_search` works", { suggsize_cond[tstsetup] <- sgg_size_valsearch >= sgg_size } } - expect_true(mean(!suggsize_cond, na.rm = TRUE) <= 0.25) + if (!all(is.na(suggsize_cond))) { + expect_true(mean(!suggsize_cond, na.rm = TRUE) <= 0.25) + } }) ## Arguments specific to K-fold CV ---------------------------------------- @@ -1913,7 +1890,7 @@ test_that(paste( skip_if_not(run_cvvs) tstsetups <- names(cvvss) if (!run_more) { - tstsetups <- head(tstsetups, 1) + tstsetups <- head(tstsetups, 2) } for (tstsetup in tstsetups) { refit_prj_crr <- !identical(args_cvvs[[tstsetup]]$validate_search, FALSE) || @@ -1922,8 +1899,9 @@ test_that(paste( # Use suppressWarnings() because test_that() somehow redirects stderr() and # so throws warnings that projpred wants to capture internally: cvvs_eval <- suppressWarnings(cv_varsel( - cvvss[[tstsetup]], refit_prj = refit_prj_crr, - nclusters_pred = nclusters_pred_crr, verbose = FALSE, seed = seed2_tst + cvvss[[tstsetup]], + refit_prj = refit_prj_crr, nclusters_pred = nclusters_pred_crr, + verbose = FALSE, seed = seed2_tst )) tstsetup_ref <- args_cvvs[[tstsetup]]$tstsetup_ref meth_exp_crr <- args_cvvs[[tstsetup]]$method %||% "forward" @@ -1940,6 +1918,7 @@ test_that(paste( prd_trms_len_expected = args_cvvs[[tstsetup]]$nterms_max, method_expected = meth_exp_crr, cv_method_expected = args_cvvs[[tstsetup]]$cv_method, + nloo_expected = args_cvvs[[tstsetup]]$nloo, valsearch_expected = args_cvvs[[tstsetup]]$validate_search, refit_prj_expected = refit_prj_crr, nprjdraws_eval_expected = if (!refit_prj_crr && meth_exp_crr == "L1") { @@ -2093,8 +2072,8 @@ test_that(paste( )) } else { cvvs_eval <- suppressWarnings(cv_varsel( - cvvss[[tstsetup]], validate_search = FALSE, refit_prj = FALSE, - verbose = FALSE, seed = seed2_tst + cvvss[[tstsetup]], nloo = nobsv, validate_search = FALSE, + refit_prj = FALSE, verbose = FALSE, seed = seed2_tst )) } tstsetup_ref <- args_cvvs[[tstsetup]]$tstsetup_ref @@ -2146,7 +2125,7 @@ test_that(paste( skip_if_not(run_cvvs) tstsetups <- names(cvvss) if (!run_more) { - tstsetups <- head(tstsetups, 8) + tstsetups <- head(tstsetups, 11) } for (tstsetup in tstsetups) { tstsetup_ref <- args_cvvs[[tstsetup]]$tstsetup_ref @@ -2206,7 +2185,7 @@ test_that(paste( skip_if_not(run_cvvs) tstsetups <- names(cvvss) if (!run_more) { - tstsetups <- head(tstsetups, 2) + tstsetups <- head(tstsetups, 3) } for (tstsetup in tstsetups) { tstsetup_ref <- args_cvvs[[tstsetup]]$tstsetup_ref @@ -2269,7 +2248,7 @@ test_that(paste( skip_if_not(run_cvvs) tstsetups <- names(cvvss) if (!run_more) { - tstsetups <- head(tstsetups, 8) + tstsetups <- head(tstsetups, 11) } for (tstsetup in tstsetups) { tstsetup_ref <- args_cvvs[[tstsetup]]$tstsetup_ref @@ -2291,12 +2270,14 @@ test_that(paste( # so throws warnings that projpred wants to capture internally: if (identical(args_cvvs[[tstsetup]]$cv_method, "kfold")) { cv_meth_crr <- "LOO" + nloo_crr <- nloo_tst[["subsmpl"]][["nloo"]] cvvs_eval <- suppressWarnings(cv_varsel( - cvvss[[tstsetup]], cv_method = cv_meth_crr, + cvvss[[tstsetup]], cv_method = cv_meth_crr, nloo = nloo_crr, nclusters_pred = nclusters_pred_crr, verbose = FALSE, seed = seed2_tst )) } else { cv_meth_crr <- "kfold" + nloo_crr <- NULL cvvs_eval <- suppressWarnings(cv_varsel( cvvss[[tstsetup]], cv_method = cv_meth_crr, K = K_tst, cvfits = cvfitss[[tstsetup_ref]], nclusters_pred = nclusters_pred_crr, @@ -2312,6 +2293,7 @@ test_that(paste( prd_trms_len_expected = args_cvvs[[tstsetup]]$nterms_max, method_expected = meth_exp_crr, cv_method_expected = cv_meth_crr, + nloo_expected = nloo_crr, valsearch_expected = TRUE, K_expected = K_tst, nprjdraws_eval_expected = nclusters_pred_crr, @@ -2341,12 +2323,13 @@ test_that(paste( } else { nclusters_pred_crr <- args_vs[[tstsetup]]$nclusters_pred } + nloo_crr <- nloo_tst[["subsmpl"]][["nloo"]] # Use suppressWarnings() because test_that() somehow redirects stderr() and # so throws warnings that projpred wants to capture internally: cvvs_eval <- try( suppressWarnings(cv_varsel( - vss[[tstsetup]], nclusters_pred = nclusters_pred_crr, verbose = FALSE, - seed = seed2_tst + vss[[tstsetup]], nloo = nloo_crr, nclusters_pred = nclusters_pred_crr, + verbose = FALSE, seed = seed2_tst )), silent = TRUE ) @@ -2375,6 +2358,7 @@ test_that(paste( prd_trms_len_expected = args_vs[[tstsetup]]$nterms_max, method_expected = meth_exp_crr, cv_method_expected = "LOO", + nloo_expected = nloo_crr, valsearch_expected = TRUE, nprjdraws_eval_expected = nclusters_pred_crr, search_terms_expected = args_vs[[tstsetup]]$search_terms, @@ -2452,7 +2436,6 @@ test_that(paste( prd_trms_len_expected = args_vs[[tstsetup]]$nterms_max, method_expected = meth_exp_crr, cv_method_expected = "kfold", - nloo_expected = NULL, valsearch_expected = TRUE, nprjdraws_eval_expected = nclusters_pred_crr, search_terms_expected = args_vs[[tstsetup]]$search_terms, @@ -2473,7 +2456,7 @@ test_that(paste( skip_if_not(run_cvvs) tstsetups <- names(cvvss) if (!run_more) { - tstsetups <- head(tstsetups, 2) + tstsetups <- head(tstsetups, 3) } for (tstsetup in tstsetups) { if (run_more && !args_cvvs[[tstsetup]]$mod_nm %in% c("glm", "gam")) { @@ -2487,14 +2470,17 @@ test_that(paste( # Use suppressWarnings() because test_that() somehow redirects stderr() and # so throws warnings that projpred wants to capture internally: if (identical(args_cvvs[[tstsetup]]$cv_method, "kfold")) { + nloo_crr <- nloo_tst[["subsmpl"]][["nloo"]] cvvs_eval <- try( suppressWarnings(cv_varsel( - cvvss[[tstsetup]], cv_method = "LOO", validate_search = TRUE, - nclusters_pred = nclusters_pred_crr, verbose = FALSE, seed = seed2_tst + cvvss[[tstsetup]], cv_method = "LOO", nloo = nloo_crr, + validate_search = TRUE, nclusters_pred = nclusters_pred_crr, + verbose = FALSE, seed = seed2_tst )), silent = TRUE ) } else { + nloo_crr <- cvvss[[tstsetup]][["nloo"]] cvvs_eval <- try( suppressWarnings(cv_varsel( cvvss[[tstsetup]], validate_search = TRUE, @@ -2543,6 +2529,7 @@ test_that(paste( prd_trms_len_expected = args_cvvs[[tstsetup]]$nterms_max, method_expected = meth_exp_crr, cv_method_expected = "LOO", + nloo_expected = nloo_crr, valsearch_expected = TRUE, nprjdraws_eval_expected = nclusters_pred_crr, K_expected = if (identical(args_cvvs[[tstsetup]]$cv_method, "kfold")) { @@ -2568,7 +2555,7 @@ test_that(paste( skip_if_not(run_cvvs) tstsetups <- names(cvvss) if (!run_more) { - tstsetups <- head(tstsetups, 2) + tstsetups <- head(tstsetups, 3) } for (tstsetup in tstsetups) { tstsetup_ref <- args_cvvs[[tstsetup]]$tstsetup_ref @@ -2643,122 +2630,6 @@ test_that(paste( } }) -test_that("cv_varsel.vsel(): `nloo` works for `vsel` objects from varsel()", { - skip_if_not(run_vs) - skip_if_not(run_cvvs) - nloo_tst <- nobsv %/% 5L - tstsetup_counter <- 0L - for (tstsetup in names(vss)) { - if (!run_more && tstsetup_counter > 0L) { - next - } else if (run_more && tstsetup_counter > length(vss) %/% 6) { - next - } else if (run_more) { - refit_prj_crr <- TRUE - nclusters_pred_crr <- args_vs[[tstsetup]]$nclusters_pred - 1L - } else { - refit_prj_crr <- FALSE - nclusters_pred_crr <- args_vs[[tstsetup]]$nclusters_pred - } - # Use suppressWarnings() because test_that() somehow redirects stderr() and - # so throws warnings that projpred wants to capture internally: - cvvs_eval_valT <- suppressWarnings(cv_varsel( - vss[[tstsetup]], nloo = nloo_tst, nclusters_pred = nclusters_pred_crr, - verbose = FALSE, seed = seed2_tst - )) - tstsetup_ref <- args_vs[[tstsetup]]$tstsetup_ref - meth_exp_crr <- args_vs[[tstsetup]]$method %||% "forward" - vsel_tester( - cvvs_eval_valT, - with_cv = TRUE, - refmod_expected = refmods[[tstsetup_ref]], - prd_trms_len_expected = args_vs[[tstsetup]]$nterms_max, - method_expected = meth_exp_crr, - cv_method_expected = "LOO", - nloo_expected = nloo_tst, - valsearch_expected = TRUE, - nprjdraws_eval_expected = nclusters_pred_crr, - search_terms_expected = args_vs[[tstsetup]]$search_terms, - search_trms_empty_size = - length(args_vs[[tstsetup]]$search_terms) && - all(grepl("\\+", args_vs[[tstsetup]]$search_terms)), - search_control_expected = args_vs[[tstsetup]][c("avoid.increase")], - info_str = tstsetup - ) - tstsetup_counter <- tstsetup_counter + 1L - } -}) - -test_that(paste( - "cv_varsel.vsel(): `nloo` works for `vsel` objects from cv_varsel()" -), { - skip_if_not(run_cvvs) - nloo_tst <- nobsv %/% 5L - tstsetups <- names(cvvss) - if (!run_more) { - tstsetups <- head(tstsetups, 1) - } else { - tstsetups <- head(grep("\\.glm\\.|\\.gam\\.", tstsetups, value = TRUE), - length(tstsetups) %/% 6) - # Make sure that in the test setups, we have `validate_search = TRUE` as - # well as `validate_search = FALSE`: - valsearches <- !unlist(lapply( - lapply(args_cvvs[tstsetups], "[[", "validate_search"), - isFALSE - )) - stopifnot(any(valsearches), any(!valsearches)) - } - for (tstsetup in tstsetups) { - # Use suppressWarnings() because test_that() somehow redirects stderr() and - # so throws warnings that projpred wants to capture internally: - if (identical(args_cvvs[[tstsetup]]$cv_method, "kfold")) { - cvvs_eval_valT <- suppressWarnings(cv_varsel( - cvvss[[tstsetup]], cv_method = "LOO", validate_search = TRUE, - nloo = nloo_tst, refit_prj = FALSE, verbose = FALSE, seed = seed2_tst - )) - } else { - cvvs_eval_valT <- suppressWarnings(cv_varsel( - cvvss[[tstsetup]], nloo = nloo_tst, validate_search = TRUE, - refit_prj = FALSE, verbose = FALSE, seed = seed2_tst - )) - } - meth_exp_crr <- args_cvvs[[tstsetup]]$method %||% "forward" - vsel_tester( - cvvs_eval_valT, - with_cv = TRUE, - refmod_expected = refmods[[args_cvvs[[tstsetup]]$tstsetup_ref]], - cvfits_expected = if (identical(args_cvvs[[tstsetup]]$cv_method, - "kfold")) { - cvfitss[[args_cvvs[[tstsetup]]$tstsetup_ref]] - } else { - refmods[[args_cvvs[[tstsetup]]$tstsetup_ref]]$cvfits - }, - prd_trms_len_expected = args_cvvs[[tstsetup]]$nterms_max, - method_expected = meth_exp_crr, - cv_method_expected = "LOO", - nloo_expected = nloo_tst, - valsearch_expected = TRUE, - refit_prj_expected = FALSE, - nprjdraws_eval_expected = if (meth_exp_crr == "L1") { - 1L - } else { - nclusters_tst - }, - K_expected = if (identical(args_cvvs[[tstsetup]]$cv_method, "kfold")) { - K_tst - } else { - NULL - }, - search_terms_expected = args_cvvs[[tstsetup]]$search_terms, - search_trms_empty_size = - length(args_cvvs[[tstsetup]]$search_terms) && - all(grepl("\\+", args_cvvs[[tstsetup]]$search_terms)), - search_control_expected = args_cvvs[[tstsetup]][c("avoid.increase")], - info_str = tstsetup - ) - } -}) - # run_cvfun() ------------------------------------------------------------- test_that("argument `folds` of run_cvfun() works", { From 27d47ea7006bf8522646c1662c74691c082a6d45 Mon Sep 17 00:00:00 2001 From: fweber144 Date: Wed, 5 Feb 2025 20:34:38 +0100 Subject: [PATCH 35/40] fix a test (only discovered this now after changing the tests and running with `run_more = TRUE`) --- tests/testthat/test_methods_vsel.R | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/tests/testthat/test_methods_vsel.R b/tests/testthat/test_methods_vsel.R index d5c0523dc..da92f85c2 100644 --- a/tests/testthat/test_methods_vsel.R +++ b/tests/testthat/test_methods_vsel.R @@ -608,6 +608,12 @@ test_that(paste( nterms_max_expected_crr <- args_rk_cvvs[[tstsetup_rk]][["nterms_max"]] if (is.null(nterms_max_expected_crr)) { nterms_max_expected_crr <- args_cvvs[[tstsetup_cvvs]][["nterms_max"]] + if (length(args_cvvs[[tstsetup_cvvs]]$search_terms) && + all(grepl("\\+", args_cvvs[[tstsetup_cvvs]]$search_terms))) { + # This is the "empty_size" setting, so we have to subtract the skipped + # model size (see issue #307): + nterms_max_expected_crr <- nterms_max_expected_crr - 1L + } } cv_proportions_tester( prs_cvvs[[tstsetup]], From 92ad48139c484e39221ab11d51614f3fa197e682 Mon Sep 17 00:00:00 2001 From: fweber144 Date: Fri, 7 Feb 2025 21:09:32 +0100 Subject: [PATCH 36/40] revert everything not related to `progressr` --- R/cv_varsel.R | 208 ++++++++++++++----------------- R/glmfun.R | 2 +- R/misc.R | 8 +- R/summary_funs.R | 20 +-- R/varsel.R | 36 ++---- man/cv_varsel.Rd | 2 +- man/varsel.Rd | 2 +- tests/testthat/helpers/testers.R | 4 +- tests/testthat/test_misc.R | 98 +++++++-------- 9 files changed, 173 insertions(+), 207 deletions(-) diff --git a/R/cv_varsel.R b/R/cv_varsel.R index 6f6f6d117..9f09600d0 100644 --- a/R/cv_varsel.R +++ b/R/cv_varsel.R @@ -191,20 +191,16 @@ cv_varsel.vsel <- function( validate_search = object$validate_search %||% TRUE, ... ) { - ## the following arguments should not change arg_nms_internal <- c("method", "ndraws", "nclusters", "nterms_max", "search_control", "penalty", "search_terms") - dots <- list(...) - arg_nms_internal_used <- intersect(arg_nms_internal, names(dots)) - for (arg in arg_nms_internal_used) { - if (!identical(object[["args_search"]][[arg]], dots[[arg]])) { - message("Argument `", arg, "` ignored. Using the argument value stored ", - "in the `vsel` object.") - } - ## remove duplicate arguments - dots[[arg]] <- NULL + arg_nms_internal_used <- intersect(arg_nms_internal, ...names()) + n_arg_nms_internal_used <- length(arg_nms_internal_used) + if (n_arg_nms_internal_used > 0) { + stop("Argument", if (n_arg_nms_internal_used > 1) "s" else "", " ", + paste(paste0("`", arg_nms_internal_used, "`"), collapse = ", "), " ", + "cannot be specified in this case because cv_varsel.vsel() specifies ", + if (n_arg_nms_internal_used > 1) "them" else "it", " ", "internally.") } - refmodel <- get_refmodel(object) rk_foldwise <- ranking(object)[["foldwise"]] if (validate_search && !is.null(rk_foldwise)) { @@ -236,26 +232,23 @@ cv_varsel.vsel <- function( "brms:::get_refmodel.brmsfit() to some non-`NULL` value.") } } - - return(do_call(cv_varsel, c( - list( - object = refmodel, - method = object[["args_search"]][["method"]], - ndraws = object[["args_search"]][["ndraws"]], - nclusters = object[["args_search"]][["nclusters"]], - nterms_max = object[["args_search"]][["nterms_max"]], - search_control = object[["args_search"]][["search_control"]], - penalty = object[["args_search"]][["penalty"]], - search_terms = object[["args_search"]][["search_terms"]], - cv_method = cv_method, - nloo = nloo, - K = K, - cvfits = cvfits, - validate_search = validate_search, - search_out = nlist(search_path = object[["search_path"]], rk_foldwise) - ), - dots - ))) + return(cv_varsel( + object = refmodel, + method = object[["args_search"]][["method"]], + ndraws = object[["args_search"]][["ndraws"]], + nclusters = object[["args_search"]][["nclusters"]], + nterms_max = object[["args_search"]][["nterms_max"]], + search_control = object[["args_search"]][["search_control"]], + penalty = object[["args_search"]][["penalty"]], + search_terms = object[["args_search"]][["search_terms"]], + cv_method = cv_method, + nloo = nloo, + K = K, + cvfits = cvfits, + validate_search = validate_search, + search_out = nlist(search_path = object[["search_path"]], rk_foldwise), + ... + )) } #' @rdname cv_varsel @@ -271,7 +264,7 @@ cv_varsel.refmodel <- function( refit_prj = !inherits(object, "datafit"), nterms_max = NULL, penalty = NULL, - verbose = getOption("projpred.verbose", interactive()), + verbose = TRUE, nloo = if (cv_method == "LOO") object$nobs else NULL, K = if (!inherits(object, "datafit")) 5 else 10, cvfits = object$cvfits, @@ -284,7 +277,8 @@ cv_varsel.refmodel <- function( search_terms = NULL, search_out = NULL, parallel = getOption("projpred.prll_cv", FALSE), - ...) { + ... +) { if (!missing(lambda_min_ratio)) { warning("Argument `lambda_min_ratio` is deprecated. Please specify ", "control arguments for the search via argument `search_control`. ", @@ -346,25 +340,16 @@ cv_varsel.refmodel <- function( if (!is.null(search_out)) { search_path_fulldata <- search_out[["search_path"]] } else { - verb_txt_search <- paste0("-----\nRunning ", method, " search ") + verb_txt_search <- "-----\nRunning the search " if (validate_search) { # Point out that this is the full-data search (if `validate_search` is # `FALSE`, this is still a full-data search, but in that case, there are # no fold-wise searches, so pointing out "full-data" could be confusing): verb_txt_search <- paste0(verb_txt_search, "using the full dataset ") } - # Note concerning the following verbose text: If `nclusters == S`, - # get_refdist() will use "thinning", not "clustering" (in that case, they - # give the same set of draws, namely the original one; hence the quotation - # marks), but here for this verbose message, we do not want to make things - # too complicated: - verb_txt_search <- paste0(verb_txt_search, "with ", - ifelse(!is.null(nclusters), - paste0(nclusters, " clusters "), - paste0(ndraws, " draws (from thinning) "))) verb_txt_search <- paste0(verb_txt_search, "...") verb_out(verb_txt_search, verbose = verbose) - search_path_fulldata <- .select( + search_path_fulldata <- select( refmodel = refmodel, ndraws = ndraws, nclusters = nclusters, method = method, nterms_max = nterms_max, penalty = penalty, verbose = verbose, search_control = search_control, @@ -766,15 +751,8 @@ loo_varsel <- function(refmodel, method, nterms_max, ndraws, # "Run" the performance evaluation for the submodels along the predictor # ranking (in fact, we only prepare the performance evaluation by computing # precursor quantities, but for users, this difference is not perceivable): - verb_out("-----\nRunning the performance evaluation with ", - ifelse(refit_prj, - ifelse(!is.null(nclusters_pred), - paste0(nclusters_pred, " clusters"), - paste0(ndraws_pred, " draws (from thinning)")), - ifelse(!is.null(nclusters), - paste0(nclusters, " clusters"), - paste0(ndraws, " draws (from thinning)"))), - " (`refit_prj = ", refit_prj, "`) ...", verbose = verbose) + verb_out("-----\nRunning the performance evaluation with `refit_prj = ", + refit_prj, "` ...", verbose = verbose) # Step 1: Re-project (using the full dataset) onto the submodels along the # full-data predictor ranking and evaluate their predictive performance. perf_eval_out <- perf_eval( @@ -824,53 +802,59 @@ loo_varsel <- function(refmodel, method, nterms_max, ndraws, )) } if (nrow(log_lik_ref) > 1) { - # Take into account that clustered draws usually have different weights: - lw_sub <- log_lik_ref + log(refdist_eval$wdraws_prj) - # This re-weighting requires a re-normalization (as.array() is applied to - # have stricter consistency checks, see `?sweep`): - lw_sub <- sweep(lw_sub, 2, as.array(apply(lw_sub, 2, log_sum_exp))) - # Internally, loo::psis() doesn't perform the Pareto smoothing if the - # number of draws is small (as indicated by object `no_psis_eval`, see - # below). In projpred, this can occur, e.g., if users request a number - # of projected draws (for performance evaluation, either after - # clustering or thinning the reference model's posterior draws) that is - # much smaller than the default of 400. In order to throw a customized - # warning message (and to avoid the calculation of Pareto k-values, see - # loo issue stan-dev/loo#227), object `no_psis_eval` indicates whether - # loo::psis() would perform the Pareto smoothing or not (for the - # decision rule, see loo:::n_pareto() and loo:::enough_tail_samples(), - # keeping in mind that we have `r_eff = 1` for all observations here). - S_for_psis_eval <- nrow(log_lik_ref) - no_psis_eval <- ceiling(min(0.2 * S_for_psis_eval, - 3 * sqrt(S_for_psis_eval))) < 5 - if (no_psis_eval) { + # Use loo::sis() if the projected draws (i.e., the draws resulting + # from the clustering or thinning) have nonconstant weights: + if (refdist_eval$const_wdraws_prj) { + # Internally, loo::psis() doesn't perform the Pareto smoothing if the + # number of draws is small (as indicated by object `no_psis_eval`, see + # below). In projpred, this can occur, e.g., if users request a number + # of projected draws (for performance evaluation, either after + # clustering or thinning the reference model's posterior draws) that is + # much smaller than the default of 400. In order to throw a customized + # warning message (and to avoid the calculation of Pareto k-values, see + # loo issue stan-dev/loo#227), object `no_psis_eval` indicates whether + # loo::psis() would perform the Pareto smoothing or not (for the + # decision rule, see loo:::n_pareto() and loo:::enough_tail_samples(), + # keeping in mind that we have `r_eff = 1` for all observations here). + S_for_psis_eval <- nrow(log_lik_ref) + no_psis_eval <- ceiling(min(0.2 * S_for_psis_eval, + 3 * sqrt(S_for_psis_eval))) < 5 + if (no_psis_eval) { + if (getOption("projpred.warn_psis", TRUE)) { + warning( + "Using standard importance sampling (SIS), as the number of ", + "draws or clusters is too small for PSIS. For improved ", + "accuracy, increase the number of draws or clusters, or use ", + "K-fold-CV." + ) + } + # Use loo::sis(). + # In principle, we could rely on loo::psis() here (because in such a + # case, it would internally switch to SIS automatically), but using + # loo::sis() explicitly is safer because if the loo package changes + # its decision rule, we would get a mismatch between our customized + # warning here and the IS method used by loo. See also loo issue + # stan-dev/loo#227. + importance_sampling_nm <- "sis" + } else { + # Use loo::psis(). + # Usually, we have a small number of projected draws here (400 by + # default), which means that the 'loo' package will automatically + # perform the regularization from Vehtari et al. (2024, + # , appendix G). + importance_sampling_nm <- "psis" + } + } else { if (getOption("projpred.warn_psis", TRUE)) { - message( - "Using standard importance sampling (SIS) due to a small number of", - ifelse(refit_prj, - ifelse(!is.null(nclusters_pred), - " clusters", - " draws (from thinning)"), - ifelse(!is.null(nclusters), - " clusters", - " draws (from thinning)")) + warning( + "The projected draws used for the performance evaluation have ", + "different (i.e., nonconstant) weights, so using standard ", + "importance sampling (SIS) instead of Pareto-smoothed importance ", + "sampling (PSIS). In general, PSIS is recommended over SIS." ) } # Use loo::sis(). - # In principle, we could rely on loo::psis() here (because in such a - # case, it would internally switch to SIS automatically), but using - # loo::sis() explicitly is safer because if the loo package changes - # its decision rule, we would get a mismatch between our customized - # warning here and the IS method used by loo. See also loo issue - # stan-dev/loo#227. importance_sampling_nm <- "sis" - } else { - # Use loo::psis(). - # Usually, we have a small number of projected draws here (400 by - # default), which means that the 'loo' package will automatically - # perform the regularization from Vehtari et al. (2022, - # , appendix G). - importance_sampling_nm <- "psis" } importance_sampling_func <- get(importance_sampling_nm, asNamespace("loo")) @@ -907,6 +891,11 @@ loo_varsel <- function(refmodel, method, nterms_max, ndraws, } else { lw_sub <- matrix(0, nrow = nrow(log_lik_ref), ncol = ncol(log_lik_ref)) } + # Take into account that clustered draws usually have different weights: + lw_sub <- lw_sub + log(refdist_eval$wdraws_prj) + # This re-weighting requires a re-normalization (as.array() is applied to + # have stricter consistency checks, see `?sweep`): + lw_sub <- sweep(lw_sub, 2, as.array(apply(lw_sub, 2, log_sum_exp))) for (k in seq_len(1 + length(search_path_fulldata$predictor_ranking))) { # TODO: For consistency, replace `k` in this `for` loop by `j`. mu_k <- perf_eval_out[["mu_by_size"]][[k]] @@ -974,23 +963,14 @@ loo_varsel <- function(refmodel, method, nterms_max, ndraws, } if (verbose) { - verb_out("-----\nRunning ", - ifelse(!search_out_rks_was_null, "", - paste0(method, " the search with ", - ifelse(!is.null(nclusters), - paste0(nclusters, " clusters"), - paste0(ndraws, " draws (from thinning)")), - " and ")), - "the performance evaluation with ", - ifelse(refit_prj, - ifelse(!is.null(nclusters_pred), - paste0(nclusters_pred, " clusters"), - paste0(ndraws_pred, " draws (from thinning)")), - ifelse(!is.null(nclusters), - paste0(nclusters, " clusters"), - paste0(ndraws, " draws (from thinning)"))), - " (`refit_prj = ", refit_prj, - "`) for each of the `nloo = ", nloo, "` ", + verb_txt_start <- "-----\nRunning " + if (!search_out_rks_was_null) { + verb_txt_mid <- "" + } else { + verb_txt_mid <- "the search and " + } + verb_out(verb_txt_start, verb_txt_mid, "the performance evaluation with ", + "`refit_prj = ", refit_prj, "` for each of the N = ", nloo, " ", "LOO-CV folds separately ...") } one_obs <- function(run_index, @@ -1007,7 +987,7 @@ loo_varsel <- function(refmodel, method, nterms_max, ndraws, if (!search_out_rks_was_null) { search_path <- list(predictor_ranking = search_out_rks[[run_index]]) } else { - search_path <- .select( + search_path <- select( refmodel = refmodel, ndraws = ndraws, nclusters = nclusters, reweighting_args = list(cl_ref = cl_sel, wdraws_ref = exp(lw[, i])), method = method, nterms_max = nterms_max, penalty = penalty, @@ -1322,7 +1302,7 @@ kfold_varsel <- function(refmodel, method, nterms_max, ndraws, nclusters, } else if (!search_out_rks_was_null) { search_path <- list(predictor_ranking = rk) } else { - search_path <- .select( + search_path <- select( refmodel = fold$refmodel, ndraws = ndraws, nclusters = nclusters, method = method, nterms_max = nterms_max, penalty = penalty, verbose = verbose_search, search_control = search_control, diff --git a/R/glmfun.R b/R/glmfun.R index 044678b9a..648b25348 100644 --- a/R/glmfun.R +++ b/R/glmfun.R @@ -16,7 +16,7 @@ standardization <- function(x, center = TRUE, scale = TRUE, weights = NULL) { mx <- rep(0, ncol(x)) } if (scale) { - sx <- apply(x, 2, .weighted_sd, w) + sx <- apply(x, 2, weighted.sd, w) } else { sx <- rep(1, ncol(x)) } diff --git a/R/misc.R b/R/misc.R index 20f7d75ad..b4e8fb451 100644 --- a/R/misc.R +++ b/R/misc.R @@ -1,8 +1,8 @@ .onAttach <- function(...) { ver <- utils::packageVersion("projpred") msg <- paste0("This is projpred version ", ver, ".") - msg <- paste0(msg, "\n", "NOTE: In projpred 2.7.0, the default search ", - "method was set to \"forward\" (for all kinds of models).") + msg <- paste0(msg, " ", "NOTE: In projpred 2.7.0, the default search method ", + "was set to \"forward\" (for all kinds of models).") packageStartupMessage(msg) } @@ -14,7 +14,7 @@ nms_y_wobs_test <- function(wobs_nm = "wobs") { c("y", "y_oscale", wobs_nm) } -.weighted_sd <- function(x, w, na.rm = FALSE) { +weighted.sd <- function(x, w, na.rm = FALSE) { if (length(w) == 1 && length(x) > 1) { w <- rep_len(w, length.out = length(x)) } @@ -66,7 +66,7 @@ ilinkfun_raw <- function(x, link_nm) { return(basic_ilink(x)) } -.auc <- function(x) { +auc <- function(x) { resp <- x[, 1] pred <- x[, 2] wobs <- x[, 3] diff --git a/R/summary_funs.R b/R/summary_funs.R index 8903e373c..e15256825 100644 --- a/R/summary_funs.R +++ b/R/summary_funs.R @@ -343,7 +343,7 @@ get_stat <- function(summaries, summaries_baseline = NULL, } else { # full LOO estimator value <- mean(wobs * (mu - y)^2) - value_se <- .weighted_sd((mu - y)^2, wobs) / sqrt(n_full) + value_se <- weighted.sd((mu - y)^2, wobs) / sqrt(n_full) } # store for later calculations mse_e <- value @@ -353,7 +353,7 @@ get_stat <- function(summaries, summaries_baseline = NULL, # limit theorem would ensure convergence -- probably slower, though -- to # a normal distribution even for MSE and RMSE) mse_b <- mean(wobs * (mu_baseline - y)^2) - var_mse_b <- .weighted_sd((mu_baseline - y)^2, wobs)^2 / n_full + var_mse_b <- weighted.sd((mu_baseline - y)^2, wobs)^2 / n_full if (n_loo < n_full) { mse_e_fast <- mean(wobs * (summaries_fast$mu - y)^2) srs_diffe <- @@ -391,7 +391,7 @@ get_stat <- function(summaries, summaries_baseline = NULL, mse_y <- mean(wobs * (y_mean_w - y)^2) value <- 1 - mse_e / mse_y - ifelse(is.null(summaries_baseline), 0, 1 - mse_b / mse_y) # the first-order Taylor approximation of the variance - var_mse_y <- .weighted_sd((y_mean_w - y)^2, wobs)^2 / n_full + var_mse_y <- weighted.sd((y_mean_w - y)^2, wobs)^2 / n_full if (n_loo < n_full) { mse_e_fast <- mean(wobs * (summaries_fast$mu - y)^2) if (is.null(summaries_baseline)) { @@ -524,7 +524,7 @@ get_stat <- function(summaries, summaries_baseline = NULL, } else { # full LOO estimator value <- mean(wobs * correct) - mean(wobs * correct_baseline) - value_se <- .weighted_sd(correct - correct_baseline, wobs) / sqrt(n_full) + value_se <- weighted.sd(correct - correct_baseline, wobs) / sqrt(n_full) } } else if (stat == "auc") { if (n_loo < n_full) { @@ -535,15 +535,15 @@ get_stat <- function(summaries, summaries_baseline = NULL, if (!is.null(mu_baseline)) { auc_data <- cbind(y, mu, wobs) auc_data_baseline <- cbind(y, mu_baseline, wobs) - value <- .auc(auc_data) - .auc(auc_data_baseline) + value <- auc(auc_data) - auc(auc_data_baseline) idxs_cols <- seq_len(ncol(auc_data)) idxs_cols_bs <- setdiff(seq_len(ncol(auc_data) + ncol(auc_data_baseline)), idxs_cols) diffvalue.bootstrap <- bootstrap( cbind(auc_data, auc_data_baseline), function(x) { - .auc(x[, idxs_cols, drop = FALSE]) - - .auc(x[, idxs_cols_bs, drop = FALSE]) + auc(x[, idxs_cols, drop = FALSE]) - + auc(x[, idxs_cols_bs, drop = FALSE]) }, ... ) @@ -558,8 +558,8 @@ get_stat <- function(summaries, summaries_baseline = NULL, } } else { auc_data <- cbind(y, mu, wobs) - value <- .auc(auc_data) - value.bootstrap <- bootstrap(auc_data, .auc, ...) + value <- auc(auc_data) + value.bootstrap <- bootstrap(auc_data, auc, ...) value_se <- sd(value.bootstrap) if (any(is.na(value.bootstrap))) { # quantile() is not able to deal with `NA`s @@ -654,7 +654,7 @@ get_nfeat_baseline <- function(object, baseline, stat, ...) { # eq (7) est_list$y_hat <- t_pi_tilde + t_e # eq (8) - var_e_i <- .weighted_sd(e_i, w = wobs_m)^2 + var_e_i <- weighted.sd(e_i, w = wobs_m)^2 est_list$v_y_hat <- N^2 * (1 - m / N) * var_e_i / m # eq (9) first row second `+` should be `-` # Supplementary material eq (6) has this correct diff --git a/R/varsel.R b/R/varsel.R index 1f7b5632a..3fadef546 100644 --- a/R/varsel.R +++ b/R/varsel.R @@ -248,27 +248,15 @@ varsel.vsel <- function(object, ...) { #' @rdname varsel #' @export -varsel.refmodel <- function( - object, - d_test = NULL, - method = "forward", - ndraws = NULL, - nclusters = 20, - ndraws_pred = 400, - nclusters_pred = NULL, - refit_prj = !inherits(object, "datafit"), - nterms_max = NULL, - verbose = getOption("projpred.verbose", interactive()), - search_control = NULL, - lambda_min_ratio = 1e-5, - nlambda = 150, - thresh = 1e-6, - penalty = NULL, - search_terms = NULL, - search_out = NULL, - seed = NA, - ... -) { +varsel.refmodel <- function(object, d_test = NULL, method = "forward", + ndraws = NULL, nclusters = 20, ndraws_pred = 400, + nclusters_pred = NULL, + refit_prj = !inherits(object, "datafit"), + nterms_max = NULL, verbose = TRUE, + search_control = NULL, lambda_min_ratio = 1e-5, + nlambda = 150, thresh = 1e-6, penalty = NULL, + search_terms = NULL, search_out = NULL, seed = NA, + ...) { if (!missing(lambda_min_ratio)) { warning("Argument `lambda_min_ratio` is deprecated. Please specify ", "control arguments for the search via argument `search_control`. ", @@ -388,7 +376,7 @@ varsel.refmodel <- function( search_path <- search_out[["search_path"]] } else { verb_out("-----\nRunning the search ...", verbose = verbose) - search_path <- .select( + search_path <- select( refmodel = refmodel, ndraws = ndraws, nclusters = nclusters, method = method, nterms_max = nterms_max, penalty = penalty, verbose = verbose, search_control = search_control, @@ -523,8 +511,8 @@ varsel.refmodel <- function( # `outdmins` (the submodel fits along the predictor ranking, with the number # of fits per model size being equal to the number of projected draws), and # `p_sel` (the output from get_refdist() for the search). -.select <- function(refmodel, ndraws, nclusters, reweighting_args = NULL, - method, nterms_max, penalty, verbose, search_control, ...) { +select <- function(refmodel, ndraws, nclusters, reweighting_args = NULL, method, + nterms_max, penalty, verbose, search_control, ...) { if (is.null(reweighting_args)) { p_sel <- get_refdist(refmodel, ndraws = ndraws, nclusters = nclusters) } else { diff --git a/man/cv_varsel.Rd b/man/cv_varsel.Rd index eafb9ab63..46196c9b5 100644 --- a/man/cv_varsel.Rd +++ b/man/cv_varsel.Rd @@ -32,7 +32,7 @@ cv_varsel(object, ...) refit_prj = !inherits(object, "datafit"), nterms_max = NULL, penalty = NULL, - verbose = getOption("projpred.verbose", interactive()), + verbose = TRUE, nloo = if (cv_method == "LOO") object$nobs else NULL, K = if (!inherits(object, "datafit")) 5 else 10, cvfits = object$cvfits, diff --git a/man/varsel.Rd b/man/varsel.Rd index f0cc84937..b9ade6fbb 100644 --- a/man/varsel.Rd +++ b/man/varsel.Rd @@ -23,7 +23,7 @@ varsel(object, ...) nclusters_pred = NULL, refit_prj = !inherits(object, "datafit"), nterms_max = NULL, - verbose = getOption("projpred.verbose", interactive()), + verbose = TRUE, search_control = NULL, lambda_min_ratio = 1e-05, nlambda = 150, diff --git a/tests/testthat/helpers/testers.R b/tests/testthat/helpers/testers.R index 4d2422fe9..b746b8f83 100644 --- a/tests/testthat/helpers/testers.R +++ b/tests/testthat/helpers/testers.R @@ -1980,8 +1980,7 @@ vsel_tester <- function( search_control_expected = NULL, extra_tol = 1.1, info_str = "" - ) { - +) { # Preparations: if (with_cv) { if (is.null(cv_method_expected)) { @@ -2308,7 +2307,6 @@ vsel_tester <- function( } return(invisible(TRUE)) } - for (j in seq_along(vs$summaries$sub)) { smmrs_sub_j_tester(vs$summaries$sub[[j]]) if (vs$refmodel$family$for_latent) { diff --git a/tests/testthat/test_misc.R b/tests/testthat/test_misc.R index ecbfbc2c8..66758b964 100644 --- a/tests/testthat/test_misc.R +++ b/tests/testthat/test_misc.R @@ -273,142 +273,142 @@ test_that(".srs_diff_est_w() propagates input `NA`s to its output", { ) }) -test_that(".weighted_sd() with `na.rm = FALSE` propagates input `NA`s to its output", { +test_that("weighted.sd() with `na.rm = FALSE` propagates input `NA`s to its output", { expect_identical( - .weighted_sd(x = rep(NA, nobsv), - w = rep(NA, nobsv)), + weighted.sd(x = rep(NA, nobsv), + w = rep(NA, nobsv)), NA_real_, info = "all inputs are all-`NA`s" ) expect_identical( - .weighted_sd(x = rep(-0.8, nobsv), - w = c(rep_len(c(2, 3), length.out = nobsv - 1L), NA)), + weighted.sd(x = rep(-0.8, nobsv), + w = c(rep_len(c(2, 3), length.out = nobsv - 1L), NA)), NA_real_, info = "`x` without `NA`s, `w` with a single `NA`" ) expect_identical( - .weighted_sd(x = c(rep(-0.8, nobsv - 1L), NA), - w = rep_len(c(2, 3), length.out = nobsv)), + weighted.sd(x = c(rep(-0.8, nobsv - 1L), NA), + w = rep_len(c(2, 3), length.out = nobsv)), NA_real_, info = "`x` with a single `NA`, `w` without `NA`s" ) expect_identical( - .weighted_sd(x = rep(-0.8, nobsv), - w = rep(NA, nobsv)), + weighted.sd(x = rep(-0.8, nobsv), + w = rep(NA, nobsv)), NA_real_, info = "`x` without `NA`s, `w` all-`NA`s" ) expect_identical( - .weighted_sd(x = rep(NA, nobsv), - w = rep_len(c(2, 3), length.out = nobsv)), + weighted.sd(x = rep(NA, nobsv), + w = rep_len(c(2, 3), length.out = nobsv)), NA_real_, info = "`x` all-`NA`s, `w` without `NA`s" ) }) -test_that(".auc() works", { +test_that("auc() works", { nobsv_auc <- 19L expect_equal( - .auc(cbind(rep_len(c(0, 1), length.out = nobsv_auc), - imperfect_alternation(c(0.3, 0.8), length.out = nobsv_auc), - rep(1, nobsv_auc))), + auc(cbind(rep_len(c(0, 1), length.out = nobsv_auc), + imperfect_alternation(c(0.3, 0.8), length.out = nobsv_auc), + rep(1, nobsv_auc))), 0.6833333333333333333333, info = "`wobs` column only `1`s" ) expect_equal( - .auc(cbind(rep_len(c(0, 1), length.out = nobsv_auc), - imperfect_alternation(c(0.3, 0.8), length.out = nobsv_auc), - imperfect_alternation(c(1, 2), n_tail = 11L, length.out = nobsv_auc))), + auc(cbind(rep_len(c(0, 1), length.out = nobsv_auc), + imperfect_alternation(c(0.3, 0.8), length.out = nobsv_auc), + imperfect_alternation(c(1, 2), n_tail = 11L, length.out = nobsv_auc))), 0.717948717948718062587, info = "`wobs` column with imperfect alternation of `1`s and `2`s" ) }) -test_that(".auc() propagates input `NA`s to its output", { +test_that("auc() propagates input `NA`s to its output", { nobsv_auc <- 19L expect_identical( - .auc(cbind(rep(NA, nobsv_auc), - rep(NA, nobsv_auc), - rep(NA, nobsv_auc))), + auc(cbind(rep(NA, nobsv_auc), + rep(NA, nobsv_auc), + rep(NA, nobsv_auc))), NA_real_, info = "all inputs are all-`NA`s" ) expect_identical( - .auc(cbind(rep(1, nobsv_auc), - rep(NA, nobsv_auc), - rep(1, nobsv_auc))), + auc(cbind(rep(1, nobsv_auc), + rep(NA, nobsv_auc), + rep(1, nobsv_auc))), NA_real_, info = paste0("`resp` column without `NA`s, ", "`pred` column only `NA`s, ", "`wobs` column without `NA`s") ) expect_identical( - .auc(cbind(rep(1, nobsv_auc), - rep(NA, nobsv_auc), - c(rep(1, nobsv_auc - 1L), NA))), + auc(cbind(rep(1, nobsv_auc), + rep(NA, nobsv_auc), + c(rep(1, nobsv_auc - 1L), NA))), NA_real_, info = paste0("`resp` column without `NA`s, ", "`pred` column only `NA`s, ", "`wobs` column with a single `NA`") ) expect_identical( - .auc(cbind(c(rep_len(c(0, 1), length.out = nobsv_auc - 1L), NA), - rep(0.7, nobsv_auc), - rep(1, nobsv_auc))), + auc(cbind(c(rep_len(c(0, 1), length.out = nobsv_auc - 1L), NA), + rep(0.7, nobsv_auc), + rep(1, nobsv_auc))), NA_real_, info = paste0("`resp` column with a single `NA`, ", "`pred` column without `NA`s, ", "`wobs` column without `NA`s") ) expect_identical( - .auc(cbind(rep(NA, nobsv_auc), - rep(0.7, nobsv_auc), - rep(1, nobsv_auc))), + auc(cbind(rep(NA, nobsv_auc), + rep(0.7, nobsv_auc), + rep(1, nobsv_auc))), NA_real_, info = paste0("`resp` column only `NA`s, ", "`pred` column without `NA`s, ", "`wobs` column without `NA`s") ) expect_identical( - .auc(cbind(c(rep_len(c(0, 1), length.out = nobsv_auc - 1L), NA), - c(rep(0.7, nobsv_auc - 1L), NA), - rep(1, nobsv_auc))), + auc(cbind(c(rep_len(c(0, 1), length.out = nobsv_auc - 1L), NA), + c(rep(0.7, nobsv_auc - 1L), NA), + rep(1, nobsv_auc))), NA_real_, info = paste0("`resp` column with a single `NA`, ", "`pred` column with a single `NA`, ", "`wobs` column without `NA`s") ) expect_identical( - .auc(cbind(c(rep_len(c(0, 1), length.out = nobsv_auc - 1L), NA), - c(NA, rep(0.7, nobsv_auc - 1L)), - rep(1, nobsv_auc))), + auc(cbind(c(rep_len(c(0, 1), length.out = nobsv_auc - 1L), NA), + c(NA, rep(0.7, nobsv_auc - 1L)), + rep(1, nobsv_auc))), NA_real_, info = paste0("`resp` column with a single `NA`, ", "`pred` column with a single `NA` (at different position), ", "`wobs` column without `NA`s") ) expect_identical( - .auc(cbind(c(rep_len(c(0, 1), length.out = nobsv_auc - 1L), NA), - rep(NA, nobsv_auc), - rep(1, nobsv_auc))), + auc(cbind(c(rep_len(c(0, 1), length.out = nobsv_auc - 1L), NA), + rep(NA, nobsv_auc), + rep(1, nobsv_auc))), NA_real_, info = paste0("`resp` column with a single `NA`, ", "`pred` column only `NA`s, ", "`wobs` column without `NA`s") ) expect_identical( - .auc(cbind(rep(NA, nobsv_auc), - rep(0.7, nobsv_auc), - c(rep(1, nobsv_auc - 1L), NA))), + auc(cbind(rep(NA, nobsv_auc), + rep(0.7, nobsv_auc), + c(rep(1, nobsv_auc - 1L), NA))), NA_real_, info = paste0("`resp` column only `NA`s, ", "`pred` column without `NA`s, ", "`wobs` column with a single `NA`") ) expect_identical( - .auc(cbind(rep_len(c(0, 1), length.out = nobsv_auc), - imperfect_alternation(c(0.3, 0.8), length.out = nobsv_auc), - rep(NA, nobsv_auc))), + auc(cbind(rep_len(c(0, 1), length.out = nobsv_auc), + imperfect_alternation(c(0.3, 0.8), length.out = nobsv_auc), + rep(NA, nobsv_auc))), NA_real_, info = paste0("`resp` column without `NA`s, ", "`pred` column without `NA`s, ", From 7bef2c2c1af26e36cc66bfba5a074a39f1b92288 Mon Sep 17 00:00:00 2001 From: fweber144 Date: Mon, 10 Feb 2025 15:02:41 +0100 Subject: [PATCH 37/40] document use of the progressr package, in particular global option `projpred.use_progressr` --- R/projpred-package.R | 12 ++++++++++++ man/projpred-package.Rd | 11 +++++++++++ 2 files changed, 23 insertions(+) diff --git a/R/projpred-package.R b/R/projpred-package.R index 7154adf68..e376a6ac7 100644 --- a/R/projpred-package.R +++ b/R/projpred-package.R @@ -135,6 +135,18 @@ #' (i.e., a parallelization of \pkg{projpred}'s cross-validation) which can be #' activated via argument `parallel`. #' +#' During parallelization (either of the projection or the CV), progression +#' updates can be received via the \pkg{progressr} package. This only works if +#' the \pkg{doFuture} backend is used for parallelization, e.g., via +#' `doFuture::registerDoFuture()` and `future::plan(future::multisession, +#' workers = 4)`. In that case, the \pkg{progressr} package can be used, e.g., +#' by calling `progressr::handlers(global = TRUE)` before running the projection +#' or the CV in parallel. The \pkg{projpred} package also offers the global +#' option `projpred.use_progressr` for controlling whether to use the +#' \pkg{progressr} package (`TRUE` or `FALSE`), but since that global option +#' defaults to `requireNamespace("progressr", quietly = TRUE) && interactive()`, +#' it usually does not need to be set by the user. +#' #' # Multilevel models: "Integrating out" group-level effects #' #' In case of multilevel models, \pkg{projpred} offers two global options for diff --git a/man/projpred-package.Rd b/man/projpred-package.Rd index b4f0627fe..c9d929bb5 100644 --- a/man/projpred-package.Rd +++ b/man/projpred-package.Rd @@ -131,6 +131,17 @@ without additive terms may be regarded as a projection onto a submodel which is a GLM). However, for \code{\link[=cv_varsel]{cv_varsel()}}, there is also a \emph{CV} parallelization (i.e., a parallelization of \pkg{projpred}'s cross-validation) which can be activated via argument \code{parallel}. + +During parallelization (either of the projection or the CV), progression +updates can be received via the \pkg{progressr} package. This only works if +the \pkg{doFuture} backend is used for parallelization, e.g., via +\code{doFuture::registerDoFuture()} and \code{future::plan(future::multisession, workers = 4)}. In that case, the \pkg{progressr} package can be used, e.g., +by calling \code{progressr::handlers(global = TRUE)} before running the projection +or the CV in parallel. The \pkg{projpred} package also offers the global +option \code{projpred.use_progressr} for controlling whether to use the +\pkg{progressr} package (\code{TRUE} or \code{FALSE}), but since that global option +defaults to \code{requireNamespace("progressr", quietly = TRUE) && interactive()}, +it usually does not need to be set by the user. } \section{Multilevel models: "Integrating out" group-level effects}{ From d512dba3d9afa2d00673abc8feed7a729a056143 Mon Sep 17 00:00:00 2001 From: fweber144 Date: Mon, 10 Feb 2025 21:25:10 +0100 Subject: [PATCH 38/40] add `identical(foreach::getDoParName(), "doFuture")` to the default for global option `projpred.use_progressr` --- R/misc.R | 4 +++- R/projpred-package.R | 5 +++-- man/projpred-package.Rd | 4 ++-- 3 files changed, 8 insertions(+), 5 deletions(-) diff --git a/R/misc.R b/R/misc.R index b4e8fb451..32b1b5fad 100644 --- a/R/misc.R +++ b/R/misc.R @@ -716,5 +716,7 @@ element_unq <- function(list_obj, nm) { use_progressr <- function() { getOption("projpred.use_progressr", - requireNamespace("progressr", quietly = TRUE) && interactive()) + requireNamespace("progressr", quietly = TRUE) && + interactive() && + identical(foreach::getDoParName(), "doFuture")) } diff --git a/R/projpred-package.R b/R/projpred-package.R index e376a6ac7..62d01db19 100644 --- a/R/projpred-package.R +++ b/R/projpred-package.R @@ -144,8 +144,9 @@ #' or the CV in parallel. The \pkg{projpred} package also offers the global #' option `projpred.use_progressr` for controlling whether to use the #' \pkg{progressr} package (`TRUE` or `FALSE`), but since that global option -#' defaults to `requireNamespace("progressr", quietly = TRUE) && interactive()`, -#' it usually does not need to be set by the user. +#' defaults to `requireNamespace("progressr", quietly = TRUE) && interactive() +#' && identical(foreach::getDoParName(), "doFuture")`, it usually does not need +#' to be set by the user. #' #' # Multilevel models: "Integrating out" group-level effects #' diff --git a/man/projpred-package.Rd b/man/projpred-package.Rd index c9d929bb5..d6dfa1638 100644 --- a/man/projpred-package.Rd +++ b/man/projpred-package.Rd @@ -140,8 +140,8 @@ by calling \code{progressr::handlers(global = TRUE)} before running the projecti or the CV in parallel. The \pkg{projpred} package also offers the global option \code{projpred.use_progressr} for controlling whether to use the \pkg{progressr} package (\code{TRUE} or \code{FALSE}), but since that global option -defaults to \code{requireNamespace("progressr", quietly = TRUE) && interactive()}, -it usually does not need to be set by the user. +defaults to \code{requireNamespace("progressr", quietly = TRUE) && interactive() && identical(foreach::getDoParName(), "doFuture")}, it usually does not need +to be set by the user. } \section{Multilevel models: "Integrating out" group-level effects}{ From 2730bcc1da4765d408a12b96968630acb5a9d268 Mon Sep 17 00:00:00 2001 From: fweber144 Date: Mon, 10 Feb 2025 21:33:12 +0100 Subject: [PATCH 39/40] update progressor *after* executing the corresponding code, not before executing it --- R/cv_varsel.R | 14 +++++++++----- R/divergence_minimizers.R | 4 ++-- 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/R/cv_varsel.R b/R/cv_varsel.R index 9f09600d0..6a72daf3d 100644 --- a/R/cv_varsel.R +++ b/R/cv_varsel.R @@ -1051,9 +1051,11 @@ loo_varsel <- function(refmodel, method, nterms_max, ndraws, "loo_ref_oscale", "validset", "loo_sub", "mu_sub", "loo_sub_oscale", "mu_sub_oscale") ) %do_projpred% { + out_one_obs <- do.call(one_obs, c(list(run_index = run_index, + verbose_search = FALSE), + dot_args)) if (!is.null(progressor_obj)) progressor_obj() - do.call(one_obs, c(list(run_index = run_index, verbose_search = FALSE), - dot_args)) + return(out_one_obs) } } # For storing the fold-wise predictor rankings: @@ -1378,10 +1380,12 @@ kfold_varsel <- function(refmodel, method, nterms_max, ndraws, nclusters, .export = c("one_fold", "dot_args", "progressor_obj"), .noexport = c("list_cv", "search_out_rks") ) %do_projpred% { + out_one_fold <- do_call(one_fold, c(list(fold = list_cv_k, + rk = search_out_rks_k, + verbose_search = FALSE), + dot_args)) if (!is.null(progressor_obj)) progressor_obj() - do_call(one_fold, c(list(fold = list_cv_k, rk = search_out_rks_k, - verbose_search = FALSE), - dot_args)) + return(out_one_fold) } } verb_out("-----", verbose = verbose) diff --git a/R/divergence_minimizers.R b/R/divergence_minimizers.R index dbab5a6cd..5611d1869 100644 --- a/R/divergence_minimizers.R +++ b/R/divergence_minimizers.R @@ -109,7 +109,6 @@ divmin <- function( "projpred_var", "projpred_ws_aug", "projpred_formulas_no_random" ) ) %do_projpred% { - if (!is.null(progressor_obj)) progressor_obj() mssgs_warns_capt <- capt_mssgs_warns( soutdmin <- do.call( sdivmin, @@ -120,6 +119,7 @@ divmin <- function( dot_args) ) ) + if (!is.null(progressor_obj)) progressor_obj() return(nlist(soutdmin, mssgs_warns_capt)) } } @@ -675,7 +675,6 @@ divmin_augdat <- function( "projpred_ws_aug", "linkobjs" ) ) %do_projpred% { - if (!is.null(progressor_obj)) progressor_obj() mssgs_warns_capt <- capt_mssgs_warns( soutdmin <- do.call( sdivmin, @@ -688,6 +687,7 @@ divmin_augdat <- function( dot_args) ) ) + if (!is.null(progressor_obj)) progressor_obj() return(nlist(soutdmin, mssgs_warns_capt)) } } From 39b662ed4d90958c905057699b7a89d312ac99e6 Mon Sep 17 00:00:00 2001 From: fweber144 Date: Wed, 12 Feb 2025 11:21:12 +0100 Subject: [PATCH 40/40] add `NEWS.md` entry --- NEWS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/NEWS.md b/NEWS.md index 38bc390f8..23fa7658a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -14,6 +14,7 @@ If you read this from a place other than ). This threshold depends on the Monte Carlo sample size and is often close to the former fixed threshold of 0.7 (a short introduction may also be found in the [LOO glossary](https://mc-stan.org/loo/reference/loo-glossary.html)). Correspondingly, the former "secondary" threshold of 0.5 is not used anymore either. (GitHub: #490, #498) +* When using the **doFuture** backend for parallelization, progression updates can now be received via the **progressr** package, see `` ?`projpred-package` `` (section "Parallelization"). (GitHub: #504) # projpred 2.8.0