Simulating data

Here we will demonstrate using a simulated dataset.

set.seed(123)
predictor_dat <- data.frame(
  X = runif(500), Y = runif(500),
  a1 = rnorm(500)
)
mesh <- make_mesh(predictor_dat, xy_cols = c("X", "Y"), cutoff = 0.1)
# plot(mesh)
# mesh$mesh$n
sim_dat <- sdmTMB_simulate(
  formula = ~a1,
  data = predictor_dat,
  mesh = mesh,
  family = gaussian(),
  range = 0.3,
  phi = 0.2,
  sigma_O = 0.2,
  seed = 123,
  B = c(0.8, -0.4) # B0 = intercept, B1 = a1 slope
)

Visualize our simulated data:

ggplot(sim_dat, aes(X, Y, colour = observed)) +
  geom_point() +
  scale_color_viridis_c()

Fitting the model with marginal likelihood

First, fit a spatial random field GLMM with maximum likelihood:

fit <- sdmTMB(
  observed ~ a1,
  data = sim_dat,
  mesh = mesh,
  family = gaussian(),
  spatial = "on"
)
fit

Adding priors

In that first model fit we did not use any priors. In that case, the priors are implied as uniform on the internal parameter space. However, sdmTMB provides the option of applying priors. Here we will show an example of applying a Normal(0, 5) (mean, SD) prior on the intercept and a Normal(0, 1) prior on the slope parameter. We could guess at the model matrix structure based on our formula, but we can verify it by looking at the internal model matrix from the previous fit (using do_fit = FALSE would save time if you didn’t want to fit it the first time).

head(fit$tmb_data$X_ij[[1]])

Each column corresponds to the order of the b priors:

fit <- sdmTMB(
  observed ~ a1,
  data = sim_dat,
  mesh = mesh,
  family = gaussian(),
  spatial = "on",
  priors = sdmTMBpriors(
    # location = vector of means; scale = vector of standard deviations:
    b = normal(location = c(0, 0), scale = c(5, 2)),
  )
)
fit

Fixing a spatial correlation parameter to improve convergence

Sometimes some of the spatial correlation parameters can be challenging to estimate with Stan. One option is to apply penalized complexity (PC) priors with sdmTMBpriors() to the Matérn parameters. Another option, which can also be used in conjunction with the priors, is to fix one or more parameters at their maximum likelihood estimate (MLE) values. Frequently, fixing the parameter ln_kappa can help convergence (e.g., Monnahan et al. 2021). This estimated parameter is transformed into the range estimate, so it controls the rate of spatial correlation decay.

Now we will rebuild the fitted object with fixed (‘mapped’) ln_kappa parameters using the update() function. We’ll use do_fit = FALSE to avoid actually fitting the updated model since it’s not necessary.

# grab the internal parameter list at estimated values:
pars <- sdmTMB::get_pars(fit)
# create a 'map' vector for TMB
# factor NA values cause TMB to fix or map the parameter at the starting value:
kappa_map <- factor(rep(NA, length(pars$ln_kappa)))

# rebuild model updating some elements:
fit_mle <- update(
  fit,
  control = sdmTMBcontrol(
    start = list(
      ln_kappa = pars$ln_kappa #<
    ),
    map = list(
      ln_kappa = kappa_map #<
    )
  ),
  do_fit = FALSE #<
)

Passing the model to tmbstan

Now we can pass the $tmb_obj element of our model to tmbstan::tmbstan(). We are only using 1000 iterations and 2 chains so this vignette builds quickly. In practice, you will likely want to use more (e.g., 2000 iterations, 4 chains).

fit_stan <- tmbstan::tmbstan(
  fit_mle$tmb_obj,
  iter = 1000, chains = 2,
  seed = 8217 # ensures repeatability
)

Sometimes you may need to adjust the sampler settings such as:

tmbstan::tmbstan(
  ...,
  control = list(adapt_delta = 0.9, max_treedepth = 12)
)

See the Details section in ?rstan::stan.

You can also ‘thin’ samples via the thin argument if working with model predictions becomes cumbersome given a large number of required samples.

We can look at the model:

fit_stan

The Rhat values look reasonable (< 1.05). The n_eff (number of effective samples) values mostly look reasonable (> 100) for inference about the mean for all parameters except the intercept (b_j[1]). Furthermore, we can see correlation in the MCMC samples for b_j[1]. We could try running for more iterations and chains and/or placing priors on this and other parameters as described below (highly recommended).

Now we can use various functions to visualize the posterior:

plot(fit_stan)
pars_plot <- c("b_j[1]", "b_j[2]", "ln_tau_O", "omega_s[1]")

bayesplot::mcmc_trace(fit_stan, pars = pars_plot)
bayesplot::mcmc_pairs(fit_stan, pars = pars_plot)

Posterior predictive checks

We can perform posterior predictive checks to assess whether our model can generate predictive data that are consistent with the observations. For this, we can make use of simulate.sdmTMB() while passing in our Stan model. simulate.sdmTMB() will take draws from the joint parameter posterior and add observation error.

set.seed(19292)
s <- simulate(fit_mle, tmbstan_model = fit_stan, nsim = 50)
bayesplot::pp_check(
  sim_dat$observed,
  yrep = t(s),
  fun = bayesplot::ppc_dens_overlay
)

See ?bayesplot::pp_check. The solid line represents the density of the observed data and the light blue lines represent the density of 50 posterior predictive simulations. In this case, the simulated data seem consistent with the observed data.

Plotting predictions

We can make predictions with our Bayesian model by supplying it to the tmbstan_model argument in predict.sdmTMB().

pred <- predict(fit_mle, tmbstan_model = fit_stan)

The output is a matrix where each row corresponds to a row of predicted data and each column corresponds to a sample.

dim(pred)

We can summarize these draws in various ways to visualize them:

sim_dat$post_mean <- apply(pred, 1, mean)
sim_dat$post_sd <- apply(pred, 1, sd)

ggplot(sim_dat, aes(X, Y, colour = post_mean)) +
  geom_point() +
  scale_color_viridis_c()

ggplot(sim_dat, aes(X, Y, colour = post_sd)) +
  geom_point() +
  scale_color_viridis_c()

Or predict on a grid for a given value of a1:

nd <- expand.grid(
  X = seq(0, 1, length.out = 70),
  Y = seq(0, 1, length.out = 70),
  a1 = 0
)
pred <- predict(fit_mle, newdata = nd, tmbstan_model = fit_stan)

nd$post_mean <- apply(pred, 1, mean)
nd$post_sd <- apply(pred, 1, sd)

ggplot(nd, aes(X, Y, fill = post_mean)) +
  geom_raster() +
  scale_fill_viridis_c() +
  coord_fixed()

ggplot(nd, aes(X, Y, fill = post_sd)) +
  geom_raster() +
  scale_fill_viridis_c() +
  coord_fixed()

Extracting parameter posterior samples

We can extract posterior samples with rstan::extract(),

post <- rstan::extract(fit_stan)

The result is a list where each element corresponds to a parameter or set of parameters:

names(post)
hist(post$b_j[, 1])

As an example of calculating a derived parameter, here we will calculate the marginal spatial random field standard deviation:

ln_kappa <- get_pars(fit_mle)$ln_kappa[1] # 2 elements since 2nd would be for spatiotemporal
ln_tau_O <- post$ln_tau_O
sigma_O <- 1 / sqrt(4 * pi * exp(2 * ln_tau_O + 2 * ln_kappa))
hist(sigma_O)

Extracting the posterior of other predicted elements

By default predict.sdmTMB() returns the overall prediction in link space when a tmbstan model is passed in. If instead we want some other element that we might find in the usual data frame returned by predict.sdmTMB() when applied to a regular sdmTMB model, we can specify that through the sims_var argument.

For example, let’s extract the spatial random field values "omega_s". Other options are documented in ?predict.sdmTMB().

fit_pred <- predict(
  fit_mle,
  newdata = nd,
  tmbstan_model = fit_stan,
  sims_var = "omega_s" #<
)

nd$spatial_rf_mean <- apply(fit_pred, 1, mean)
nd$spatial_rf_sd <- apply(fit_pred, 1, sd)

ggplot(nd, aes(X, Y, fill = spatial_rf_mean)) +
  geom_raster() +
  scale_fill_gradient2() +
  coord_fixed()

ggplot(nd, aes(X, Y, fill = spatial_rf_sd)) +
  geom_raster() +
  scale_fill_viridis_c() +
  coord_fixed()