Example of simulating with sdmTMB

sdmTMB-RCA-sim.R

library(sdmTMB)
library(dplyr)
library(ggplot2)

mesh <- make_mesh(pcod_2011, c("X", "Y"), cutoff = 10)
fit <- sdmTMB(
  density ~ poly(depth, 2),
  data = pcod_2011, mesh = mesh,
  family = tweedie(link = "log")
)

b1 <- tidy(fit)
b2 <- tidy(fit, "ran_pars")
b <- bind_rows(b1, b2)

pars <- get_pars(fit)

# this would be new locations and depths; for now, just use what we have with
# some fake RCAs
dat <- select(pcod_2011, X, Y, depth)
dat <- bind_rows(mutate(dat, later_survey = 0L), mutate(dat, later_survey = 1L))
random_rcas <- sample(seq_len(nrow(pcod_2011)), 100)
dat$rca <- 0L
dat[random_rcas, "rca"] <- 1L
head(dat)

# figure out model matrix for below:
model.matrix(~ later_survey * rca + poly(depth, 2), data = dat) |> head()

sim_mesh <- make_mesh(dat, c("X", "Y"), mesh = mesh$mesh)

s <- sdmTMB_simulate(
  ~ later_survey * rca + poly(depth, 2),
  data = dat,
  mesh = sim_mesh,
  family = tweedie(),
  sigma_O = b$estimate[b$term == "sigma_O"],
  phi = b$estimate[b$term == "phi"],
  tweedie_p = b$estimate[b$term == "tweedie_p"],
  range = b$estimate[b$term == "range"],
  fixed_re = list(omega_s = pars$omega_s),
  B = c( # see 'model.matrix' part above
    b$estimate[b$term == "(Intercept)"], # intercept (no RCA, early survey expected value)
    0, # difference from early to late survey non-RCA
    0, # difference from non-RCA to RCA in early survey
    b$estimate[b$term == "poly(depth, 2)1"],
    b$estimate[b$term == "poly(depth, 2)2"],
    0.25 # difference (in log space) from non-RCA to RCA accounting for the early-late survey effect
  ),
  seed = 123
)

head(s)

# expected values:
ggplot(s, aes(X, Y, colour = mu)) +
  scale_colour_viridis_c(trans = "sqrt") +
  facet_wrap(vars(later_survey)) +
  geom_point()

# observed values including zeros:
# (fit to these)
ggplot(s, aes(X, Y, colour = observed)) +
  scale_colour_viridis_c(trans = "sqrt") +
  facet_wrap(vars(later_survey)) +
  geom_point()