If the code in this vignette has not been evaluated, a rendered version is available on the documentation site under 'Articles'.
dplyr_installed <- require("dplyr", quietly = TRUE) ggplot_installed <- require("ggplot2", quietly = TRUE) tmbstan <- requireNamespace("tmbstan", quietly = TRUE) rstan <- requireNamespace("rstan", quietly = TRUE) bayesplot <- requireNamespace("bayesplot", quietly = TRUE) pkgs <- dplyr_installed && ggplot_installed && tmbstan && rstan && bayesplot EVAL <- identical(Sys.getenv("NOT_CRAN"), "true") && pkgs knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.asp = 0.618, eval = EVAL, purl = EVAL )
library(ggplot2) library(dplyr) library(sdmTMB) library(rstan) # for plot() method options(mc.cores = parallel::detectCores()) # use rstan parallel processing
Bayesian estimation is possible with sdmTMB by passing fitted models to tmbstan::tmbstan()
[@monnahan2018].
All sampling is then done using Stan [@standev_2021], and output is returned as a stanfit
object.
Why might you want to pass an sdmTMB model to Stan?
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()
First, fit a spatial random field GLMM with maximum likelihood:
fit <- sdmTMB( observed ~ a1, data = sim_dat, mesh = mesh, family = gaussian(), spatial = "on" ) fit
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
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_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 #< )
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)
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.
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.
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 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, 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()
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)
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()
# fit <- sdmTMB( # present ~ 0 + as.factor(year) + depth_scaled + depth_scaled2, # data = pcod, # mesh = pcod_mesh, # family = binomial(link = "logit"), # priors = sdmTMBpriors( # matern_s = pc_matern(range_gt = 5, sigma_lt = 1) # ) # )
# fit_stan <- tmbstan::tmbstan(fit$tmb_obj, iter = 1000, chains = 4)
# plot(fit_stan, pars = c("ln_tau_O", "ln_tau_E", "ln_kappa")) # ggsave("vignettes/bayes_fig_2.png", height = 7, width = 7, dpi = 72)
# knitr::include_graphics("vignettes/bayes_fig_2.png")
# pars <- rstan::extract(fit_stan)
# p <- bayesplot::mcmc_hex(fit_stan, pars = c("ln_tau_O", "ln_tau_E")) + # bayesplot::plot_bg(fill = "gray95") + # bayesplot::panel_bg(fill = "gray70")
# p <- bayesplot::mcmc_hex(fit_stan, pars = c("ln_tau_O", "ln_tau_E")) + # bayesplot::plot_bg(fill = "gray95") + # bayesplot::panel_bg(fill = "gray70") # p # ggsave("vignettes/bayes_fig_3.png", height = 7, width = 7)
# knitr::include_graphics("vignettes/bayes_fig_3.png")
# p <- bayesplot::mcmc_trace(fit_stan, pars = c("ln_tau_O", "ln_tau_E")) # p
# p <- bayesplot::mcmc_trace(fit_stan, pars = c("ln_tau_O", "ln_tau_E")) # p # ggsave("vignettes/bayes_fig_4.png", height = 7, width = 7)
# knitr::include_graphics("vignettes/bayes_fig_4.png")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.