diff --git a/R/posterior-samples.R b/R/posterior-samples.R
index c1698227..d622880b 100644
--- a/R/posterior-samples.R
+++ b/R/posterior-samples.R
@@ -492,14 +492,49 @@
   )
 }
 
+#' @param n_vals numeric; how many locations to evaluate the smooth at if
+#'   `data` not supplied
+#' @param term character; select which smooth's posterior to draw from.
+#'   The default (`NULL`) means the posteriors of all smooths in `model`
+#'   wil be sampled from. If supplied, a character vector of requested terms.
+#' @param rng_per_smooth logical; if TRUE, the behaviour of gratia version
+#' 0.8.1 or earlier is used, whereby a separate call the the random number
+#' generator (RNG) is performed for each smooth. If FALSE, a single call to the
+#' RNG is performed for all model parameters
+#'
+#' @section Warning:
+#' The set of variables returned and their order in the tibble is subject to
+#' change in future versions. Don't rely on position.
+#'
 #' @export
-`smooth_samples.scam` <- function(
+#'
+#' @rdname smooth_samples
+#'
+#' @inheritParams posterior_samples
+#'
+#' @importFrom mvnfast rmvn
+#' @importFrom dplyr bind_rows starts_with
+#' @importFrom tibble as_tibble add_column
+#' @importFrom tidyr pivot_longer
+#' @importFrom mgcv PredictMat
+#' @importFrom tidyselect any_of
+`smooth_samples.gam` <- function(
     model, term = NULL, n = 1, data = newdata,
-    method = c("gaussian", "mh", "inla", "user"), seed = NULL,
-    freq = FALSE, unconditional = FALSE, n_cores = 1L, n_vals = 200,
-    burnin = 1000, thin = 1, t_df = 40, rw_scale = 0.25, rng_per_smooth = FALSE,
-    draws = NULL, mvn_method = c("mvnfast", "mgcv"), ...,
-    newdata = NULL, ncores = NULL) {
+    method = c("gaussian", "mh", "inla", "user"),
+    seed = NULL,
+    freq = FALSE,
+    unconditional = FALSE,
+    n_cores = 1L,
+    n_vals = 200,
+    burnin = 1000,
+    thin = 1,
+    t_df = 40,
+    rw_scale = 0.25,
+    rng_per_smooth = FALSE,
+    draws = NULL,
+    ...,
+    newdata = NULL,
+    ncores = NULL) {
   # smooth_samples begins
   if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) {
     runif(1)
@@ -512,12 +547,14 @@
     RNGstate <- structure(seed, kind = as.list(RNGkind()))
     on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
   }
+
   S <- smooths(model) # vector of smooth labels - "s(x)"
   take <- seq_along(S) # in default case 1,2,3,..,n smooths
   if (!is.null(term)) {
     take <- which_smooths(model, term)
     S <- S[take]
   }
+
   # At least for now, don't work for random effect smooths
   # do we have ranef smooths?
   re_sms <- vapply(get_smooths_by_id(model, take),
@@ -532,24 +569,29 @@
     S <- S[!re_sms]
     take <- take[!re_sms]
   }
+
   # do we have any remaining terms?
   if (length(S) < 1L) {
     stop("No smooths left that can be sampled from.")
   }
+
   if (!is.null(newdata)) {
     newdata_deprecated()
   }
+
   if (!is.null(ncores)) {
     message("Argument `ncores` is deprecated. Use `n_cores` instead.")
     n_cores <- ncores
   }
+
   need_data <- FALSE
   if (is.null(data)) {
     need_data <- TRUE
   }
+
   # what posterior sampling are we using
   method <- match.arg(method)
-  mvn_method <- match.arg(mvn_method)
+
   # get posterior draws - call this once for all parameters
   if (isFALSE(rng_per_smooth)) {
     betas <- post_draws(
@@ -557,9 +599,10 @@
       n_cores = n_cores, model = model,
       burnin = burnin, thin = thin, t_df = t_df, rw_scale = rw_scale,
       index = NULL, frequentist = freq, unconditional = unconditional,
-      draws = draws, mvn_method = mvn_method, seed = seed
+      draws = draws, seed = seed,
     )
   }
+
   sims <- data_names <- vector("list", length = length(S))
   for (i in seq_along(S)) {
     if (need_data) {
@@ -585,13 +628,7 @@
       )
       Xp %*% t(betas)
     } else {
-      # In a scam model, the intercept seems to be parametrised into the
-      # smooth as PredictMat returns a constant column
-      if (isTRUE(length(idx) < ncol(Xp))) {
-        Xp[, -1, drop = FALSE] %*% t(betas[, idx, drop = FALSE])
-      } else {
-        Xp %*% t(betas[, idx, drop = FALSE])
-      }
+      Xp %*% t(betas[, idx, drop = FALSE])
     }
     colnames(simu) <- paste0("..V", seq_len(NCOL(simu)))
     simu <- as_tibble(simu)
@@ -618,7 +655,7 @@
       names_to = ".draw", values_to = ".value",
       names_transform = \(x) as.integer(sub("\\.\\.V", "", x))
     )
-    simu <- relocate(simu, names(data), .after = last_col()) |>
+    simu <- relocate(simu, any_of(names(data)), .after = last_col()) |>
       relocate(".by", .after = 3L)
     ## nest all columns with varying data
     simu <- nest(simu, data = all_of(c(
@@ -636,48 +673,20 @@
   sims
 }
 
-#' @param n_vals numeric; how many locations to evaluate the smooth at if
-#'   `data` not supplied
-#' @param term character; select which smooth's posterior to draw from.
-#'   The default (`NULL`) means the posteriors of all smooths in `model`
-#'   wil be sampled from. If supplied, a character vector of requested terms.
-#' @param rng_per_smooth logical; if TRUE, the behaviour of gratia version
-#' 0.8.1 or earlier is used, whereby a separate call the the random number
-#' generator (RNG) is performed for each smooth. If FALSE, a single call to the
-#' RNG is performed for all model parameters
-#'
-#' @section Warning:
-#' The set of variables returned and their order in the tibble is subject to
-#' change in future versions. Don't rely on position.
-#'
 #' @export
-#'
-#' @rdname smooth_samples
-#'
-#' @inheritParams posterior_samples
-#'
 #' @importFrom mvnfast rmvn
 #' @importFrom dplyr bind_rows starts_with
 #' @importFrom tibble as_tibble add_column
 #' @importFrom tidyr pivot_longer
 #' @importFrom mgcv PredictMat
-`smooth_samples.gam` <- function(
+#' @importFrom tidyselect any_of
+`smooth_samples.scam` <- function(
     model, term = NULL, n = 1, data = newdata,
-    method = c("gaussian", "mh", "inla", "user"),
-    seed = NULL,
-    freq = FALSE,
-    unconditional = FALSE,
-    n_cores = 1L,
-    n_vals = 200,
-    burnin = 1000,
-    thin = 1,
-    t_df = 40,
-    rw_scale = 0.25,
-    rng_per_smooth = FALSE,
-    draws = NULL,
-    ...,
-    newdata = NULL,
-    ncores = NULL) {
+    method = c("gaussian", "mh", "inla", "user"), seed = NULL,
+    freq = FALSE, unconditional = FALSE, n_cores = 1L, n_vals = 200,
+    burnin = 1000, thin = 1, t_df = 40, rw_scale = 0.25, rng_per_smooth = FALSE,
+    draws = NULL, mvn_method = c("mvnfast", "mgcv"), ...,
+    newdata = NULL, ncores = NULL) {
   # smooth_samples begins
   if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) {
     runif(1)
@@ -690,14 +699,12 @@
     RNGstate <- structure(seed, kind = as.list(RNGkind()))
     on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
   }
-
   S <- smooths(model) # vector of smooth labels - "s(x)"
   take <- seq_along(S) # in default case 1,2,3,..,n smooths
   if (!is.null(term)) {
     take <- which_smooths(model, term)
     S <- S[take]
   }
-
   # At least for now, don't work for random effect smooths
   # do we have ranef smooths?
   re_sms <- vapply(get_smooths_by_id(model, take),
@@ -712,29 +719,24 @@
     S <- S[!re_sms]
     take <- take[!re_sms]
   }
-
   # do we have any remaining terms?
   if (length(S) < 1L) {
     stop("No smooths left that can be sampled from.")
   }
-
   if (!is.null(newdata)) {
     newdata_deprecated()
   }
-
   if (!is.null(ncores)) {
     message("Argument `ncores` is deprecated. Use `n_cores` instead.")
     n_cores <- ncores
   }
-
   need_data <- FALSE
   if (is.null(data)) {
     need_data <- TRUE
   }
-
   # what posterior sampling are we using
   method <- match.arg(method)
-
+  mvn_method <- match.arg(mvn_method)
   # get posterior draws - call this once for all parameters
   if (isFALSE(rng_per_smooth)) {
     betas <- post_draws(
@@ -742,10 +744,9 @@
       n_cores = n_cores, model = model,
       burnin = burnin, thin = thin, t_df = t_df, rw_scale = rw_scale,
       index = NULL, frequentist = freq, unconditional = unconditional,
-      draws = draws, seed = seed,
+      draws = draws, mvn_method = mvn_method, seed = seed
     )
   }
-
   sims <- data_names <- vector("list", length = length(S))
   for (i in seq_along(S)) {
     if (need_data) {
@@ -771,20 +772,17 @@
       )
       Xp %*% t(betas)
     } else {
-      Xp %*% t(betas[, idx, drop = FALSE])
+      # In a scam model, the intercept seems to be parametrised into the
+      # smooth as PredictMat returns a constant column
+      if (isTRUE(length(idx) < ncol(Xp))) {
+        Xp[, -1, drop = FALSE] %*% t(betas[, idx, drop = FALSE])
+      } else {
+        Xp %*% t(betas[, idx, drop = FALSE])
+      }
     }
     colnames(simu) <- paste0("..V", seq_len(NCOL(simu)))
     simu <- as_tibble(simu)
     nr_simu <- nrow(simu)
-    # is_fac_by <- is_factor_by_smooth(sm)
-    # if (is_fac_by)) {
-    #     simu <- add_factor_by_data(simu, n = n_vals,
-    #                                by_name = by_variable(sm),
-    #                                by_data = data, before = 1L)
-    # } else {
-    #     simu <- add_by_var_column(simu, n = nr_simu,
-    #         by_var = by_variable(sm))
-    # }
     simu <- add_by_data(simu,
       by_name = by_variable(sm), by_data = data,
       before = 1L
@@ -807,7 +805,7 @@
       names_to = ".draw", values_to = ".value",
       names_transform = \(x) as.integer(sub("\\.\\.V", "", x))
     )
-    simu <- relocate(simu, names(data), .after = last_col()) |>
+    simu <- relocate(simu, any_of(names(data)), .after = last_col()) |>
       relocate(".by", .after = 3L)
     ## nest all columns with varying data
     simu <- nest(simu, data = all_of(c(