diff --git a/R/simulate-methods.R b/R/simulate-methods.R
index 5f75d90b1..8d033d825 100644
--- a/R/simulate-methods.R
+++ b/R/simulate-methods.R
@@ -1,9 +1,13 @@
 ##' Simulate from the posterior distribution of a GAM
 ##'
 ##' Simulations from the posterior distribution of a fitted GAM model involve
-##'   making random draws from a multivariate normal with mean vector equal to
-##'   the estimated model coefficients and covariance matrix equal to the
-##'   covariance matrix of the coefficients.
+##' computing predicted values for the observation data for which simulated
+##' data are required, then generating random draws from the probability
+##' distribution used when fitting the model.
+##'
+##' For `simulate.gam()` to function, the `family` component of the fitted
+##' model must contain, or be updateable to contain, the required random
+##' number generator. See [mgcv::fix.family.rd()].
 ##'
 ##' @param object a fitted GAM, typically the result of a call to [mgcv::gam]`
 ##'   or [mgcv::gamm()].
@@ -12,15 +16,13 @@
 ##' @param newdata data frame; new observations at which the posterior draws
 ##'   from the model should be evaluated. If not supplied, the data used to fit
 ##'   the model will be used for `newdata`, if available in `object`.
-##' @param freq logical; `TRUE` to return the frequentist covariance matrix of
-##'   the parameter estimators, `FALSE` to return the Bayesian posterior
-##'   covariance matrix of the parameters.
-##' @param unconditional logical; if `TRUE` (and `freq == FALSE`) then the
-##'   Bayesian smoothing parameter uncertainty corrected covariance matrix is
-##'   returned, if available.
 ##' @param weights numeric; a vector of prior weights. If `newdata` is null
 ##'   then defaults to `object[["prior.weights"]]`, otherwise a vector of ones.
-##' @param ... arguments passed to methods
+##' @param ... arguments passed to methods. `simulate.gam()` and
+##'   `simulate.scam()` pass `...` on to `predict.gam()`. As such you can pass
+##'   additional arguments such as `terms`, `exclude`, to select which model
+##'   terms are included in the predictions. This may be useful, for example,
+##'   for excluding the effects of random effect terms.
 ##'
 ##' @return (Currently) A matrix with `nsim` columns.
 ##'
@@ -42,7 +44,6 @@
 ##' sims <- simulate(m1, nsim = 5, seed = 42)
 ##' head(sims)
 `simulate.gam` <- function(object, nsim = 1, seed = NULL, newdata = NULL,
-                           freq = FALSE, unconditional = FALSE,
                            weights = NULL, ...) {
     if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) {
         runif(1)
@@ -73,7 +74,7 @@
         }
     }
 
-    mu <- predict(object, newdata = newdata, type = "response")
+    mu <- predict(object, newdata = newdata, type = "response", ...)
     sims <- replicate(nsim, rd_fun(mu = mu, wt = weights, scale = scale))
 
     attr(sims, "seed") <- RNGstate
@@ -84,10 +85,8 @@
 ##'
 ##' @export
 `simulate.gamm` <- function(object, nsim = 1, seed = NULL, newdata = NULL,
-                            freq = FALSE, unconditional = FALSE, weights = NULL,
-                            ...) {
+                            weights = NULL, ...) {
     simulate(object$gam, nsim = nsim, seed = seed, newdata = newdata,
-             freq = freq, unconditional = unconditional,
              weights = weights, ...)
 }
 
@@ -96,7 +95,7 @@
 ##'
 ##' @export
 `simulate.scam` <- function(object, nsim = 1, seed = NULL, newdata = NULL,
-                            freq = FALSE, weights = NULL, ...) {
+                            weights = NULL, ...) {
     if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) {
         runif(1)
     }
@@ -126,7 +125,7 @@
         }
     }
 
-    mu <- predict(object, newdata = newdata, type = "response")
+    mu <- predict(object, newdata = newdata, type = "response", ...)
 
     ## sims <- mu
     sims <- replicate(nsim, rd_fun(mu = mu, wt = weights, scale = scale))
diff --git a/man/simulate.Rd b/man/simulate.Rd
index 42eb34ab1..56e18eaae 100644
--- a/man/simulate.Rd
+++ b/man/simulate.Rd
@@ -6,37 +6,11 @@
 \alias{simulate.scam}
 \title{Simulate from the posterior distribution of a GAM}
 \usage{
-\method{simulate}{gam}(
-  object,
-  nsim = 1,
-  seed = NULL,
-  newdata = NULL,
-  freq = FALSE,
-  unconditional = FALSE,
-  weights = NULL,
-  ...
-)
+\method{simulate}{gam}(object, nsim = 1, seed = NULL, newdata = NULL, weights = NULL, ...)
 
-\method{simulate}{gamm}(
-  object,
-  nsim = 1,
-  seed = NULL,
-  newdata = NULL,
-  freq = FALSE,
-  unconditional = FALSE,
-  weights = NULL,
-  ...
-)
+\method{simulate}{gamm}(object, nsim = 1, seed = NULL, newdata = NULL, weights = NULL, ...)
 
-\method{simulate}{scam}(
-  object,
-  nsim = 1,
-  seed = NULL,
-  newdata = NULL,
-  freq = FALSE,
-  weights = NULL,
-  ...
-)
+\method{simulate}{scam}(object, nsim = 1, seed = NULL, newdata = NULL, weights = NULL, ...)
 }
 \arguments{
 \item{object}{a fitted GAM, typically the result of a call to \link[mgcv:gam]{mgcv::gam}`
@@ -50,27 +24,28 @@ or \code{\link[mgcv:gamm]{mgcv::gamm()}}.}
 from the model should be evaluated. If not supplied, the data used to fit
 the model will be used for \code{newdata}, if available in \code{object}.}
 
-\item{freq}{logical; \code{TRUE} to return the frequentist covariance matrix of
-the parameter estimators, \code{FALSE} to return the Bayesian posterior
-covariance matrix of the parameters.}
-
-\item{unconditional}{logical; if \code{TRUE} (and \code{freq == FALSE}) then the
-Bayesian smoothing parameter uncertainty corrected covariance matrix is
-returned, if available.}
-
 \item{weights}{numeric; a vector of prior weights. If \code{newdata} is null
 then defaults to \code{object[["prior.weights"]]}, otherwise a vector of ones.}
 
-\item{...}{arguments passed to methods}
+\item{...}{arguments passed to methods. \code{simulate.gam()} and
+\code{simulate.scam()} pass \code{...} on to \code{predict.gam()}. As such you can pass
+additional arguments such as \code{terms}, \code{exclude}, to select which model
+terms are included in the predictions. This may be useful, for example,
+for excluding the effects of random effect terms.}
 }
 \value{
 (Currently) A matrix with \code{nsim} columns.
 }
 \description{
 Simulations from the posterior distribution of a fitted GAM model involve
-making random draws from a multivariate normal with mean vector equal to
-the estimated model coefficients and covariance matrix equal to the
-covariance matrix of the coefficients.
+computing predicted values for the observation data for which simulated
+data are required, then generating random draws from the probability
+distribution used when fitting the model.
+}
+\details{
+For \code{simulate.gam()} to function, the \code{family} component of the fitted
+model must contain, or be updateable to contain, the required random
+number generator. See \code{\link[mgcv:fix.family.rd]{mgcv::fix.family.rd()}}.
 }
 \examples{
 load_mgcv()