If the code in this vignette has not been evaluated, a rendered version is available on the documentation site under ‘Articles’.
library(ggplot2)
library(dplyr)
library(sdmTMB)
library(rstan) # for plot() method
options(mc.cores = parallel::detectCores()) # use rstan parallel processingBayesian estimation is possible with sdmTMB by passing fitted models
to tmbstan::tmbstan() (Monnahan
& Kristensen 2018). All sampling is then done using Stan
(Stan Development Team 2021), and output
is returned as a stanfit object.
Why might you want to pass an sdmTMB model to Stan?
- to obtain full posterior inference on parameters
- to avoid the Laplace approximation on the random effects
- to quantify uncertainty in derived quantities that are not already calculated in the model
- in some cases, models that struggle to converge with maximum likelihood can be sampled adequately with MCMC given carefully chosen priors (e.g., Monnahan et al. 2021)
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:
Adding priors
In that first model fit we did not use any priors (with maximum
likelihood estimation, these can also be thought of as penalties on the
likelihood). For this first model, the priors are implied as uniform on
the internal parameter space. However, sdmTMB provides the option of
applying additional priors. Here we will show an example of applying a
Normal(0, 5) prior (mean, SD) on the intercept and a Normal(0, 2) prior
on the slope parameter. We could guess the model-matrix structure from
our formula, but we can verify it by looking at the internal model
matrix from the previous fit (do_fit = FALSE would save
time if you did not want to fit it the first time).
head(fit$tmb_data$X_ij[[1]])Each column corresponds to the order of the b
priors:
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 #<
)Jacobian adjustments
Adding priors/penalties and fixing spatial parameters are strategies that can help with successful maximum likelihood estimation. If we want to do Bayesian sampling, we need to make one more adjustment to our function call: accounting for non-linear parameter transformations with Jacobian adjustments.
What are Jacobian adjustments and why do we need them? Jacobian
adjustments are necessary when model parameters are transformed in a way
that changes their scale or implied distribution. A good example is the
estimation of variance parameters. Whether we are interested in spatial,
spatiotemporal, or residual variation, the quantity of interest is
usually the variance or standard deviation
.
These quantities are constrained to be greater than 0, so a widely used
estimation strategy is to estimate them in log space, which is
unconstrained. With ln_sigma estimated,
sigma = exp(ln_sigma) can be calculated internally and used
to calculate the likelihood. There are a number of helpful references
detailing the mathematics behind this, including the Stan
manual. Without equations, the Jacobian adjustment can be thought of
as properly stretching the posterior distribution to account for the
change of variables.
In sdmTMB, we can turn these Jacobian adjustments on
with the flag bayesian = TRUE. Applying this to our
fit_mle object,
fit_bayes <- update(fit_mle,
bayesian = TRUE
)It is important to emphasize that this bayesian flag
needs to be enabled for any model passed to Stan; MCMC estimation
without it will generally lead to biased parameter estimates on the
transformed scale.
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_bayes$tmb_obj,
iter = 1000, chains = 2,
seed = 8217 # ensures repeatability
)Sometimes you may need to adjust the sampler settings such as:
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_stanThe Rhat values look reasonable (< 1.05). The
n_eff (number of effective samples) values also 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:
Posterior predictive checks
We can perform posterior predictive checks to assess whether our
model can generate 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. We need to ensure
nsim is less than or equal to the total number of
post-warmup samples.
set.seed(19292)
samps <- sdmTMBextra::extract_mcmc(fit_stan)
s <- simulate(fit_mle, mcmc_samples = samps, 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 the
posterior samples to the mcmc_samples argument in
predict.sdmTMB().
pred <- predict(fit_mle, mcmc_samples = samps)The output is a matrix where each row corresponds to a predicted data point and each column corresponds to a posterior 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, mcmc_samples = samps)
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:
As an example of calculating a derived parameter, here we will calculate the marginal spatial random-field standard deviation:
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,
mcmc_samples = samps,
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()