This vignette describes gratia’s tools for posterior simulation, a powerful way of doing inference using estimated generalized additive models.
We’ll need the following packages for this vignette:
pkgs <- c(
"mgcv", "gratia", "dplyr", "tidyr", "ggplot2", "ggdist",
"distributional", "tibble", "withr", "patchwork"
)
vapply(pkgs, library, logical(1L), logical.return = TRUE, character.only = TRUE)
#> mgcv gratia dplyr tidyr ggplot2
#> TRUE TRUE TRUE TRUE TRUE
#> ggdist distributional tibble withr patchwork
#> TRUE TRUE TRUE TRUE TRUEA generalized additive model (GAM) has the following form
g(\mu_i) = \eta_i = A_i + \mathbf{X}_i \boldsymbol{\gamma} + \sum_{j=1}^{J} f_j(x_{ij}), \; y_i \sim \mathcal{D}(\mu_i, \phi)
where g() is the link function, A is an offset, \mathbf{X}_i is the ith row of a parametric model matrix, \boldsymbol{\gamma} is a vector of parameters for the parametric terms, f_j is a smooth function of covariate x_j. y_i \sim \mathcal{D}(\mu_i, \phi) denotes that the observations y_i are distributed as some member of the exponential family of distributions with mean \mu_i and scale parameter \phi.
The smooth functions f_j are represented in the model via penalised splines basis expansions of the covariates, that are a weighted sum of basis functions
f_j(x_{ij}) = \sum_{k=1}^{K} \beta_{jk} b_{jk}(x_{ij})
where \beta_{jk} is the weight (coefficient) associated with the kth basis function b_{jk}() evaluated at the covariate value x_{ij} for the jth smooth function f_j. Wiggliness penalties \sum_j \lambda_j \boldsymbol{\beta}^{\mathsf{T}} \mathbf{S}_j \boldsymbol{\beta} controls the degree of smoothing applied to the f_j through the smoothing parameters \lambda_j.
Having fitted the GAM in mgcv using REML or ML smoothness selection, we obtain a vector of coefficients \hat{\boldsymbol{\beta}}, including the parametric coefficients \boldsymbol{\gamma}. Under the Bayesian view of smoothing (Miller, 2025), we use a Gaussian approximation to their posterior:
\boldsymbol{\beta} \mid \mathbf{y}, \mathbf{B}, \hat{\boldsymbol{\lambda}}, \hat{\phi} \;\overset{\mathrm{approx}}{\sim}\; \mathrm{MVN}(\hat{\boldsymbol{\beta}}, \mathbf{V}_{\mathrm{b}}).
Here \mathbf{y} is the observed response vector and \mathbf{B} is the full model matrix, including the parametric terms and spline bases. The chosen smoothing parameters \hat{\boldsymbol{\lambda}} and the response scale \hat{\phi}, when estimated, are held fixed. \mathbf{V}_{\mathrm{b}} is the Bayesian coefficient covariance matrix, rather than the frequentist covariance matrix also available from mgcv.
For a Gaussian identity-link GAM with a quadratic smoothing penalty and no coefficient inequality constraints, the Gaussian likelihood and the corresponding Gaussian smoothing prior give an exactly multivariate-normal coefficient posterior when the response variance \sigma^2 = \phi and all smoothing parameters are fixed, provided the posterior is proper. Using fitted values of these quantities gives this conditional posterior; it does not integrate over their uncertainty. For other response distributions or links, normality is generally an approximation. The Metropolis Hastings example below shows a case where that approximation performs poorly.
What are we simulating?
Posterior simulation involves randomly sampling from \text{MVN}(\hat{\boldsymbol{\beta}}, \mathbf{V}_{\text{b}}) or \mathcal{D}(\mu_i, \phi), or both.
We might simulate from the posterior distribution of a single estimated smooth function to see the uncertainty in the estimate of that function. To do this we simulate for just a subset of \boldsymbol{\beta}, \beta_{j}, associated with the f_j of interest. Instead, we might be interested in the uncertainty in the expectation (expected value) of the model at some given values of the covariates, in which case we can simulate for all \boldsymbol{\beta} to sample from the posterior of \mathbb{E}(y_i), the fitted values of the model. Or we might want to generate new values of response variable via draws from the conditional distribution of the response, by simulating new response data \mathbf{y}^{\ast}, at either the observed \mathbf{x} or new values \mathbf{x}^{\ast}, from y^{\ast}_i | \boldsymbol{\eta}, \mathbf{x} \sim \mathcal{D}(\hat{\mu_i}, \phi). Finally, we can combine posterior simulation from both distributions to generate posterior draws for new data \mathbb{y}^{\ast} that also include the uncertainty in the expected values.
gratia has functionality for all these options through the following functions
-
smooth_samples()generates draws from the posterior distribution of single estimated smooth functions, -
fitted_samples()generates joint draws of distribution parameters; usequantity = "expectation"for draws of \mathbb{E}(y_i | \mathbf{X}_i = x_i), the expected response, -
predicted_samples()generates new response data given supplied values of covariates y^{\ast}_i | \mathbf{X}_i = x^{\ast}_i -
posterior_samples(), generates draws from the posterior distribution of the model, including the uncertainty in the estimated parameters of that model.
With quantity = "expectation",
fitted_samples() describes uncertainty in the expected
response at the supplied covariates. It adds no response noise.
posterior_samples() draws observations after drawing model
coefficients, so it includes response variation as well as coefficient
uncertainty. predicted_samples() draws observations with
the coefficients held at their fitted values.
By default, coefficient draws use \mathbf{V}_{\mathrm{b}}, holding the selected
smoothing parameters fixed. For many mgcv models fitted with REML or ML,
a corrected covariance matrix \mathbf{V}_{\mathrm{c}} is also available.
Setting unconditional = TRUE uses this correction to
approximately account for smoothing-parameter uncertainty. It does not
draw smoothing parameters from their posterior or include uncertainty in
the response scale.
There are two additional functions available in gratia that do posterior simulation:
gratia provides simulate() methods for models estimated
using gam(), bam(), and gamm(),
as well as those fitted via scam() in the scam package.
simulate() is a base R convention that does the same thing
as predicted_samples(), returning simulated responses as a
matrix. derivative_samples() draws derivatives of the
expected response with respect to a focal covariate. It evaluates
expectations at neighbouring covariate values using the same coefficient
draws, preserving their dependence. We’ll reserve discussion of
derivatives for a separate vignette focused on estimating derivatives
from GAMs.
In the following sections we’ll look at each of the four main posterior simulation functions in turn.
Posterior smooths and smooth_samples()
For a selected smooth f_j, we retain
its coefficients \boldsymbol{\beta}_j^{\ast} from each joint
coefficient draw and evaluate \mathbf{f}_j^{\ast} =
\mathbf{B}_j\boldsymbol{\beta}_j^{\ast}, where \mathbf{B}_j contains the smooth’s basis
functions at the requested covariate values. Under the Gaussian
approximation, the marginal distribution of these coefficients is \mathrm{MVN}(\hat{\boldsymbol{\beta}}_j,
\mathbf{V}_{\mathrm{b},jj}), where \mathbf{V}_{\mathrm{b},jj} is their block of
the full covariance matrix. smooth_samples() uses the same
joint coefficient draws for all selected smooths, preserving dependence
between their sampled contributions.
To illustrate smooth uncertainty, we simulate independent covariates
x0 to x3 on the unit interval. The Gaussian
response combines nonlinear effects of x0, x1
and x2 with independent Gaussian noise; x3 has
no effect. This is Gu & Wahba’s example, available through
data_sim().
ss_df <- data_sim("eg1", seed = 42)We fit a GAM with a smooth for each covariate, selecting smoothing parameters by REML.
The contribution of f(x_0) uses coefficients 2 through 10 of this model. The other coefficients are not included when evaluating this smooth contribution.
s_x0 <- get_smooth(m_ss, "s(x0)")
smooth_coef_indices(s_x0)
#> [1] 2 3 4 5 6 7 8 9 10To sample from the posterior distribution of these coefficients we
use smooth_samples() choosing the particular smooth we’re
interested in using the select argument; if we want to
sample smooths from the posteriors of all smooths in a model, then
select can be left at its default value.
sm_grid <- data_slice(m_ss, x0 = evenly(x0, n = 100))
sm_samp <- smooth_samples(m_ss, select = "s(x0)", data = sm_grid,
n_draws = 100, seed = 21)Here we supply 100 evenly spaced values of x0 through
data so that we can reuse exactly the same grid for the
interval below. Without data, n_grid controls
the number of covariate evaluation points constructed by
smooth_samples(). n_draws controls the number
of coefficient draws and hence sampled smooths; here we request 100.
Objects returned by smooth_samples() have a
draw() method available for them
To only draw some of the posterior smooths you can set
n_samples which will randomly select that many smooths to
draw (a seed can be provided via argument seed to make the
set of chosen smooths repeatable.)
We now compare the sampled contribution of s(x0) with
its pointwise 95% credible interval on the same grid.
smooth_samples() draws only the selected smooth
contribution. We therefore set overall_uncertainty = FALSE
in smooth_estimates() to omit overall-mean uncertainty from
the interval too.
# evaluate the fitted smooth over x0 and add on a credible interval
sm_est <- smooth_estimates(m_ss, select = "s(x0)", data = sm_grid,
overall_uncertainty = FALSE) |>
add_confint()
# plot the smooth, credible interval, and posterior smooths
sm_est |>
ggplot(aes(x = x0)) +
geom_lineribbon(aes(ymin = .lower_ci, ymax = .upper_ci),
orientation = "vertical", fill = "#56B4E9", alpha = 0.5
) +
geom_line(
data = sm_samp,
aes(y = .value, group = .draw), alpha = 0.2
) +
geom_line(aes(y = .estimate), linewidth = 1, colour = "#E69F00") +
labs(y = smooth_label(s_x0))
At a fixed value of x_0, the blue interval has approximately 95% posterior probability for the smooth contribution, conditional on the observed data, model design and fitted smoothing and variance parameters. With many posterior draws, about 95% of the sampled values at that particular x_0 would fall inside it. This does not mean that 95% of complete sampled curves stay inside the interval over the whole grid.
Marra & Wood (2012) study a different property: repeated-sampling coverage averaged over the function. Imagine repeatedly generating data from the same true function, refitting the GAM and constructing intervals. For each fit, record the proportion of covariate values at which the interval contains the true function, then average that proportion over repeated data sets. Their theory and simulations show coverage close to the nominal level in many settings, provided heavy oversmoothing is avoided. Coverage can be poor when a centred component is close to a straight line or plane; including overall-mean uncertainty can improve it. These results do not promise 95% frequentist coverage at every covariate value. The component-only interval here is chosen to match the sampled smooth contribution.
Posterior fitted values via fitted_samples()
In this section we request quantity = "expectation"
explicitly. For an ordinary mean-based GAM, these are draws of the same
expected response estimated by predict(). Distributional
models can have several predictors, and their first parameter need not
be the expected response. The default
quantity = "parameters" returns those parameters instead;
family_parameters() describes their meanings.
In this example, using data_sim() we simulate data from
example 6 of Luo & Wahba (1997)
\sin(2 \cdots(4x - 2)) + 2 \cdots\exp(-256 \cdots(x - 0.5)^2)
f <- function(x) {
sin(2 * ((4 * x) - 2)) + (2 * exp(-256 * (x - 0.5)^2))
}
df <- data_sim("lwf6", dist = "normal", scale = 0.3, seed = 2)
plt <- df |>
ggplot(aes(x = x, y = y)) +
geom_point(alpha = 0.5) +
geom_function(fun = f)
plt
To these data we fit an adaptive smoother
Next we create a data slice of 200 values over the interval (0,1) at which we’ll predict from the model and generate posterior fitted values for
new_df <- data_slice(m, x = evenly(x, lower = 0, upper = 1, n = 200)) |>
mutate(.row = row_number())then we compute the fitted values for the new data
fv <- fitted_values(m, quantity = "expectation", data = new_df)The posterior fitted values are drawn with
fitted_samples() using a Gaussian approximation to the
posterior. Here we just take 10 draws from the posterior for each
observation in new_df and merge the posterior draws with
the data
fs <- fitted_samples(m, quantity = "expectation", data = new_df, n_draws = 10, seed = 4) |>
left_join(new_df |> select(.row, x), by = join_by(.row == .row))Adding the posterior fitted samples to the plot of the data, superimposing the Bayesian credible interval on the fitted values
plt +
geom_ribbon(data = fv, aes(y = .fitted, ymin = .lower_ci, ymax = .upper_ci),
fill = "red", alpha = 0.3) +
geom_line(data = fs, aes(group = .draw, x = x, y = .fitted),
colour = "yellow", alpha = 0.4)
we see the posterior draws are largely contained the credible interval.
The difference between what we did here and what we did with
smooth_samples() is that now we’re including the effects of
all the other model terms. In this simple model with a single smooth and
an identity link, the only difference is that the model constant term
and its uncertainty is included in the samples.
Additional examples
Distributional GAMs
Where possible, predicted_samples() and
posterior_samples() also work for distributional GAMs. This
is possible when a suitable random number generator is available in the
family object stored within the model. To illustrate, I
reuse an example of Matteo Fasiolo, author of the mgcViz
package, based on the classic mcycle data set from the
MASS package.
We being by loading the data and adding on a row number variable for use late
data(mcycle, package = "MASS")
mcycle <- mcycle |>
mutate(
.row = row_number()
) |>
relocate(.row, .before = 1L)To these data we fit a standard Gaussian GAM for the conditional mean
of accel.
Next, we simulate n_sim new response values for each of
the observed data using predicted_samples()
n_sim <- 10
n_data <- nrow(mcycle)
sim_gau <- predicted_samples(m_gau, n_draws = n_sim, seed = 10) |>
left_join(mcycle |> select(-accel), # join on the observed data for times
by = ".row"
) |>
rename(accel = .response) |> # rename
bind_rows(mcycle |>
relocate(.after = last_col())) |> # bind on observer data
mutate( # add indicator: simulated or observed
type = rep(c("simulated", "observed"),
times = c(n_sim * n_data, n_data)
),
.alpha = rep( # set alpha values for sims & observed
c(0.2, 1), time = c(n_sim * n_data, n_data)
)
)The comments briefly indicate what the dplyr code is doing. Now we can plot the observed and simulated data
plt_labs <- labs(
x = "Time after impact [ms]",
y = "Acceleration [g]"
)
plt_gau <- sim_gau |>
ggplot(aes(x = times)) +
geom_point(aes(y = accel, colour = type, alpha = .alpha)) +
plt_labs +
scale_colour_manual(values = gratia_palette()(9)[c(6, 5)]) +
scale_alpha_identity()The resulting plot is shown in the left-hand panel of the figure
below. There is clearly a problem here; the simulated data don’t look
much like the observations in the 15ms immediately after the impact and
again at ~45ms after impact. This is due to the model we fitted only
being for the conditional mean of accel.
Instead, we model both the conditional mean and the conditional variance of the data, through linear predictors for both parameters of the Gaussian distribution
m_gaulss <- gam(
list(accel ~ s(times, k = 20, bs = "tp"),
~ s(times, bs = "tp")
), data = mcycle, family = gaulss()
)Simulating new data follows using the same code as earlier
sim_gaulss <- predicted_samples(m_gaulss, n_draws = n_sim, seed = 20) |>
left_join(mcycle |> select(-accel),
by = ".row"
) |>
rename(accel = .response) |>
bind_rows(mcycle |> relocate(.after = last_col())) |>
mutate(
type = rep(c("simulated", "observed"),
times = c(n_sim * n_data, n_data)
),
.alpha = rep(c(0.2, 1), time = c(n_sim * n_data, n_data))
)The benefit of all that data wrangling is now realised as we can replace the data in the plot we created earlier with the simulations from the distribution GAM, and then plot it
plt_gaulss <- plt_gau + sim_gaulss
plt_gau + plt_gaulss +
plot_annotation(tag_levels = "a") +
plot_layout(guides = "collect", ncol = 2)
The plot of the simulated response data for the distributional GAM is shown in the right hand panel of the plot. Now, there is much less disagreement between the observed data and that which we can produce from the fitted model.
Prediction intervals
One use for posterior simulation is to generate prediction intervals for a fitted model. Prediction intervals include two sources of uncertainty; that from the estimated model itself, plus the sampling uncertainty or error that arises from drawing observations from the conditional distribution of the response.
For example, in a Gaussian GAM, the first source of uncertainty comes from the uncertainty in the estimates of \beta_j, the model coefficients. This uncertainty is in the mean or expected value of the response. The second source of uncertainty stems from the error term, the estimated variance of the response. These two parameters define the conditional distribution of Y_i. For any value of the covariate(s) \mathbf{X}, our estimated model defines the entire distribution of the response values we might expect to observe at those covariate values.
To illustrate, we’ll fit a simple GAM with a single smooth function
to data simulate from Gu & Wahba’s function f_2 using data_sim(). We
simulate 400 values from a Gaussian distribution with variance \sigma^2 = 1.
df <- data_sim("gwf2", n = 400, scale = 1, dist = "normal", seed = 8)The simulated data, and the true function from which they were generated are shown below
df |>
ggplot(aes(x = x, y = y)) +
geom_point() +
geom_function(fun = gw_f2, colour = "#0072B2", linewidth = 1.5)
A GAM for these data contains a single smooth function of
x
If we consider a new value of the covariate x, x^{\ast} = 0.5, the expected value of the
response given our model, \mathbb{E}(y^{*} | x
= x^{*}), is ~2.92, which we obtain using
predict()
mu <- predict(m, newdata = data.frame(x = 0.5))
mu
#> 1
#> 2.919094This value is the mean of a Gaussian distribution that, if
our model is a correct description of the data, describes the
distribution of the values that Y might
take when x = 0.5. The Gaussian
distribution is defined by two parameters; the mean, \mu, which describes the middle of the
distribution, and the variance, \sigma^2, which describes how spread out the
distribution is about the mean. To fully describe the Gaussian
distribution of the response when x =
0.5, we need an estimate of the variance. We obtain an estimate
from the model’s scale parameter, \phi.
This variance is stored as the element scale in the model
object. We take its square root to obtain the standard deviation
required by dist_normal() below
sigma <- sqrt(m$scale)
sigma
#> [1] 1.009666We can visualise what this distribution looks like with some magic from the ggdist package
df |>
ggplot(aes(x = x, y = y)) +
stat_halfeye(aes(ydist = dist_normal(mean = mu, sd = sigma)),
x = 0.5, scale = 0.2, slab_fill = "#E69F00", slab_alpha = 0.7
) +
geom_point() +
geom_function(fun = gw_f2, colour = "#0072B2", linewidth = 1.5) +
geom_point(x = 0.5, y = mu, colour = "red")
The orange region shows the expected density of response values at
x^{\ast} = 0.5 that our model predicts
we could expect to observe. This region holds the estimated mean and
variance fixed. Prediction intervals combine variation about the
expected value with uncertainty in that expected value. We first use
fitted_values() to obtain a 95% credible interval for the
expected value alone
fitted_values(m, quantity = "expectation", data = data.frame(x = 0.5))
#> # A tibble: 1 × 7
#> .row x .parameter .fitted .se .lower_ci .upper_ci
#> <int> <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0.5 expectation 2.92 0.161 2.60 3.23The .se column is the standard error (standard
deviation) of the estimated value (.fitted), while
.lower_ci and .upper_ci are lower and upper
uncertainty bounds (at the 95% level) on the estimated value
respectively. These coefficient-based calculations hold the response
variance \phi = \sigma^2 at its fitted
value; they do not include uncertainty in that variance.
While it would be pretty easy to compute upper and lower tail
quantiles of the fitted Gaussian distribution for a range of values of
x to get a prediction interval, we’d be ignoring the
uncertainty in the model estimates of the mean. Posterior simulation
provides a simple and convenient way to generate a prediction interval
that includes the model uncertainty, and which works in principle for
any of the families available in mgcv (although in practice, not all
families are currently supported by gratia).
To compute a prediction interval over x for our GAM, we
being by creating a set of data evenly over the range of x
observed in the data used to fit the model
ds <- data_slice(m, x = evenly(x, n = 200)) |>
mutate(.row = row_number())The added variable .row will be used later to match
simulated responses to the rows in ds. We also compute
fitted means and credible intervals at these covariate values. We set
unconditional = TRUE here and in
posterior_samples() below so that both use the same
covariance correction for smoothing-parameter uncertainty.
fv <- fitted_values(m, quantity = "expectation", data = ds, unconditional = TRUE)That step isn’t required in order to do posterior simulation with gratia, but we’ll use the fitted values later to show the model estimated values and their uncertainty in contrast to the prediction interval.
We use posterior_samples() to generate new response data
for each of the new x values in ds and use a
join to add the prediction data to each draw
ps <- posterior_samples(m, n_draws = 10000, data = ds, seed = 24,
unconditional = TRUE) |>
left_join(ds, by = join_by(.row == .row))
ps
#> # A tibble: 2,000,000 × 4
#> .row .draw .response x
#> <int> <int> <dbl> <dbl>
#> 1 1 1 -0.589 0.00129
#> 2 2 1 -0.483 0.00629
#> 3 3 1 -1.27 0.0113
#> 4 4 1 0.399 0.0163
#> 5 5 1 0.0368 0.0213
#> 6 6 1 0.443 0.0263
#> 7 7 1 -0.249 0.0313
#> 8 8 1 -0.406 0.0363
#> 9 9 1 0.683 0.0413
#> 10 10 1 -0.417 0.0463
#> # ℹ 1,999,990 more rowsHere we asked for 10000 posterior draws for each new value of
x. Ideally we’d generate at least three or four times this
many draws to get a more precise estimate of the prediction interval,
but we keep the number low in this vignette to avoid excessive
computation time. We’re also using the smoothness parameter selection
corrected version of the Bayesian covariance matrix; this matrix has
been adjusted to account for us not knowing the value of the smoothing
parameter for f(x_i).
ps is a tibble, with n_draws * nrow(ds)
(here, 10000 * nrow(ds)) rows. The .draw
variable groups the simulated values by posterior draw, while
.row groups posterior draws for the same value of
x. To summarise the posterior draws using {dplyr} we need a
function that will compute quantiles of the posterior predictive
distribution for each value of x (each .row).
The following function is a simple wrapper around the
quantile() function from base R, which arranges the output
from quantile() as a data frame.
quantile_fun <- function(x, probs = c(0.025, 0.5, 0.975), ...) {
tibble::tibble(
.value = quantile(x, probs = probs, ...),
.q = probs * 100
)
}We apply this function to the simulated responses, grouping by
.row to summarise the posterior predictive distribution
separately at each new value of x. reframe()
applies quantile_fun() within each group. We then pivot the
summary from long to wide format and join the covariate values using
.row.
p_int <- ps |>
group_by(.row) |>
reframe(quantile_fun(.response)) |>
pivot_wider(
id_cols = .row, names_from = .q, values_from = .value,
names_prefix = ".q"
) |>
left_join(ds, by = join_by(.row == .row))
p_int
#> # A tibble: 200 × 5
#> .row .q2.5 .q50 .q97.5 x
#> <int> <dbl> <dbl> <dbl> <dbl>
#> 1 1 -2.93 -0.859 1.17 0.00129
#> 2 2 -2.68 -0.644 1.38 0.00629
#> 3 3 -2.47 -0.414 1.62 0.0113
#> 4 4 -2.24 -0.236 1.80 0.0163
#> 5 5 -2.03 -0.0245 2.02 0.0213
#> 6 6 -1.83 0.195 2.20 0.0263
#> 7 7 -1.58 0.376 2.41 0.0313
#> 8 8 -1.43 0.601 2.60 0.0363
#> 9 9 -1.16 0.843 2.80 0.0413
#> 10 10 -0.975 1.04 3.05 0.0463
#> # ℹ 190 more rowsThe printed tibble shows the 95% pointwise prediction interval for
the first 10 rows of ds. .q50 is the median of
the posterior predictive distribution at each value of
x.
We can now use the various objects we have produced to plot the
fitted values from the model (and their uncertainties), as well as the
prediction intervals we just generated. We add the observed data used to
fit the model as black points, and summarise the posterior samples (from
ps) using a hexagonal binning (to avoid plotting all 2
million posterior samples)
fv |>
ggplot(aes(x = x, y = .fitted)) +
# summarise the posterior samples
geom_hex(
data = ps, aes(x = x, y = .response, fill = after_stat(count)),
bins = 50, alpha = 0.7
) +
# add the lower and upper prediction intervals
geom_line(data = p_int, aes(y = .q2.5), colour = "#56B4E9",
linewidth = 1.5) +
geom_line(data = p_int, aes(y = .q97.5), colour = "#56B4E9",
linewidth = 1.5) +
# add the lower and upper credible intervals
geom_line(aes(y = .lower_ci), colour = "#56B4E9", linewidth = 1) +
geom_line(aes(y = .upper_ci), colour = "#56B4E9", linewidth = 1) +
# add the fitted model
geom_line() +
# add the observed data
geom_point(data = df, aes(x = x, y = y)) +
scale_fill_viridis_c(option = "plasma") +
theme(legend.position = "none") +
labs(y = "Response")
The outer pair of blue lines shows the 95% pointwise posterior
predictive interval for a new response at each value of x.
It combines coefficient uncertainty with response variation, while
holding the response variance at its fitted value. The inner pair shows
uncertainty in the fitted mean alone. The hexagonal bins summarise the
simulated responses, with warmer colours indicating larger counts.
Metropolis Hastings sampler
In some cases, the Gaussian approximation to the posterior
distribution of the model coefficients can fail. Simon Wood shows an
example in ?gam.mh, where it performs poorly for a binomial
GAM with large numbers of zeroes. mgcv::gam.mh() alternates
independent multivariate t proposals centred on the fitted
coefficients with Gaussian random walk proposals scaled from the
approximate coefficient covariance. Its acceptance step evaluates the
conditional posterior using the penalised likelihood, with smoothing and
other family parameters held fixed.
In this section, I rework Simon’s example of the failure of the
Gaussian approximation from ?gam.mh to show how to use
gratia to generate posterior draws using the Metropolis
Hastings sampler provided by gam.mh().
We begin by defining a function that will simulate data for the example.
ga_fail <- function(seed) {
df <- tibble(y = c(
rep(0, 89), 1, 0, 1, 0, 0, 1, rep(0, 13), 1, 0, 0, 1,
rep(0, 10), 1, 0, 0, 1, 1, 0, 1, rep(0, 4), 1, rep(0, 3),
1, rep(0, 3), 1, rep(0, 10), 1, rep(0, 4), 1, 0, 1, 0, 0,
rep(1, 4), 0, rep(1, 5), rep(0, 4), 1, 1, rep(0, 46)
)) |>
mutate(
x = withr::with_seed(
seed,
sort(c(0:10 * 5, rnorm(length(y) - 11) * 20 + 100))
),
.row = row_number()
) |>
relocate(.row, .before = 1L)
df
}Which we use to simulate a data set and plot it
df <- ga_fail(3)
df |>
ggplot(aes(x = x, y = y)) +
geom_point()
Note how there are very few ones (successes) and that for large parts of the covariate space the response consists of only zeroes (failures).
We fit a binomial (logistic) GAM to the data
and then generate sample from the posterior distribution using the default Gaussian approximation and subsequently using the simpler Metropolis Hastings sampler.
fs_ga <- fitted_samples(m_logit, quantity = "expectation", n_draws = 2000, seed = 2)
fs_mh <- fitted_samples(m_logit, quantity = "expectation",
n_draws = 2000, seed = 2, method = "mh", thin = 2,
rw_scale = 0.4
)The method argument is used to select the Metropolis
Hastings sampler, and we specify two additional arguments:
-
thin, which controls how many draws are skipped between each retained sample, and -
rw_scale, which is the scaling factor by which the posterior covariance matrix is shrunk for the random walk proposals.
We leave the two other important arguments at their defaults:
-
burnin = 1000, the number of samples to discard prior to sampling, and -
t_df = 40, the degrees of freedom for the t proposals.
With t_df = 40, the independent t proposals are
close to Gaussian. The acceptance step still targets the conditional
posterior, so nearly Gaussian proposals do not require the accepted
draws to follow the Gaussian approximation.
Having collected the posterior draws, we summarise each set into 50%,
80%, and 95% intervals using ggdist::median_qi(), and add
on the data locations with a left join
excl_col <- c(".draw", ".parameter", ".row")
int_ga <- fs_ga |>
group_by(.row) |>
median_qi(.width = c(0.5, 0.8, 0.95), .exclude = excl_col) |>
left_join(df, by = join_by(.row == .row))
int_mh <- fs_mh |>
group_by(.row) |>
median_qi(.width = c(0.5, 0.8, 0.95), .exclude = excl_col) |>
left_join(df, by = join_by(.row == .row))First we plot the intervals for the Gaussian approximation to the posterior, and then we repeat the plot using the intervals derived from the Metropolis Hastings sampler, arranging the two plots using patchwork
plt_ga <- df |>
ggplot(aes(x = x, y = y)) +
geom_point() +
geom_lineribbon(
data = int_ga,
aes(x = x, y = .fitted, ymin = .lower, ymax = .upper)
) +
scale_fill_brewer() +
labs(title = "Gaussian approximation")
plt_mh <- df |>
ggplot(aes(x = x, y = y)) +
geom_point() +
geom_lineribbon(
data = int_mh,
aes(x = x, y = .fitted, ymin = .lower, ymax = .upper)
) +
scale_fill_brewer() +
labs(title = "Metropolis Hastings sampler")
plt_ga + plt_mh + plot_layout(guides = "collect")
The Gaussian approximation-based intervals are shown on the left of
the figure, which for most of the range of x are largely
useless, covering the entire range of the response, despite the fact
that we only observed zeroes for large parts of the covariate space.
Contrast those intervals with the ones obtained using the Metropolis
Hastings sampler; these intervals much better reflect the uncertainty in
the estimated response as a function of x where the data
are all zeroes.
Multiple chains and sampling diagnostics
n_draws is the total number of retained coefficient
draws across chains. For example, four chains with
n_draws = 1000 retain 250 draws each. Each MH chain has its
own burn-in period. Use n_cores to run chains concurrently
with mirai. gratia reuses an active mirai pool; if none is available, it
creates and removes a separate temporary pool. An existing pool is never
resized or stopped by the sampler.
response_draws <- posterior_samples(
m_ss, n_draws = 1000, chains = 4, n_cores = 4,
method = "mh", burnin = 1000, seed = 281,
diagnostics = TRUE
)
sampling_diagnostics(response_draws)The returned table includes .chain and
.iteration. .draw identifies a draw across all
chains. diagnostics = TRUE requires the posterior package
and attaches coefficient summaries, including R-hat, bulk and tail
effective sample sizes, and Monte Carlo standard errors. The effective
sample sizes measure how much information remains after allowing for
correlation between successive draws. A total of 1,000 retained draws
need not provide the information of 1,000 independent draws.
These summaries describe the coefficient chains before response simulation. Response noise can conceal poor mixing of the coefficient chains. To examine coefficient trace plots and then reuse the same draws for response sampling, generate the coefficients once:
coefficients <- post_draws(
m_ss, n_draws = 1000, chains = 4, n_cores = 4,
method = "mh", burnin = 1000, seed = 281,
diagnostics = TRUE
)
bayesplot::mcmc_trace(as_mcmc_draws(coefficients, format = "array"))
response_draws <- posterior_samples(
m_ss, method = "user", draws = coefficients, seed = 282
)as_mcmc_draws() also converts fitted values, smooths,
derivatives, and response samples. Select the observations or smooths of
interest before converting a large table: each sampled quantity becomes
one variable. Choose "draws_array" or
"draws_df" for posterior, "mcmc.list" for
coda, or "array" for an ordinary iteration by chain by
variable array accepted by bayesplot. Conversions require the same
consecutive retained iterations in each chain. Keep the chain columns
when manipulating sampled tables; ordinary matrix subsetting can remove
the chain mapping attached to coefficient matrices.
With four available workers, elapsed sampling time is approximately the time for one chain’s burn-in plus 250 retained iterations, rather than burn-in plus 1,000 iterations. Worker startup, model transfer, and prediction also take time. If burn-in is 1,000 iterations, the corresponding iteration counts are 1,250 and 2,000, so the ideal sampling speedup is about 1.6 times.
Sampling coefficients from bam fits
mgcv::gam.mh() does not support bam fits. Use
mh_backend = "gratia" for independent-observation bam fits,
or to use the compiled evaluator for a gam:
bam_data <- data_sim("eg1", n = 1000, dist = "poisson", seed = 281)
b <- mgcv::bam(y ~ s(x0) + s(x1), data = bam_data, family = poisson())
response_draws <- posterior_samples(
b, n_draws = 1000, chains = 4, n_cores = 4, method = "mh",
mh_backend = "gratia", diagnostics = TRUE, seed = 281
)This backend currently supports Gaussian identity, Poisson log,
binomial logit, Gamma log, and negative binomial log models. It holds
smoothing parameters, dispersion, and any negative binomial shape at
their fitted values. It rejects AR1 dependence, additional penalties
specified through paraPen or H, and fits made
from pre-fit G objects. For discrete bam fits it evaluates
the fitted bases at the original covariates; it does not use a
compressed, discretised likelihood.
Small model matrices are cached. Large designs are evaluated in
blocks, controlled by block_size, without retaining the
full matrix. This reduces memory use but repeatedly constructs the basis
matrices. The dense_limit argument sets the maximum number
of elements in a cached design matrix; zero forces blockwise
evaluation.
For Gaussian identity-link GAMs with no coefficient inequality
constraints, a quadratic smoothing penalty, fixed response variance and
fixed smoothing parameters, the proper coefficient posterior is
Gaussian, as described above. method = "gaussian" therefore
provides direct independent draws from that conditional posterior. Its
unconditional = TRUE covariance correction can
approximately account for smoothing-parameter uncertainty when mgcv
supplies it; it does not integrate over response-variance uncertainty.
The MH backends also hold smoothing parameters and dispersion fixed.