[avehtari]

This PR was assisted by Claude.

Summary

predict.subfit() builds the model matrix for new data with model.matrix(delete.response(terms(object$formula)), data = newdata, xlev = object$xlvls). These terms do not have the predvars attribute of the training model frame, so model.matrix() computes terms with a data-dependent basis, such as poly() or splines::ns(), again from the new data. The coefficients of the submodel then apply to other basis columns.

Consequences:

1. The predictions of projected submodels for new data are wrong, without a warning. This affects proj_linpred() and proj_predict() with newdata, varsel() with d_test, and K-fold CV, where the submodels predict the left-out fold. In the example below, the predictions for 10 training rows given as newdata differ from the training-data predictions of the same rows by up to 5.8 on the linear-predictor scale (the outcome has an SD of about 1.4).
2. cv_varsel(validate_search = TRUE) with LOO CV fails: for one left-out observation, poly(x3, 2) stops with 'degree' must be less than number of unique points.

Reproduction on master (custom reference model, no Stan fit):

library(projpred)
set.seed(6283)
n <- 80
dat <- data.frame(x1 = rnorm(n), x3 = rnorm(n))
dat$y <- 0.8 * dat$x1 + 0.5 * dat$x3^2 + rnorm(n)
fml <- y ~ x1 + poly(x3, 2)
mf <- model.frame(fml, dat)
mm <- model.matrix(attr(mf, "terms"), mf)
beta <- drop(solve(crossprod(mm), crossprod(mm, dat$y))) +
  matrix(rnorm(ncol(mm) * 200, sd = 0.05), ncol(mm))
fit <- list(beta = beta, terms = attr(mf, "terms"), data = dat)
ref_predfun <- function(fit, newdata = NULL) {
  nd <- if (is.null(newdata)) fit$data else newdata
  tt <- delete.response(fit$terms)
  model.matrix(tt, model.frame(tt, nd)) %*% fit$beta
}
refm <- init_refmodel(fit, data = dat, formula = fml, family = gaussian(),
                      ref_predfun = ref_predfun, dis = rep(1, 200))
prj <- project(refm, predictor_terms = c("x1", "poly(x3, 2)"), ndraws = 20)
max(abs(proj_linpred(prj)$pred[, 1:10] -
        proj_linpred(prj, newdata = dat[1:10, ])$pred))
# 5.8 instead of 0

cv_varsel(refm, validate_search = TRUE, nloo = 5, nclusters = 3, nterms_max = 2)
# Error in poly(x3, 2): 'degree' must be less than number of unique points

Changes

• fit_glm_ridge_callback() stores the terms of its training model frame (with predvars) as the new element terms of the subfit object. The L1 search stores the terms of its model frame in the same element.
• The new subfit_terms_newdata() builds the terms of object$formula without the response and copies the predvars entries of its variables from object$terms. predict.subfit() uses these terms for new data. If object$terms is missing or does not contain all variables (for example, a subfit object from an older projpred version), the function returns the plain terms, as before.
• The test helper outdmin_tester_trad() expects the new element.
• NEWS.md entry.

The other submodel fitters (lme4, mgcv, gamm4, MASS, nnet, ordinal, mclogit) use their own predict() methods, which keep predvars; I did not check them. The reference model is not affected, because its predictions come from its own predict() method.

Results

For the example above, the predictions for new data now equal the training predictions of the same rows (difference about 1e-16). The predictions for the training data, the projected draws and all results for models without such terms are identical to master.

Tests

The new test_subfit_predvars.R uses a Gaussian reference model with a poly() or splines::ns() term. It checks that the predictions for new data equal the predictions for the same training observations (forward projection with poly() and with ns(), and L1-search submodels with poly()), and that cv_varsel(validate_search = TRUE) runs. All four tests fail on master (three with wrong predictions, one with the error above). The full suite passes with NOT_CRAN=true: 257 tests, 55 810 expectations, 0 failures, the same 4 warnings as master. Tested with this change applied to master 73cf0ee (projpred 2.11.0.9000), R 4.6.1, brms 2.23.0, rstanarm 2.32.2, cmdstanr 0.9.0.9001, CmdStan 2.39.0, on CSC Roihu (Linux, AMD EPYC).