Merge branch 'issue-#74': close #74

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))