diff --git a/NEWS.md b/NEWS.md
index 94ffe21cd..3b032c1a1 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -7,6 +7,12 @@ If you read this from a place other than <https://mc-stan.org/projpred/news/inde
 
 ## Major changes
 
+* Faster predictions of multilevel submodels (`lme4::lmer()`, `lme4::glmer()`) for data without new group levels, for example the training data: one model matrix and one random-effects design serve all projected draws. The predictions are unchanged. (GitHub: #584)
+
+* Faster projections onto Gaussian multilevel submodels (`lme4::lmer()`): all projected draws share one `lme4::lFormula()` result. The fits are unchanged. If a check at the first use finds that the installed lme4 version gives other fits, projpred calls `lme4::lmer()` for each draw. Together with the faster predictions above, `cv_varsel()` with `validate_search = TRUE` for a Gaussian multilevel model is about 3 times faster. (GitHub: #584)
+
+* Faster projections onto additive submodels without group-level terms (`mgcv::gam()`): all projected draws share one setup of `mgcv::gam()`, and one linear-predictor matrix gives the predictions of all draws. The fits and predictions are unchanged. For example, `varsel()` for an additive model is about 1.5 times faster. (GitHub: #584)
+
 ## Minor changes
 
 ## Bug fixes and maintenance
diff --git a/R/divergence_minimizers.R b/R/divergence_minimizers.R
index 36345eda8..e1ac419f9 100644
--- a/R/divergence_minimizers.R
+++ b/R/divergence_minimizers.R
@@ -70,6 +70,23 @@ divmin <- function(
     projpred_random <- formula_random$random
   }
   formulas <- validate_response_formula(formula)
+  # Gaussian multilevel submodels: one lFormula() for all draws (see
+  # lmer_shared_callback()). The options and the lme4 check run here, because
+  # the draw-wise calls may run in parallel workers:
+  if (identical(sdivmin, fit_glmer_callback) && length(formulas) > 1 &&
+      !getOption("projpred.PQL", FALSE) &&
+      getOption("projpred.gaussian_not_as_generalized", TRUE) &&
+      identical(list(...)$family$family, "gaussian") &&
+      identical(list(...)$family$link, "identity") &&
+      getOption("projpred.lmer_shared", TRUE) && lmer_shared_ok()) {
+    sdivmin <- shared_setup_callback(lmer_shared_callback, formulas, sdivmin)
+  }
+  # Additive submodels without group-level terms: one mgcv setup for all draws
+  # (see gam_shared_callback()):
+  if (identical(sdivmin, fit_gam_callback) && length(formulas) > 1 &&
+      getOption("projpred.gam_shared", TRUE)) {
+    sdivmin <- shared_setup_callback(gam_shared_callback, formulas, sdivmin)
+  }
   if (is.list(projpred_formulas_no_random)) {
     stopifnot(length(projpred_formulas_no_random) == length(formulas))
   } else {
@@ -169,6 +186,324 @@ divmin <- function(
   return(outdmin)
 }
 
+# Draw-wise divergence minimizer for divmin() that builds the shared setup
+# (`make_shared`: lmer_shared_callback() or gam_shared_callback())
+# from the arguments of its first call. The setup then uses the data of the
+# draw-wise fits, also if divmin() changes the data before these fits. A call
+# with other data, or a setup that does not apply, uses `sdivmin_plain`. The
+# setup is built once, or once per task in parallel runs.
+shared_setup_callback <- function(make_shared, formulas, sdivmin_plain) {
+  # Force the arguments: divmin() assigns the result to `sdivmin`, so a
+  # promise for `sdivmin_plain` would return this function (infinite
+  # recursion):
+  force(sdivmin_plain)
+  force(make_shared)
+  force(formulas)
+  sdivmin_shared <- NULL
+  data_shared <- NULL
+  done <- FALSE
+  return(function(formula, ...) {
+    if (!done) {
+      done <<- TRUE
+      dot_args <- list(...)
+      sdivmin_shared <<- make_shared(formulas, dot_args = dot_args)
+      data_shared <<- dot_args[["data"]]
+    }
+    # identical() returns fast if both are the same object (the usual case):
+    if (is.null(sdivmin_shared) ||
+        !identical(list(...)[["data"]], data_shared)) {
+      return(sdivmin_plain(formula = formula, ...))
+    }
+    return(sdivmin_shared(formula = formula, ...))
+  })
+}
+
+# Additive submodels (mgcv::gam()) of all draws of one projection
+# ---------------------------------------------------------------
+#
+# The draw-wise formulas differ only in the response, so all draws can share the
+# setup of mgcv::gam() (model frame, smooth bases and penalties;
+# mgcv::gam(fit = FALSE)). gam_shared_callback() fits each draw with
+# mgcv::gam(G = ), the public interface for a prepared setup, with the response
+# of the draw. It sets the elements of the fit that contain the response name
+# (`formula`, `terms`, `pterms`, `model`, `pred.formula`) and `call` as in a fit
+# by fit_gam_callback(), so the fits are equal (all.equal()) to those of
+# fit_gam_callback().
+
+# Draw-wise divergence minimizer for divmin() that uses one mgcv::gam() setup
+# for all draws, or `NULL` if this does not apply.
+gam_shared_callback <- function(formulas, dot_args) {
+  if (!getOption("projpred.gam_shared", TRUE)) {
+    return(NULL)
+  }
+  data <- dot_args$data
+  resp_nms <- vapply(formulas, function(fml) {
+    if (is.name(fml[[2]])) as.character(fml[[2]]) else NA_character_
+  }, character(1))
+  if (!is.data.frame(data) || anyNA(resp_nms) ||
+      !all(resp_nms %in% names(data)) || anyNA(data[resp_nms])) {
+    return(NULL)
+  }
+  # As in fit_gam_callback():
+  gam_dot_args <- dot_args[intersect(
+    names(dot_args),
+    c(methods::formalArgs(mgcv::gam), methods::formalArgs(mgcv::gam.fit))
+  )]
+  if (any(c("G", "fit") %in% names(gam_dot_args))) {
+    return(NULL)
+  }
+  g_setup <- tryCatch(
+    do.call(mgcv::gam, c(list(formula = formulas[[1]]), gam_dot_args,
+                         list(fit = FALSE))),
+    error = function(e) NULL
+  )
+  if (is.null(g_setup) || length(g_setup$y) != nrow(data)) {
+    return(NULL)
+  }
+  resp_old <- resp_nms[1]
+  return(function(formula, ...) {
+    resp_new <- as.character(formula[[2]])
+    tryCatch({
+      g_s <- g_setup
+      g_s$y <- data[[resp_new]]
+      g_s$mf[[resp_old]] <- data[[resp_new]]
+      fit <- mgcv::gam(G = g_s)
+      if (resp_new != resp_old) {
+        for (el in c("formula", "terms", "pterms", "model", "pred.formula")) {
+          if (!is.null(fit[[el]])) {
+            fit[[el]] <- rename_sym(fit[[el]], resp_old, resp_new)
+          }
+        }
+      }
+      fit$call <- match.call(mgcv::gam, as.call(c(
+        list(mgcv::gam), list(formula = formula), gam_dot_args
+      )))
+      fit
+    }, error = function(e) fit_gam_callback(formula = formula, ...))
+  })
+}
+
+# Linear predictors of `gam` fits (without group-level terms) of the same
+# submodel for `newdata`, from one linear-predictor matrix (predict.gam() with
+# `type = "lpmatrix"`) for all draws, instead of one predict.gam() call per
+# draw. Returns `NULL` for other cases.
+subprd_gam_direct <- function(fits, newdata) {
+  if (is.null(newdata) || length(fits) < 2 ||
+      !all(vapply(fits, function(fit) {
+        inherits(fit, "gam") && !inherits(fit, "gamm4")
+      }, logical(1)))) {
+    return(NULL)
+  }
+  fit_1 <- fits[[1]]
+  terms_noenv <- function(trms) {
+    trms <- delete.response(trms)
+    attr(trms, ".Environment") <- NULL
+    trms
+  }
+  trms_1 <- terms_noenv(fit_1$pterms)
+  smooth_1 <- fit_1$smooth
+  if (!all(vapply(fits, function(fit) {
+    identical(terms_noenv(fit$pterms), trms_1) &&
+      identical(fit$smooth, smooth_1) &&
+      identical(names(stats::coef(fit)), names(stats::coef(fit_1)))
+  }, logical(1)))) {
+    return(NULL)
+  }
+  if (length(attr(fit_1$terms, "offset")) ||
+      length(attr(fit_1$pterms, "offset"))) {
+    return(NULL)
+  }
+  x_lp <- stats::predict(fit_1, newdata = newdata, type = "lpmatrix")
+  prd <- x_lp %*% do.call(cbind, lapply(fits, stats::coef))
+  dimnames(prd) <- list(rownames(x_lp), NULL)
+  return(prd)
+}
+
+# Gaussian multilevel submodels (lme4::lmer()) of all draws of one projection
+# --------------------------------------------------------------------------
+#
+# lme4::lmer() is lme4::lFormula() (model frame, fixed-effects model matrix,
+# random-effects structure), followed by mkLmerDevfun(), optimizeLmer(),
+# checkConv() and mkMerMod(). The draw-wise formulas of one submodel differ only
+# in the response, so all draws can share one lFormula() result.
+# lmer_shared_fit() repeats the other steps of lmer() for one draw: it
+# transcribes the body of lmer() in lme4 2.0-6 to 2.1-1 (R package lme4,
+# GPL-2 | GPL-3). lme4's C++ objects change some vectors of the lFormula()
+# result in place, so each draw gets fresh copies of them. The fits are equal
+# (all.equal()) to those of lmer(). lmer_shared_ok() checks this once per R
+# session on a small fixed problem, and the test suite checks it on CRAN. If
+# the check fails (for example after a change in lme4), divmin() fits each draw
+# with lmer().
+
+# Replaces the name `old` by `new` in all language objects of `x`, including
+# its attributes and names (used for the response of a shared model frame).
+rename_sym <- function(x, old, new) {
+  if (is.name(x)) {
+    if (identical(x, as.name(old))) {
+      return(as.name(new))
+    }
+    return(x)
+  }
+  if (is.call(x) || is.pairlist(x) || (is.list(x) && !is.data.frame(x))) {
+    at <- attributes(x)
+    for (i in seq_along(x)) {
+      if (!is.null(x[[i]])) {
+        x[[i]] <- rename_sym(x[[i]], old, new)
+      }
+    }
+    attributes(x) <- at
+  }
+  at <- attributes(x)
+  if (!is.null(at)) {
+    for (nm in setdiff(names(at), c("names", "dim", "dimnames", "row.names",
+                                    "class"))) {
+      attr(x, nm) <- rename_sym(at[[nm]], old, new)
+    }
+    if (!is.null(at$names)) {
+      names(x)[names(x) == old] <- new
+    }
+    if (!is.null(at$dimnames)) {
+      attr(x, "dimnames") <- lapply(at$dimnames, function(dnms) {
+        if (!is.null(dnms)) {
+          dnms[dnms == old] <- new
+        }
+        dnms
+      })
+    }
+  }
+  return(x)
+}
+
+# Fits the lmer() submodel `formula` (response: a column of `data`) from the
+# shared lFormula() result `lmod0` (of the same submodel for another draw).
+# `lmer_args` are the arguments of the lmer() call that this fit replaces.
+lmer_shared_fit <- function(lmod0, formula, data, control, lmer_args) {
+  resp_old <- as.character(lmod0$formula[[2]])
+  resp_new <- as.character(formula[[2]])
+  lmod <- lmod0
+  # Fresh copies of the vectors that lme4's C++ objects change in place:
+  lmod$reTrms$Lambdat@x <- lmod0$reTrms$Lambdat@x + 0
+  lmod$reTrms$theta <- lmod0$reTrms$theta + 0
+  lmod$reTrms$Zt@x <- lmod0$reTrms$Zt@x + 0
+  lmod$X <- lmod0$X + 0
+  fr <- lmod$fr
+  fr[[resp_old]] <- data[[resp_new]]
+  if (resp_new != resp_old) {
+    fr <- rename_sym(fr, resp_old, resp_new)
+  }
+  lmod$fr <- fr
+  # The following steps are those of lme4::lmer():
+  mcout <- match.call(lme4::lmer, as.call(c(list(lme4::lmer), lmer_args)))
+  mcout$formula <- formula
+  lmod$formula <- NULL
+  devfun <- do.call(lme4::mkLmerDevfun,
+                    c(lmod, list(start = NULL, verbose = 0L, control = control)))
+  rho <- environment(devfun)
+  nobs <- nrow(lmod$fr)
+  npar <- length(rho$lower)
+  calc.derivs <- control$calc.derivs %||%
+    (nobs < control$checkConv$check.conv.nobsmax &&
+       npar < control$checkConv$check.conv.nparmax)
+  opt <- lme4::optimizeLmer(
+    devfun, optimizer = control$optimizer,
+    restart_edge = control$restart_edge, boundary.tol = control$boundary.tol,
+    control = control$optCtrl, verbose = 0L, start = NULL,
+    calc.derivs = calc.derivs,
+    force.calc.derivs = control$calc.derivs %||% FALSE,
+    use.last.params = control$use.last.params
+  )
+  # Same name as in lmer(), so that its warnings are the same:
+  checkConv <- "lme4" %:::% "checkConv"
+  cc <- checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,
+                  lbound = rho$lower, ubound = rho$upper, nobs = nobs,
+                  ndim = npar)
+  return(lme4::mkMerMod(rho, opt, lmod$reTrms, fr = lmod$fr, mc = mcout,
+                        lme4conv = cc))
+}
+
+# Draw-wise divergence minimizer for divmin() that uses lmer_shared_fit() with
+# the lFormula() result of the first draw, or `NULL` if this does not apply.
+lmer_shared_callback <- function(formulas, dot_args) {
+  if (!getOption("projpred.lmer_shared", TRUE) || !lmer_shared_ok()) {
+    return(NULL)
+  }
+  data <- dot_args$data
+  resp_nms <- vapply(formulas, function(fml) {
+    if (is.name(fml[[2]])) as.character(fml[[2]]) else NA_character_
+  }, character(1))
+  if (!is.data.frame(data) || anyNA(resp_nms) ||
+      !all(resp_nms %in% names(data)) || anyNA(data[resp_nms])) {
+    return(NULL)
+  }
+  control <- dot_args$control %||% control_callback(dot_args$family)
+  lmer_dot_args <- dot_args[intersect(names(dot_args),
+                                      methods::formalArgs(lme4::lmer))]
+  lmer_dot_args$control <- NULL
+  lmod0 <- tryCatch(
+    do.call(lme4::lFormula, c(list(formula = formulas[[1]], control = control),
+                              lmer_dot_args)),
+    error = function(e) NULL
+  )
+  if (is.null(lmod0) || nrow(lmod0$fr) != nrow(data)) {
+    return(NULL)
+  }
+  force(lmod0)
+  return(function(formula, ...) {
+    lmer_args <- c(list(formula = formula, control = control), lmer_dot_args)
+    tryCatch(
+      lmer_shared_fit(lmod0, formula = formula, data = data, control = control,
+                      lmer_args = lmer_args),
+      error = function(e) fit_glmer_callback(formula = formula, ...)
+    )
+  })
+}
+
+# Whether lmer_shared_fit() gives the fits of lme4::lmer() with the installed
+# lme4 (checked once per R session on a small fixed problem, without random
+# numbers).
+lmer_shared_state <- new.env(parent = emptyenv())
+lmer_shared_ok <- function() {
+  if (is.null(lmer_shared_state$ok)) {
+    lmer_shared_state$ok <- isTRUE(tryCatch(lmer_shared_check(),
+                                            error = function(e) FALSE))
+  }
+  return(lmer_shared_state$ok)
+}
+lmer_shared_check <- function() {
+  n <- 60L
+  idx <- seq_len(n)
+  dat <- data.frame(x1 = sin(idx), x2 = cos(1.7 * idx),
+                    g = factor(paste0("g", idx %% 6L)),
+                    h = factor(paste0("h", idx %% 4L)))
+  u <- c(0.8, -0.5, 0.3, -0.9, 0.6, -0.2)
+  for (s in 1:3) {
+    dat[[paste0(".y.", s)]] <- dat$x1 + 0.5 * dat$x2 +
+      u[as.integer(dat$g)] * (1 + 0.3 * s) + 0.4 * sin(s * idx) +
+      0.2 * dat$x1 * u[as.integer(dat$g)]
+  }
+  wts <- 1 + 0.5 * cos(idx)
+  control <- lme4::lmerControl()
+  fmls <- lapply(1:3, function(s) {
+    as.formula(paste0(".y.", s, " ~ x1 + x2 + (1 + x1 | g) + (1 | h)"))
+  })
+  lmod0 <- lme4::lFormula(fmls[[1]], data = dat, weights = wts,
+                          control = control)
+  same <- logical(3)
+  for (s in 1:3) {
+    lmer_args <- list(formula = fmls[[s]], control = control, data = dat,
+                      weights = wts)
+    suppressMessages(suppressWarnings({
+      fit_ref <- do.call(lme4::lmer, lmer_args)
+      fit_shared <- lmer_shared_fit(lmod0, formula = fmls[[s]], data = dat,
+                                    control = control, lmer_args = lmer_args)
+    }))
+    same[s] <- isTRUE(all.equal(fit_ref, fit_shared, tolerance = 0)) &&
+      identical(lme4::getME(fit_ref, "theta"),
+                lme4::getME(fit_shared, "theta"))
+  }
+  return(all(same))
+}
+
 # Use projpred's own implementation to fit non-multilevel non-additive
 # submodels:
 fit_glm_ridge_callback <- function(formula, data,
@@ -1159,6 +1494,21 @@ check_conv_s <- function(fit_s) {
 # Prediction functions for submodels --------------------------------------
 
 subprd <- function(fits, newdata) {
+  # lme4 fits of one submodel and `newdata` without new group levels: all draws
+  # at once (see subprd_mer_direct()):
+  prd_direct <- subprd_mer_direct(fits, newdata)
+  if (!is.null(prd_direct)) {
+    return(prd_direct)
+  }
+  # mgcv fits of one submodel: all draws at once (see subprd_gam_direct()):
+  if (!is.null(newdata) && length(fits) > 1 && inherits(fits[[1]], "gam")) {
+    prd_direct <- subprd_gam_direct(
+      fits, cbind(`(Intercept)` = rep(1, NROW(newdata)), newdata)
+    )
+    if (!is.null(prd_direct)) {
+      return(prd_direct)
+    }
+  }
   prd_list <- lapply(fits, function(fit) {
     is_glmmPQL <- inherits(fit, "glmmPQL")
     is_glmm <- inherits(fit, c("lmerMod", "glmerMod"))
@@ -1213,6 +1563,97 @@ subprd <- function(fits, newdata) {
   return(do.call(cbind, prd_list))
 }
 
+# Linear predictors of lme4 fits (`lmerMod` or `glmerMod`) of the same submodel
+# for `newdata` whose group levels all occur in the fits. In this case,
+# repair_re() adds zero and draws no random numbers, and predict.merMod() with
+# `re.form = NULL` computes X beta_s + Z b_s. This function builds the
+# fixed-effects model matrix X and the random-effects design Z for `newdata`
+# once, with the steps of predict.merMod() (R package lme4, GPL-2 | GPL-3). It
+# maps the conditional modes b_s of each draw (in the row order of the fitted
+# Zt) onto the rows of Z by the names of the group levels. These names are
+# unique across the grouping variables (see repair_re.merMod()), and within a
+# grouping variable they have the same order of terms and coefficients in both.
+# The results differ from those of subprd() only by rounding. Returns `NULL`
+# for other cases.
+subprd_mer_direct <- function(fits, newdata) {
+  if (is.null(newdata) || length(fits) == 0 ||
+      getOption("projpred.mlvl_pred_new", FALSE) ||
+      !all(vapply(fits, inherits, logical(1), c("lmerMod", "glmerMod")))) {
+    return(NULL)
+  }
+  fit_1 <- fits[[1]]
+  flist_1 <- lapply(fit_1@flist, levels)
+  cnms_1 <- lme4::getME(fit_1, "cnms")
+  trms_fix <- delete.response(terms(fit_1, fixed.only = TRUE))
+  # The terms are compared without their environment:
+  terms_noenv <- function(trms) {
+    attr(trms, ".Environment") <- NULL
+    trms
+  }
+  trms_fix_noenv <- terms_noenv(trms_fix)
+  same_structure <- vapply(fits, function(fit) {
+    trms_s <- delete.response(terms(fit, fixed.only = TRUE))
+    identical(lapply(fit@flist, levels), flist_1) &&
+      identical(lme4::getME(fit, "cnms"), cnms_1) &&
+      identical(terms_noenv(trms_s), trms_fix_noenv) &&
+      identical(names(lme4::fixef(fit)), names(lme4::fixef(fit_1)))
+  }, logical(1))
+  if (!all(same_structure) || !is.null(attr(trms_fix, "offset"))) {
+    return(NULL)
+  }
+  # Group levels: all have to occur in the fits, and the level sets have to be
+  # disjoint (otherwise repair_re.merMod() handles the case):
+  vnms <- names(flist_1)
+  if (!all(vnms %in% names(newdata))) {
+    return(NULL)
+  }
+  for (vnm in vnms) {
+    if (!all(as.character(unique(newdata[[vnm]])) %in% flist_1[[vnm]])) {
+      return(NULL)
+    }
+  }
+  lvls_all <- unlist(flist_1, use.names = FALSE)
+  if (anyDuplicated(lvls_all)) {
+    return(NULL)
+  }
+  # Fixed effects, as in predict.merMod():
+  x_contr <- attr(lme4::getME(fit_1, "X"), "contrasts")
+  mfr <- model.frame(trms_fix, newdata, na.action = na.pass,
+                     xlev = .getXlevels(trms_fix, fit_1@frame))
+  x_new <- model.matrix(trms_fix, mfr, contrasts.arg = x_contr)
+  fix_nms <- names(lme4::fixef(fit_1))
+  if (!all(fix_nms %in% colnames(x_new)) || anyNA(x_new) ||
+      nrow(x_new) != nrow(newdata)) {
+    return(NULL)
+  }
+  x_new <- x_new[, fix_nms, drop = FALSE]
+  # Random effects:
+  re_new <- lme4::mkNewReTrms(fit_1, newdata = newdata,
+                              re.form = ("lme4" %:::% "reOnly")(formula(fit_1)),
+                              allow.new.levels = TRUE)
+  zt_new <- re_new$Zt
+  occ_key <- function(nms) paste(nms, stats::ave(seq_along(nms), nms,
+                                                 FUN = seq_along), sep = "\r")
+  idx_b <- match(occ_key(rownames(zt_new)),
+                 occ_key(rownames(lme4::getME(fit_1, "Zt"))))
+  if (anyNA(idx_b)) {
+    return(NULL)
+  }
+  b_1 <- as.vector(lme4::getME(fit_1, "b"))
+  if (!isTRUE(all.equal(unname(b_1[idx_b]), unname(as.vector(re_new$b)),
+                        tolerance = 1e-12))) {
+    return(NULL)
+  }
+  betas <- do.call(cbind, lapply(fits, lme4::fixef))
+  bs <- do.call(cbind, lapply(fits, function(fit) {
+    as.vector(lme4::getME(fit, "b"))[idx_b]
+  }))
+  prd <- x_new %*% betas + t(as(t(bs) %*% zt_new, "matrix"))
+  # As cbind() of the named vectors from predict.merMod():
+  dimnames(prd) <- list(rownames(newdata), NULL)
+  return(prd)
+}
+
 subprd_augdat <- function(fits, newdata) {
   if (inherits(fits[[1]], "clmm")) {
     y_nm <- extract_terms_response(formula(fits[[1]]))$response
diff --git a/tests/testthat/test_gam_shared.R b/tests/testthat/test_gam_shared.R
new file mode 100644
index 000000000..dc2282686
--- /dev/null
+++ b/tests/testthat/test_gam_shared.R
@@ -0,0 +1,118 @@
+context("gam_shared")
+
+n_gs <- 120L
+idx_gs <- seq_len(n_gs)
+dat_gs <- data.frame(x1 = sin(idx_gs), x2 = cos(1.3 * idx_gs),
+                     x3 = sin(0.7 * idx_gs + 1),
+                     f = factor(c("a", "b", "c")[idx_gs %% 3L + 1L]))
+n_draws_gs <- 4L
+for (s in seq_len(n_draws_gs)) {
+  dat_gs[[paste0(".y.", s)]] <- sin(1.5 * dat_gs$x3) + 0.5 * dat_gs$x1 +
+    0.3 * cos(s * idx_gs)
+}
+wts_gs <- 1 + 0.4 * sin(2 * idx_gs)
+fml_gs <- as.formula(paste0(
+  "cbind(", paste0(".y.", seq_len(n_draws_gs), collapse = ", "), ") ~ ",
+  "s(x3) + s(x2) + x1 + f"
+))
+divmin_gs <- function(shared) {
+  opts_old <- options(projpred.gam_shared = shared)
+  on.exit(options(opts_old))
+  suppressWarnings(divmin(
+    formula = fml_gs, data = dat_gs, family = extend_family(gaussian()),
+    weights = wts_gs, projpred_var = matrix(NA, n_gs, n_draws_gs),
+    projpred_ws_aug = matrix(NA, n_gs, n_draws_gs), verbose_divmin = FALSE
+  ))
+}
+
+test_that("divmin() gives the same gam() fits with and without sharing", {
+  fits_shared <- divmin_gs(TRUE)
+  fits_drawwise <- divmin_gs(FALSE)
+  expect_length(fits_shared, n_draws_gs)
+  for (s in seq_len(n_draws_gs)) {
+    expect_s3_class(fits_shared[[s]], "gam")
+    expect_equal(fits_shared[[s]], fits_drawwise[[s]], tolerance = 0,
+                 info = s)
+    expect_identical(formula(fits_shared[[s]]),
+                     formula(fits_drawwise[[s]]), info = s)
+  }
+})
+
+test_that("subprd() gives the predictions of predict.gam() for gam fits", {
+  fits <- divmin_gs(FALSE)
+  for (nd in list(dat_gs, dat_gs[c(7, 3, 50), ])) {
+    nd_icpt <- cbind(`(Intercept)` = rep(1, nrow(nd)), nd)
+    prd_expected <- do.call(cbind, lapply(fits, function(fit) {
+      predict(fit, newdata = nd_icpt)
+    }))
+    expect_identical(subprd(fits, nd), prd_expected)
+  }
+  # Fits of different submodels:
+  fit_other <- mgcv::gam(.y.1 ~ s(x3) + x1, data = dat_gs)
+  expect_null(subprd_gam_direct(list(fits[[1]], fit_other), dat_gs))
+})
+
+test_that("gam_shared_callback() returns `NULL` for other cases", {
+  fmls <- validate_response_formula(fml_gs)
+  dot_args <- list(data = dat_gs, family = extend_family(gaussian()),
+                   weights = wts_gs)
+  expect_true(is.function(gam_shared_callback(fmls, dot_args)))
+  opts_old <- options(projpred.gam_shared = FALSE)
+  expect_null(gam_shared_callback(fmls, dot_args))
+  options(opts_old)
+  dat_na <- dat_gs
+  dat_na$.y.3[5] <- NA
+  expect_null(gam_shared_callback(fmls, c(list(data = dat_na),
+                                          dot_args[-1])))
+})
+
+test_that("shared_setup_callback() builds the setup once, at its first call", {
+  n_built <- 0L
+  make_shared <- function(formulas, dot_args) {
+    n_built <<- n_built + 1L
+    n_setup <- nrow(dot_args[["data"]])
+    function(formula, ...) c("shared", n_setup)
+  }
+  make_null <- function(formulas, dot_args) {
+    n_built <<- n_built + 1L
+    NULL
+  }
+  plain <- function(formula, ...) c("plain", nrow(list(...)[["data"]]))
+  fmls <- validate_response_formula(fml_gs)
+  dat_red <- dat_gs[-1, ]
+  cb <- shared_setup_callback(make_shared, fmls, plain)
+  expect_identical(n_built, 0L)
+  expect_identical(cb(fmls[[1]], data = dat_red), c("shared", n_gs - 1L))
+  expect_identical(cb(fmls[[2]], data = dat_red), c("shared", n_gs - 1L))
+  # Other data than at the first call:
+  expect_identical(cb(fmls[[3]], data = dat_gs), c("plain", n_gs))
+  expect_identical(n_built, 1L)
+  n_built <- 0L
+  cb_null <- shared_setup_callback(make_null, fmls, plain)
+  expect_identical(cb_null(fmls[[1]], data = dat_red), c("plain", n_gs - 1L))
+  expect_identical(cb_null(fmls[[2]], data = dat_red), c("plain", n_gs - 1L))
+  expect_identical(n_built, 1L)
+})
+
+test_that(paste(
+  "divmin() gives the same gam() fits with and without sharing if an",
+  "observation has weight 0"
+), {
+  opts_old <- options(projpred.gam_shared = TRUE)
+  on.exit(options(opts_old))
+  wts_0 <- replace(wts_gs, 1, 0)
+  divmin_w0 <- function(shared) {
+    options(projpred.gam_shared = shared)
+    suppressWarnings(divmin(
+      formula = fml_gs, data = dat_gs, family = extend_family(gaussian()),
+      weights = wts_0, projpred_var = matrix(NA, n_gs, n_draws_gs),
+      projpred_ws_aug = matrix(NA, n_gs, n_draws_gs), verbose_divmin = FALSE
+    ))
+  }
+  fits_shared <- divmin_w0(TRUE)
+  fits_drawwise <- divmin_w0(FALSE)
+  for (s in seq_len(n_draws_gs)) {
+    expect_equal(fits_shared[[s]], fits_drawwise[[s]], tolerance = 0,
+                 info = s)
+  }
+})
diff --git a/tests/testthat/test_lmer_shared.R b/tests/testthat/test_lmer_shared.R
new file mode 100644
index 000000000..c714c8f15
--- /dev/null
+++ b/tests/testthat/test_lmer_shared.R
@@ -0,0 +1,85 @@
+context("lmer_shared")
+
+# These tests run on CRAN on purpose: lmer_shared_fit() repeats the steps of
+# lme4::lmer() after a shared lme4::lFormula(), so a change in lme4 has to make
+# them fail (also in the reverse-dependency checks of lme4).
+
+test_that("lmer_shared_fit() gives the fits of lme4::lmer()", {
+  expect_true(lmer_shared_check())
+  lmer_shared_state$ok <- NULL
+  expect_true(lmer_shared_ok())
+})
+
+n_ls <- 80L
+idx_ls <- seq_len(n_ls)
+dat_ls <- data.frame(x1 = sin(idx_ls), x2 = cos(1.3 * idx_ls),
+                     f = factor(c("a", "b", "c")[idx_ls %% 3L + 1L]),
+                     g = factor(paste0("g", idx_ls %% 8L)),
+                     h = factor(paste0("h", idx_ls %% 5L)))
+u_ls <- sin(1:8)
+n_draws_ls <- 4L
+for (s in seq_len(n_draws_ls)) {
+  dat_ls[[paste0(".y.", s)]] <- dat_ls$x1 - 0.5 * dat_ls$x2 +
+    u_ls[as.integer(dat_ls$g)] * (1 + 0.2 * s) + 0.3 * cos(s * idx_ls) +
+    0.3 * dat_ls$x1 * u_ls[as.integer(dat_ls$g)]
+}
+wts_ls <- 1 + 0.4 * sin(2 * idx_ls)
+fml_ls <- as.formula(paste0(
+  "cbind(", paste0(".y.", seq_len(n_draws_ls), collapse = ", "), ") ~ ",
+  "x1 + poly(x2, 2, raw = TRUE) + f + (1 + x1 | g) + (1 | h)"
+))
+divmin_ls <- function(shared) {
+  opts_old <- options(projpred.lmer_shared = shared)
+  on.exit(options(opts_old))
+  suppressWarnings(divmin(
+    formula = fml_ls, data = dat_ls, family = extend_family(gaussian()),
+    weights = wts_ls, projpred_var = matrix(NA, n_ls, n_draws_ls),
+    projpred_ws_aug = matrix(NA, n_ls, n_draws_ls), verbose_divmin = FALSE
+  ))
+}
+
+test_that("divmin() gives the same lmer() fits with and without sharing", {
+  fits_shared <- divmin_ls(TRUE)
+  fits_drawwise <- divmin_ls(FALSE)
+  expect_length(fits_shared, n_draws_ls)
+  for (s in seq_len(n_draws_ls)) {
+    expect_s4_class(fits_shared[[s]], "lmerMod")
+    expect_equal(fits_shared[[s]], fits_drawwise[[s]], tolerance = 0,
+                 info = s)
+    expect_identical(lme4::getME(fits_shared[[s]], "theta"),
+                     lme4::getME(fits_drawwise[[s]], "theta"), info = s)
+    expect_identical(formula(fits_shared[[s]]),
+                     formula(fits_drawwise[[s]]), info = s)
+  }
+})
+
+test_that(paste(
+  "divmin() gives the same lmer() fits with and without sharing if an",
+  "observation has weight 0"
+), {
+  wts_bu <- wts_ls
+  on.exit(wts_ls <<- wts_bu)
+  wts_ls <<- replace(wts_ls, 1, 0)
+  fits_shared <- divmin_ls(TRUE)
+  fits_drawwise <- divmin_ls(FALSE)
+  for (s in seq_len(n_draws_ls)) {
+    expect_equal(fits_shared[[s]], fits_drawwise[[s]], tolerance = 0,
+                 info = s)
+    expect_identical(lme4::getME(fits_shared[[s]], "theta"),
+                     lme4::getME(fits_drawwise[[s]], "theta"), info = s)
+  }
+})
+
+test_that("lmer_shared_callback() returns `NULL` for other cases", {
+  fmls <- validate_response_formula(fml_ls)
+  dot_args <- list(data = dat_ls, family = extend_family(gaussian()),
+                   weights = wts_ls)
+  expect_true(is.function(lmer_shared_callback(fmls, dot_args)))
+  opts_old <- options(projpred.lmer_shared = FALSE)
+  expect_null(lmer_shared_callback(fmls, dot_args))
+  options(opts_old)
+  dat_na <- dat_ls
+  dat_na$.y.2[3] <- NA
+  expect_null(lmer_shared_callback(fmls, c(list(data = dat_na),
+                                           dot_args[-1])))
+})
diff --git a/tests/testthat/test_mer_predict.R b/tests/testthat/test_mer_predict.R
new file mode 100644
index 000000000..ab637a789
--- /dev/null
+++ b/tests/testthat/test_mer_predict.R
@@ -0,0 +1,83 @@
+context("mer_predict")
+
+if (exists(".Random.seed", envir = .GlobalEnv)) {
+  rng_old <- get(".Random.seed", envir = .GlobalEnv)
+}
+set.seed(13579)
+n_mp <- 150
+dat_mp <- data.frame(
+  x1 = rnorm(n_mp), x2 = rnorm(n_mp),
+  f = factor(sample(c("a", "b", "c"), n_mp, replace = TRUE)),
+  g = factor(paste0("g", sample(8, n_mp, replace = TRUE))),
+  h = factor(paste0("h", sample(4, n_mp, replace = TRUE)))
+)
+eta_mp <- with(dat_mp, 0.5 * x1 - 0.3 * x2 + rnorm(8)[as.integer(g)] +
+                 0.5 * rnorm(4)[as.integer(h)])
+# Draw-wise fits of one submodel (as from a projection):
+fits_mp <- function(fml, family = NULL, n_draws = 4) {
+  lapply(seq_len(n_draws), function(s) {
+    dat_s <- dat_mp
+    dat_s$y <- eta_mp + rnorm(n_mp, 0, 0.3)
+    if (is.null(family)) {
+      suppressMessages(lme4::lmer(fml, data = dat_s))
+    } else {
+      dat_s$y <- plogis(dat_s$y)
+      suppressWarnings(lme4::glmer(fml, data = dat_s, family = family))
+    }
+  })
+}
+# The per-draw predictions of subprd() on master:
+prd_per_draw <- function(fits, newdata) {
+  do.call(cbind, lapply(fits, function(fit) {
+    predict(fit, newdata = newdata, re.form = NULL, allow.new.levels = TRUE) +
+      repair_re(fit, newdata = newdata)
+  }))
+}
+
+test_that("subprd_mer_direct() gives the predictions of the draw-wise path", {
+  cases <- list(
+    lmer_slope = fits_mp(y ~ x1 + x2 + f + (1 + x1 | g) + (1 | h)),
+    lmer_icpt = fits_mp(y ~ x1 + (1 | g)),
+    lmer_nofix = fits_mp(y ~ 1 + (1 | g) + (1 | h)),
+    glmer = fits_mp(y ~ x1 + (1 | g) + (1 | h), family = binomial()),
+    # Special formula (functions in the fixed-effects terms):
+    lmer_spcl = fits_mp(y ~ poly(x1, 2, raw = TRUE) +
+                          sqrt(abs(x2)^2) * I(!as.logical(x2 > 0)) + f +
+                          (x1 | g))
+  )
+  newdat <- dat_mp[sample(n_mp, 30), ]
+  for (nm in names(cases)) {
+    for (nd in list(dat_mp, newdat)) {
+      set.seed(1)
+      prd_expected <- prd_per_draw(cases[[nm]], nd)
+      seed_expected <- .Random.seed
+      set.seed(1)
+      prd <- subprd_mer_direct(cases[[nm]], nd)
+      expect_false(is.null(prd), info = nm)
+      expect_equal(prd, prd_expected, tolerance = 1e-12, info = nm)
+      # No random numbers are drawn:
+      expect_identical(.Random.seed, seed_expected, info = nm)
+      expect_equal(subprd(cases[[nm]], nd), prd_expected, tolerance = 1e-12,
+                   info = nm)
+    }
+  }
+})
+
+test_that("subprd_mer_direct() returns `NULL` for other cases", {
+  fits <- fits_mp(y ~ x1 + (1 | g))
+  # New group level:
+  nd_new <- dat_mp[1:10, ]
+  nd_new$g <- factor(c("gnew", as.character(nd_new$g[-1])))
+  expect_null(subprd_mer_direct(fits, nd_new))
+  # All group levels are treated as new:
+  opts_old <- options(projpred.mlvl_pred_new = TRUE)
+  expect_null(subprd_mer_direct(fits, dat_mp))
+  options(opts_old)
+  # Other fits, no `newdata`, fits of different submodels:
+  expect_null(subprd_mer_direct(list(lm(x1 ~ x2, data = dat_mp)), dat_mp))
+  expect_null(subprd_mer_direct(fits, NULL))
+  fits_2 <- c(fits[1], fits_mp(y ~ x2 + (1 | g), n_draws = 1))
+  expect_null(subprd_mer_direct(fits_2, dat_mp))
+})
+
+if (exists("rng_old")) assign(".Random.seed", rng_old, envir = .GlobalEnv)