tests/testthat/test_project.R

In paasim/glmproj: Projection Predictive Feature Selection

context("project()")

# object, predictor_terms, ndraws / nclusters, nterms ---------------------

test_that(paste(
  "`object` of class `refmodel` and arguments `predictor_terms` and `ndraws`",
  "(or `nclusters`) work"
), {
  skip_if_not(run_prj)
  for (tstsetup in names(prjs)) {
    args_prj_i <- args_prj[[tstsetup]]
    ndr_ncl <- ndr_ncl_dtls(args_prj_i)
    projection_tester(prjs[[tstsetup]],
                      refmod_expected = refmods[[args_prj_i$tstsetup_ref]],
                      prd_trms_expected = args_prj_i$predictor_terms,
                      nprjdraws_expected = ndr_ncl$nprjdraws,
                      with_clusters = ndr_ncl$clust_used,
                      info_str = tstsetup)
  }
})

test_that("invalid `predictor_terms` warns or fails", {
  skip_if_not(run_prj)
  tstsetups <- grep("\\.glm\\.gauss.*\\.prd_trms_x\\.clust$", names(prjs),
                    value = TRUE)
  for (tstsetup in tstsetups) {
    args_prj_i <- args_prj[[tstsetup]]

    # Non-`vsel` object combined with `predictor_terms = NULL`:
    expect_error(
      do.call(project, c(
        list(object = refmods[[args_prj_i$tstsetup_ref]],
             predictor_terms = NULL),
        excl_nonargs(args_prj_i, nms_excl_add = "predictor_terms")
      )),
      paste("^Please provide an `object` of class `vsel` or use argument",
            "`predictor_terms`\\.$"),
      info = tstsetup
    )

    # Repeating `predictor_terms`:
    p_long <- do.call(project, c(
      list(object = refmods[[args_prj_i$tstsetup_ref]],
           predictor_terms = rep_len(args_prj_i$predictor_terms,
                                     length.out = 1e4)),
      excl_nonargs(args_prj_i, nms_excl_add = "predictor_terms")
    ))
    expect_equal(p_long, prjs[[tstsetup]], info = tstsetup)

    # Invalid type:
    for (prd_trms_crr in list(2, 1:3, list(prd_trms_x, prd_trms_x))) {
      tstsetup_crr <- paste(tstsetup, paste(prd_trms_crr, collapse = ","),
                            sep = "__")
      expect_error(
        do.call(project, c(
          list(object = refmods[[args_prj_i$tstsetup_ref]],
               predictor_terms = prd_trms_crr),
          excl_nonargs(args_prj_i, nms_excl_add = "predictor_terms")
        )),
        paste(
          "^is\\.null\\(predictor_terms\\) \\|\\|",
          "is\\.vector\\(predictor_terms, \"character\"\\) is not TRUE$"
        ),
        info = tstsetup_crr
      )
    }

    # Should be working, but result in a projection onto the intercept-only
    # submodel:
    for (prd_trms_crr in list("1",
                              c("some_dummy_string", "another_dummy_string"))) {
      tstsetup_crr <- paste(tstsetup, paste(prd_trms_crr, collapse = ","),
                            sep = "__")
      expect_warning(
        p <- do.call(project, c(
          list(object = refmods[[args_prj_i$tstsetup_ref]],
               predictor_terms = prd_trms_crr),
          excl_nonargs(args_prj_i, nms_excl_add = "predictor_terms")
        )),
        paste("The following element\\(s\\) of `predictor_terms` could not be",
              "found in the table of possible predictor terms"),
        info = tstsetup_crr
      )
      projection_tester(p,
                        refmod_expected = refmods[[args_prj_i$tstsetup_ref]],
                        prd_trms_expected = character(),
                        nprjdraws_expected = nclusters_pred_tst,
                        with_clusters = TRUE,
                        info_str = tstsetup_crr)
    }
  }
})

test_that("`object` of class `stanreg` or `brmsfit` works", {
  skip_if_not(run_prj)
  tstsetups <- grep("\\.brnll\\..*\\.prd_trms_x\\.clust$", names(prjs),
                    value = TRUE)
  for (tstsetup in tstsetups) {
    args_prj_i <- args_prj[[tstsetup]]
    p_fit <- do.call(project, c(
      list(object = fits[[args_prj_i$tstsetup_fit]]),
      excl_nonargs(args_prj_i),
      excl_nonargs(args_ref[[args_prj_i$tstsetup_ref]])
    ))
    expect_identical(p_fit, prjs[[tstsetup]], ignore.environment = TRUE,
                     info = tstsetup)
  }
})

test_that(paste(
  "`object` of class `vsel` (created by varsel()) and arguments `nclusters`",
  "and `nterms` work"
), {
  skip_if_not(run_vs)
  for (tstsetup in names(prjs_vs)) {
    tstsetup_vs <- args_prj_vs[[tstsetup]]$tstsetup_vsel
    stopifnot(length(tstsetup_vs) > 0)
    mod_crr <- args_vs[[tstsetup_vs]]$mod_nm
    fam_crr <- args_vs[[tstsetup_vs]]$fam_nm
    nterms_crr <- args_prj_vs[[tstsetup]]$nterms
    if (is.null(nterms_crr)) {
      nterms_crr <- suggest_size(vss[[tstsetup_vs]], warnings = FALSE)
    }
    if (length(nterms_crr) == 1) {
      prd_trms_expected_crr <- vss[[tstsetup_vs]]$predictor_ranking[
        seq_len(nterms_crr)
      ]
      projection_tester(
        prjs_vs[[tstsetup]],
        refmod_expected = refmods[[args_prj_vs[[tstsetup]]$tstsetup_ref]],
        prd_trms_expected = prd_trms_expected_crr,
        nprjdraws_expected = args_prj_vs[[tstsetup]]$nclusters,
        with_clusters = TRUE,
        info_str = tstsetup
      )
      # Check that projecting from the `vsel` object onto a single submodel
      # gives the same output as projecting the reference model onto that
      # submodel directly:
      tstsetup_tries <- grep(
        paste0("^",
               gsub("\\.", "\\\\.",
                    sub("(with_offs|without_offs).*", "\\1", tstsetup)),
               ".*\\.clust$"),
        names(prjs),
        value = TRUE
      )
      match_prj <- sapply(tstsetup_tries, function(tstsetup_try) {
        setequal(prd_trms_expected_crr, prjs[[tstsetup_try]]$predictor_terms) &&
          args_prj_vs[[tstsetup]]$prj_nm == args_prj[[tstsetup_try]]$prj_nm
      })
      tstsetup_match_prj <- tstsetup_tries[match_prj]
      if (length(tstsetup_match_prj) == 1) {
        if (identical(prjs_vs[[tstsetup]]$predictor_terms,
                      prjs[[tstsetup_match_prj]]$predictor_terms)) {
          expect_identical(prjs_vs[[tstsetup]], prjs[[tstsetup_match_prj]],
                           ignore.environment = TRUE, info = tstsetup)
        } else {
          expect_setequal(prjs_vs[[tstsetup]]$predictor_terms,
                          prjs[[tstsetup_match_prj]]$predictor_terms)
          expect_equal(
            lapply(seq_along(prjs_vs[[tstsetup]]$outdmin), function(s_idx) {
              outdmin_s <- prjs_vs[[tstsetup]]$outdmin[[s_idx]]
              outdmin_s$beta <- outdmin_s$beta[
                rownames(prjs[[tstsetup_match_prj]]$outdmin[[s_idx]]$beta), ,
                drop = FALSE
              ]
              outdmin_s$formula <- prjs[[tstsetup_match_prj]]$outdmin[[s_idx]][[
                "formula"
              ]]
              outdmin_s$x <- outdmin_s$x[
                , colnames(prjs[[tstsetup_match_prj]]$outdmin[[s_idx]]$x),
                drop = FALSE
              ]
              return(outdmin_s)
            }),
            prjs[[tstsetup_match_prj]]$outdmin,
            tolerance = 1e1 * .Machine$double.eps, info = tstsetup
          )
          prj_nms <- names(prjs_vs[[tstsetup]])
          expect_identical(prj_nms, names(prjs[[tstsetup_match_prj]]),
                           info = tstsetup)
          prj_el_excl <- !prj_nms %in% c("predictor_terms", "outdmin")
          expect_equal(prjs_vs[[tstsetup]][prj_el_excl],
                       prjs[[tstsetup_match_prj]][prj_el_excl],
                       tolerance = .Machine$double.eps, info = tstsetup)
        }
      } else if (length(tstsetup_match_prj) > 1) {
        stop("Unexpected number of matches.")
      }
    } else {
      proj_list_tester(
        prjs_vs[[tstsetup]],
        len_expected = length(nterms_crr),
        is_seq = all(diff(nterms_crr) == 1),
        info_str = tstsetup,
        refmod_expected = refmods[[args_prj_vs[[tstsetup]]$tstsetup_ref]],
        nprjdraws_expected = args_prj_vs[[tstsetup]]$nclusters,
        with_clusters = TRUE,
        prjdraw_weights_expected = prjs_vs[[tstsetup]][[1]]$wdraws_prj
      )
    }
  }
})

test_that(paste(
  "`object` of class `vsel` (created by cv_varsel()) and arguments",
  "`nclusters` and `nterms` work"
), {
  skip_if_not(run_cvvs)
  for (tstsetup in names(prjs_cvvs)) {
    tstsetup_cvvs <- args_prj_cvvs[[tstsetup]]$tstsetup_vsel
    stopifnot(length(tstsetup_cvvs) > 0)
    mod_crr <- args_cvvs[[tstsetup_cvvs]]$mod_nm
    fam_crr <- args_cvvs[[tstsetup_cvvs]]$fam_nm
    nterms_crr <- args_prj_cvvs[[tstsetup]]$nterms
    if (is.null(nterms_crr)) {
      nterms_crr <- suggest_size(cvvss[[tstsetup_cvvs]], warnings = FALSE)
    }
    if (length(nterms_crr) == 1) {
      prd_trms_expected_crr <- cvvss[[tstsetup_cvvs]]$predictor_ranking[
        seq_len(nterms_crr)
      ]
      projection_tester(
        prjs_cvvs[[tstsetup]],
        refmod_expected = refmods[[args_prj_cvvs[[tstsetup]]$tstsetup_ref]],
        prd_trms_expected = prd_trms_expected_crr,
        nprjdraws_expected = args_prj_cvvs[[tstsetup]]$nclusters,
        with_clusters = TRUE,
        info_str = tstsetup
      )
      # Check that projecting from the `vsel` object onto a single submodel
      # gives the same output as projecting the reference model onto that
      # submodel directly:
      tstsetup_tries <- grep(
        paste0("^",
               gsub("\\.", "\\\\.",
                    sub("(with_offs|without_offs).*", "\\1", tstsetup)),
               ".*\\.clust$"),
        names(prjs),
        value = TRUE
      )
      match_prj <- sapply(tstsetup_tries, function(tstsetup_try) {
        setequal(prd_trms_expected_crr, prjs[[tstsetup_try]]$predictor_terms) &&
          args_prj_cvvs[[tstsetup]]$prj_nm == args_prj[[tstsetup_try]]$prj_nm
      })
      tstsetup_match_prj <- tstsetup_tries[match_prj]
      if (length(tstsetup_match_prj) == 1) {
        if (identical(prjs_cvvs[[tstsetup]]$predictor_terms,
                      prjs[[tstsetup_match_prj]]$predictor_terms)) {
          expect_identical(prjs_cvvs[[tstsetup]], prjs[[tstsetup_match_prj]],
                           ignore.environment = TRUE, info = tstsetup)
        } else {
          expect_setequal(prjs_cvvs[[tstsetup]]$predictor_terms,
                          prjs[[tstsetup_match_prj]]$predictor_terms)
          expect_equal(
            lapply(seq_along(prjs_cvvs[[tstsetup]]$outdmin), function(s_idx) {
              outdmin_s <- prjs_cvvs[[tstsetup]]$outdmin[[s_idx]]
              outdmin_s$beta <- outdmin_s$beta[
                rownames(prjs[[tstsetup_match_prj]]$outdmin[[s_idx]]$beta), ,
                drop = FALSE
              ]
              outdmin_s$formula <- prjs[[tstsetup_match_prj]]$outdmin[[s_idx]][[
                "formula"
              ]]
              outdmin_s$x <- outdmin_s$x[
                , colnames(prjs[[tstsetup_match_prj]]$outdmin[[s_idx]]$x),
                drop = FALSE
              ]
              return(outdmin_s)
            }),
            prjs[[tstsetup_match_prj]]$outdmin,
            tolerance = 1e1 * .Machine$double.eps, info = tstsetup
          )
          prj_nms <- names(prjs_cvvs[[tstsetup]])
          expect_identical(prj_nms, names(prjs[[tstsetup_match_prj]]),
                           info = tstsetup)
          prj_el_excl <- !prj_nms %in% c("predictor_terms", "outdmin")
          expect_equal(prjs_cvvs[[tstsetup]][prj_el_excl],
                       prjs[[tstsetup_match_prj]][prj_el_excl],
                       tolerance = .Machine$double.eps, info = tstsetup)
        }
      } else if (length(tstsetup_match_prj) > 1) {
        stop("Unexpected number of matches.")
      }
    } else {
      proj_list_tester(
        prjs_cvvs[[tstsetup]],
        len_expected = length(nterms_crr),
        is_seq = all(diff(nterms_crr) == 1),
        info_str = tstsetup,
        refmod_expected = refmods[[args_prj_cvvs[[tstsetup]]$tstsetup_ref]],
        nprjdraws_expected = args_prj_cvvs[[tstsetup]]$nclusters,
        with_clusters = TRUE,
        prjdraw_weights_expected = prjs_cvvs[[tstsetup]][[1]]$wdraws_prj
      )
    }
  }
})

test_that("`object` not of class `vsel` and missing `predictor_terms` fails", {
  skip_if_not(length(fits) > 0)
  expect_error(
    project(fits[[1]]),
    paste("^Please provide an `object` of class `vsel` or use argument",
          "`predictor_terms`\\.$")
  )
  expect_error(
    project(refmods[[1]]),
    paste("^Please provide an `object` of class `vsel` or use argument",
          "`predictor_terms`\\.$")
  )
})

# nterms ------------------------------------------------------------------

test_that("invalid `nterms` fails", {
  skip_if_not(run_vs)
  tstsetups <- head(grep("\\.glm\\.gauss", names(vss), value = TRUE), 1)
  for (tstsetup in tstsetups) {
    for (nterms_crr in nterms_unavail) {
      expect_error(project(vss[[!!tstsetup]], nterms = !!nterms_crr),
                   paste("Cannot perform the projection with", max(nterms_crr),
                         "variables"))
    }
    for (nterms_crr in list(neg = -1, char = "a", dafr = dat)) {
      expect_error(project(vss[[!!tstsetup]], nterms = !!nterms_crr),
                   "must contain non-negative values")
    }
  }
})

# seed --------------------------------------------------------------------

test_that("non-clustered projection does not require a seed", {
  skip_if_not(run_prj)
  # This test is important to ensure that we don't have to set a seed where we
  # don't expect it to be necessary.
  tstsetups <- grep("\\.noclust$|\\.default_ndr_ncl$", names(prjs),
                    value = TRUE)
  for (tstsetup in tstsetups) {
    args_prj_i <- args_prj[[tstsetup]]
    p_orig <- prjs[[tstsetup]]
    rand_new1 <- runif(1) # Just to advance `.Random.seed[2]`.
    # Use suppressWarnings() because test_that() somehow redirects stderr() and
    # so throws warnings that projpred wants to capture internally:
    p_new <- suppressWarnings(do.call(project, c(
      list(object = refmods[[args_prj_i$tstsetup_ref]]),
      excl_nonargs(args_prj_i, nms_excl_add = "seed")
    )))
    if (args_prj_i$mod_nm %in% c("glmm", "gamm") &&
        any(grepl("\\|", args_prj_i$predictor_terms))) {
      if (getOption("projpred.mlvl_pred_new", FALSE)) {
        # In this case, the multilevel submodel fitters (fit_glmer_callback(),
        # fit_gamm_callback(), fit_cumul_mlvl(), fit_categ_mlvl()) should still
        # be deterministic, but the prediction from the fitted submodels is not
        # (because of the group-level effects drawn randomly by repair_re() (for
        # all group levels; here, only the existing ones are relevant)). Thus,
        # we cannot test the whole project() output, but need to restrict
        # ourselves to the output of as.matrix.projection().
        if (args_prj_i$mod_nm == "gamm") {
          # Skipping GAMMs because of issue #131.
          # TODO (GAMMs): Fix this.
          next
        }
        prjmat_orig <- as.matrix(p_orig)
        prjmat_new <- as.matrix(p_new)
        if (args_prj_i$fam_nm == "gauss" || args_prj_i$prj_nm == "latent") {
          # The projected dispersion parameter is affected by the group-level
          # effects drawn randomly by repair_re() (for all group levels):
          prjmat_new[, "sigma"] <- prjmat_orig[, "sigma"]
        }
        expect_equal(prjmat_new, prjmat_orig, info = tstsetup,
                     tolerance = .Machine$double.eps)
        # To facilitate the `if` conditions here:
        p_new <- p_orig
      } else if (args_prj_i$prj_nm == "augdat" &&
                 args_prj_i$fam_nm == "cumul" && args_prj_i$mod_nm == "glmm") {
        for (idx_s in seq_along(p_new$outdmin)) {
          if (!is.null(p_new$outdmin[[idx_s]][["L"]])) {
            # We could also use `sparseMatrix` instead of `Matrix`:
            expect_equal(as(p_new$outdmin[[idx_s]][["L"]], "Matrix"),
                         as(p_orig$outdmin[[idx_s]][["L"]], "Matrix"),
                         info = tstsetup)
            p_new$outdmin[[idx_s]][["L"]] <- p_orig$outdmin[[idx_s]][["L"]]
          }
        }
      }
    }
    expect_equal(p_new, p_orig, info = tstsetup)
  }
})

test_that("`seed` works (and restores the RNG state afterwards)", {
  skip_if_not(run_prj)
  tstsetups <- grep("\\.glm\\.gauss.*\\.prd_trms_x\\.clust$", names(prjs),
                    value = TRUE)
  for (tstsetup in tstsetups) {
    args_prj_i <- args_prj[[tstsetup]]
    p_orig <- prjs[[tstsetup]]
    rand_orig <- runif(1) # Just to advance `.Random.seed[2]`.
    .Random.seed_new1 <- .Random.seed
    p_new <- do.call(project, c(
      list(object = refmods[[args_prj_i$tstsetup_ref]],
           seed = args_prj_i$seed + 10L),
      excl_nonargs(args_prj_i, nms_excl_add = "seed")
    ))
    .Random.seed_new2 <- .Random.seed
    rand_new <- runif(1) # Just to advance `.Random.seed[2]`.
    .Random.seed_repr1 <- .Random.seed
    p_repr <- do.call(project, c(
      list(object = refmods[[args_prj_i$tstsetup_ref]]),
      excl_nonargs(args_prj_i)
    ))
    .Random.seed_repr2 <- .Random.seed
    # Expected equality:
    expect_equal(p_repr, p_orig, info = tstsetup)
    expect_equal(.Random.seed_new2, .Random.seed_new1, info = tstsetup)
    expect_equal(.Random.seed_repr2, .Random.seed_repr1, info = tstsetup)
    # Expected inequality:
    expect_false(isTRUE(all.equal(p_new, p_orig)), info = tstsetup)
    expect_false(isTRUE(all.equal(rand_new, rand_orig)), info = tstsetup)
    expect_false(isTRUE(all.equal(.Random.seed_repr2, .Random.seed_new2)),
                 info = tstsetup)
  }
})

# regul -------------------------------------------------------------------

test_that("for GLMs, `regul` has an expected effect", {
  skip_if_not(run_prj)
  regul_tst <- c(regul_default, 1e-1, 1e2)
  stopifnot(regul_tst[1] == regul_default)
  stopifnot(all(diff(regul_tst) > 0))
  tstsetups <- grep("\\.glm\\..*\\.clust$", names(prjs), value = TRUE)
  tstsetups <- grep(fam_nms_aug_regex, tstsetups, value = TRUE, invert = TRUE)
  for (tstsetup in tstsetups) {
    args_prj_i <- args_prj[[tstsetup]]
    ndr_ncl <- ndr_ncl_dtls(args_prj_i)

    # Calculate the objects for which to run checks:
    ssq_regul_alpha <- rep(NA, length(regul_tst))
    ssq_regul_beta <- rep(NA, length(regul_tst))
    for (j in seq_along(regul_tst)) {
      # Run project() if necessary:
      if (regul_tst[j] == regul_default) {
        prj_regul <- prjs[[tstsetup]]
      } else {
        prj_regul <- do.call(project, c(
          list(object = refmods[[args_prj_i$tstsetup_ref]],
               regul = regul_tst[j]),
          excl_nonargs(args_prj_i)
        ))
        projection_tester(
          prj_regul,
          refmod_expected = refmods[[args_prj_i$tstsetup_ref]],
          prd_trms_expected = args_prj_i$predictor_terms,
          nprjdraws_expected = ndr_ncl$nprjdraws,
          with_clusters = ndr_ncl$clust_used,
          info_str = paste(tstsetup, j, sep = "__")
        )
      }

      # Run as.matrix.projection():
      prjmat <- as.matrix(prj_regul, nm_scheme = "brms",
                          allow_nonconst_wdraws_prj = ndr_ncl$clust_used_gt1)

      # Reduce to only those columns which are necessary here:
      prjmat <- prjmat[, grep("^b_", colnames(prjmat)), drop = FALSE]

      # Posterior means:
      prjmat_mean <- colMeans(prjmat)

      # Calculate the Euclidean norm for intercept and coefficients:
      ssq_regul_alpha[j] <- prjmat_mean["b_Intercept"]^2
      coef_colnms <- setdiff(names(prjmat_mean), "b_Intercept")
      if (length(coef_colnms) > 0) {
        ssq_regul_beta[j] <- sum(prjmat_mean[coef_colnms]^2)
      }
    }

    # Checks:
    if (length(args_prj_i$predictor_terms) == 0) {
      # For an intercept-only model:
      expect_length(unique(ssq_regul_alpha), 1)
      stopifnot(all(is.na(ssq_regul_beta)))
    } else {
      # All other (i.e., not intercept-only) models (note: as discussed at issue
      # #169, the intercept is not tested here to stay the same):
      for (j in seq_along(ssq_regul_beta)[-1]) {
        expect_lt(ssq_regul_beta[!!j], ssq_regul_beta[j - 1])
      }
    }
  }
})