Simulate from a fitted sdmTMB model

simulate.sdmTMB is an S3 method for producing a matrix of simulations from a fitted model. This is similar to lme4::simulate.merMod() and glmmTMB::simulate.glmmTMB(). It can be used with the DHARMa package among other uses.

Usage

# S3 method for sdmTMB
simulate(
  object,
  nsim = 1L,
  seed = sample.int(1e+06, 1L),
  params = c("mle", "mvn"),
  model = c(NA, 1, 2),
  re_form = NULL,
  mcmc_samples = NULL,
  ...
)

Arguments

object

sdmTMB model

nsim

Number of response lists to simulate. Defaults to 1.

seed

Random number seed

params

Whether the parameters used in the simulation should come from the Maximum Likelihood Estimate ("mle") or from new draws from the joint precision matrix assuming they are multivariate normal distributed ("mvn").

model

If a delta/hurdle model, which model to simulate from? NA = combined, 1 = first model, 2 = second mdoel.

re_form

NULL to specify a simulation conditional on fitted random effects (this only simulates observation error). ~0 or NA to simulate new random affects (smoothers, which internally are random effects, will not be simulated as new).

mcmc_samples

An optional matrix of MCMC samples. See extract_mcmc() in the sdmTMBextra package.

...

Extra arguments (not used)

Value

Returns a matrix; number of columns is nsim.

See also

sdmTMB_simulate()

Examples

if (inla_installed()) {

# start with some data simulated from scratch:
set.seed(1)
predictor_dat <- data.frame(X = runif(300), Y = runif(300), a1 = rnorm(300))
mesh <- make_mesh(predictor_dat, xy_cols = c("X", "Y"), cutoff = 0.1)
dat <- sdmTMB_simulate(
  formula = ~ 1 + a1,
  data = predictor_dat,
  mesh = mesh,
  family = poisson(),
  range = 0.5,
  sigma_O = 0.2,
  seed = 42,
  B = c(0.2, -0.4) # B0 = intercept, B1 = a1 slope
)
fit <- sdmTMB(observed ~ 1 + a1, data = dat, family = poisson(), mesh = mesh)

# simulate from the model:
s1 <- simulate(fit, nsim = 300)
dim(s1)

# test whether fitted models are consistent with the observed number of zeros:
sum(s1 == 0)/length(s1)
sum(dat$observed == 0) / length(dat$observed)

# simulate with the parameters drawn from the joint precision matrix:
s2 <- simulate(fit, nsim = 1, params = "MVN")

# simulate with new random fields:
s3 <- simulate(fit, nsim = 1, re_form = ~ 0)

# simulate with new random fields and new parameter draws:
s3 <- simulate(fit, nsim = 1, params = "MVN", re_form = ~ 0)
}