Skip to contents

sdmTMB_simulate() uses TMB to simulate new data given specified parameter values. simulate.sdmTMB(), on the other hand, takes an existing model fit and simulates new observations and optionally new random effects.

Usage

sdmTMB_simulate(
  formula,
  data,
  mesh,
  family = gaussian(link = "identity"),
  time = NULL,
  B = NULL,
  range = NULL,
  rho = NULL,
  sigma_O = NULL,
  sigma_E = NULL,
  sigma_Z = NULL,
  phi = NULL,
  tweedie_p = NULL,
  df = NULL,
  threshold_coefs = NULL,
  fixed_re = list(omega_s = NULL, epsilon_st = NULL, zeta_s = NULL),
  previous_fit = NULL,
  seed = sample.int(1e+06, 1),
  ...
)

Arguments

formula

A one-sided formula describing the fixed-effect structure. Random intercepts are not (yet) supported. Fixed effects should match the corresponding B argument vector of coefficient values.

data

A data frame containing the predictors described in formula and the time column if time is specified.

mesh

Output from make_mesh().

family

Family as in sdmTMB(). Delta families are not supported. Instead, simulate the two component models separately and combine.

time

The time column name.

B

A vector of beta values (fixed-effect coefficient values).

range

Parameter that controls the decay of spatial correlation. If a vector of length 2, share_range will be set to FALSE and the spatial and spatiotemporal ranges will be unique.

rho

Spatiotemporal correlation between years; should be between -1 and 1.

sigma_O

SD of spatial process (Omega).

sigma_E

SD of spatiotemporal process (Epsilon).

sigma_Z

SD of spatially varying coefficient field (Zeta).

phi

Observation error scale parameter (e.g., SD in Gaussian).

tweedie_p

Tweedie p (power) parameter; between 1 and 2.

df

Student-t degrees of freedom.

threshold_coefs

An optional vector of threshold coefficient values if the formula includes breakpt() or logistic(). If breakpt(), these are slope and cut values. If logistic(), these are the threshold at which the function is 50% of the maximum, the threshold at which the function is 95% of the maximum, and the maximum. See the model description vignette for details.

fixed_re

A list of optional random effects to fix at specified (e.g., previously estimated) values. Values of NULL will result in the random effects being simulated.

previous_fit

(Deprecated; please use simulate.sdmTMB()). An optional previous sdmTMB() fit to pull parameter values. Will be over-ruled by any non-NULL specified parameter arguments.

seed

Seed number.

...

Any other arguments to pass to sdmTMB().

Value

A data frame where:

  • The 1st column is the time variable (if present).

  • The 2nd and 3rd columns are the spatial coordinates.

  • omega_s represents the simulated spatial random effects (only if present).

  • zeta_s represents the simulated spatial varying covariate field (only if present).

  • epsilon_st represents the simulated spatiotemporal random effects (only if present).

  • eta is the true value in link space

  • mu is the true value in inverse link space.

  • observed represents the simulated process with observation error.

  • The remaining columns are the fixed-effect model matrix.

Examples

if (inla_installed()) {
  set.seed(123)

  # make fake predictor(s) (a1) and sampling locations:
  predictor_dat <- data.frame(
    X = runif(300), Y = runif(300),
    a1 = rnorm(300), year = rep(1:6, each = 50)
  )
  mesh <- make_mesh(predictor_dat, xy_cols = c("X", "Y"), cutoff = 0.1)

  sim_dat <- sdmTMB_simulate(
    formula = ~ 1 + a1,
    data = predictor_dat,
    time = "year",
    mesh = mesh,
    family = gaussian(),
    range = 0.5,
    sigma_E = 0.1,
    phi = 0.1,
    sigma_O = 0.2,
    seed = 42,
    B = c(0.2, -0.4) # B0 = intercept, B1 = a1 slope
  )
  head(sim_dat)

  if (require("ggplot2", quietly = TRUE)) {
    ggplot(sim_dat, aes(X, Y, colour = observed)) +
      geom_point() +
      facet_wrap(~year) +
      scale_color_gradient2()
  }

  # fit to the simulated data:
  fit <- sdmTMB(observed ~ a1, data = sim_dat, mesh = mesh, time = "year")
  fit
}
#> Spatiotemporal model fit by ML ['sdmTMB']
#> Formula: observed ~ a1
#> Mesh: mesh
#> Time column: year
#> Data: sim_dat
#> Family: gaussian(link = 'identity')
#>  
#>             coef.est coef.se
#> (Intercept)     0.23    0.09
#> a1             -0.39    0.01
#> 
#> Dispersion parameter: 0.09
#> Matern range: 0.40
#> Spatial SD: 0.21
#> Spatiotemporal SD: 0.11
#> ML criterion at convergence: -162.527
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.