Migrating posterior sampling calls • gratia

The multiple-chain release gives sampling arguments the same meaning across posterior_samples(), fitted_samples(), smooth_samples(), derivative_samples(), their add_*() helpers, and coefficient draw functions. Obsolete arguments now produce errors with migration guidance. Controls must be named; only the model, or the data and model in add_*() helpers, can be passed by position. This prevents old positional calls from changing meaning.

Draw counts and evaluation points

n_draws always means the total retained draws across all chains. It must be divisible by chains. Generating draws requires an explicit count: the previous defaults ranged from one draw to 10,000 draws. chains = NULL uses one chain, or infers the chains from supplied coefficient draws.

Previous argument	Replacement
`n` in fitted, posterior, predicted, smooth and coefficient sampling	`n_draws`
`n_sim` in derivative sampling	`n_draws`
`n` in derivative sampling	`n_grid`
`n_vals` in smooth sampling	`n_grid`
`object` in derivative sampling	`model`
`object` in `add_*()` sampling helpers	`data`
`freq`	`frequentist`
`newdata`	`data`
`ncores`	`n_cores` for MH chains; `n_threads` for Gaussian draws
`n_cores` for Gaussian sampling	`n_threads`
`term` in smooth sampling	`select`

For example, the previous derivative call

# Previous interface; now produces an informative error.
derivative_samples(model, focal = "x", n = 100, n_sim = 1000)

becomes

# 1,000 retained coefficient draws, evaluated at 100 points.
derivative_samples(model, focal = "x", n_draws = 1000, n_grid = 100)

n_grid = 200 is the common default for smooth and derivative sampling. For smooths, it supplies the resolution to smooth_data(); multidimensional smooths can contain more rows than n_grid. Derivative sampling constructs a slice along focal, holding other variables at representative values. Explicit data determine the evaluation locations instead of n_grid. Use explicit data when the exact grid matters. Fitted and response sampling use fitting observations when data = NULL and do not construct a grid.

Multiple chains and parallel work

samples <- fitted_samples(
  model, n_draws = 1000, chains = 4, n_cores = 4,
  method = "mh", burnin = 1000, thin = 1,
  diagnostics = TRUE, seed = 281
)
diagnostics <- sampling_diagnostics(samples)
diagnostics
# These plots describe the coefficient chains used for the fitted samples.
draw(diagnostics, select = 1:2)
draw(diagnostics, type = "rank", select = 1:2)

This retains 250 draws per chain. Every chain discards its own 1,000 burn-in iterations; parallel chains repeat this warmup work. thin = 2 retains every second iteration after burn-in and still returns n_draws draws in total. n_cores limits concurrent MH chains using mirai, capped at chains. An existing worker pool is reused; otherwise temporary workers are cleaned up automatically. Given the same seed and sampler backend, retained draws agree between serial and parallel runs.

Gaussian sampling uses n_threads to control threads inside mvnfast. It does not run concurrent MH chains. A control explicitly supplied to a method that cannot use it produces an error, including MH controls for Gaussian draws or Gaussian covariance controls for MH. block_size and dense_limit apply only to method = "mh", mh_backend = "gratia".

The default MH backend remains "mgcv". The three options are:

`mh_backend`	Evaluation	Model support
`"mgcv"`	`mgcv::gam.mh()`	Its existing `gam` support
`"gratia"`	Compiled likelihoods and MH loop	Gaussian identity, Poisson log, binomial logit, Gamma log, and negative binomial log `gam` and independent-observation `bam` fits
`"family"`	Compiled MH loop with callbacks to the fitted family	Compatible standard, extended, and general mgcv families, including independent-observation `bam` fits

For example, a multinomial fit can use the family evaluator without a separate compiled multinomial likelihood:

samples <- posterior_samples(
  model, n_draws = 1000, chains = 4, n_cores = 4,
  method = "mh", mh_backend = "family",
  burnin = 1000, diagnostics = TRUE, seed = 281
)

The family backend prepares its full design matrix, fitted smoothing penalty, and proposal factorisation once and reuses them across chains. It calls the fitted family’s dev.resids() and linkinv(), or ll() for multiple predictors, preserving response attributes and coefficient mappings. It retains the full design matrix even for bam: use the gratia backend’s block evaluation when that matrix is too large and its family is supported. Both evaluate bases at the original covariates for discrete bam fits.

The gratia and family evaluators reject GAMM components and AR1 dependence; they do not evaluate fitted observation correlations. They also reject paraPen, H, pre-fit G objects, and singular coefficient covariances. The family backend additionally rejects quasi families, which do not provide a likelihood, and Cox proportional hazards models. All three MH backends condition on fitted smoothing, dispersion, and any family parameters outside the coefficient vector. Parameters represented by coefficients are sampled, including scale or covariance coefficients in general families. unconditional = TRUE applies to Gaussian sampling only.

Reusing coefficient draws

betas <- post_draws(model, n_draws = 1000, chains = 4, seed = 281)
fitted <- fitted_samples(model, method = "user", draws = betas)
smooth <- smooth_samples(model, method = "user", draws = betas)
derivative <- derivative_samples(
  model, focal = "x", method = "user", draws = betas
)

The coefficient matrix determines the count and chain layout. Explicit n_draws or chains must agree with it. A plain matrix without chain metadata represents one chain; chain boundaries are never guessed from the row count. Smooths always transform the same joint coefficient draws. The old rng_per_smooth option is removed because independently redrawing each smooth loses their posterior dependence. The previously ignored mu, sigma, and derivative level arguments are also rejected.

Multiple-chain tables retain .chain, .iteration, and globally unique .draw columns. Single-chain tables keep their previous columns. Attached diagnostics refer to coefficients before transformation or response simulation. as_mcmc_draws() converts coefficient matrices or sampled quantities for posterior or coda diagnostics. Optional packages are required only when their formats or diagnostics are requested.

draw(sampling_diagnostics(x)) uses native gratia plotting code. Trace plots show retained iteration within each chain; type = "rank" compares histogram outlines of ranks pooled across chains, with a reference line for equal counts per bin. Ranks use average ranks for ties and no random jitter. Use select to choose coefficient names or positions. Rank plots require at least two chains. Discarded burn-in is unavailable, and thinning is represented by the retained iteration count. Coefficient matrices provide their own draws; fitted, response, smooth and derivative outputs require diagnostics = TRUE when sampling to retain coefficient draws for these plots. Plotting itself needs neither posterior nor bayesplot. Numerical summaries still use posterior.

Which response sampling function to use

fitted_samples() draws the expected response, smooth_samples() draws individual smooth contributions, and derivative_samples() draws derivatives. posterior_samples() also simulates response variation for each coefficient draw. predicted_samples() simulates responses with fitted parameters held fixed. It uses n_draws for the number of response replicates and rejects coefficient-chain controls; use posterior_samples() to include coefficient uncertainty.

Model-specific capabilities still apply. For example, the native diagnostic samplers for lm and mixed models support Gaussian approximations and fitting observations only, and reject unsupported chain controls. See help("diagnostic_support") for model-specific support. Functions that summarise draws, such as response_derivatives() and conditional_differences(), keep their existing grid and simulation-count arguments; their internal sampling calls use this interface.

Likelihoods and external diagnostic wrappers

log_lik(), loo::loo(), loo::waic() and bayesplot::pp_check() use the same n_draws convention. Likelihood and LOO/WAIC wrappers default to 4,000 draws; bayesplot uses 50 for ordinary checks and 4,000 for leave-one-out plots. Supplied coefficient matrices determine their own counts. Pass chains, diagnostics, method and backend controls directly to log_lik(), or through sampling_args for LOO, WAIC and bayesplot. Likelihood matrices and these diagnostic results retain the coefficient sampling metadata. LOO reads stored chain iterations when estimating relative efficiency, while keeping likelihood and response draws in the same original order. Unlabelled supplied matrices still require explicit chain_id or r_eff for LOO. The old n argument and explicit controls for a different sampler produce errors.

Location-scale-shape models

posterior_samples() now accepts evaluation data for gaulss(), gammals(), gumbls(), gevlss(), shash(), twlss() and ziplss() fits. Supply covariates and formula-offset variables; observed responses are unnecessary. Every predictor uses the same joint coefficient draw. For example:

# Evaluate the first 200 observations, retaining their joint parameter draws.
b <- post_draws(model, n_draws = 1000, chains = 4, seed = 281)
y <- posterior_samples(model, data = evaluation_data[1:200, ],
  method = "user", draws = b, seed = 282)

The functions return different quantities:

Function	Quantity for an LSS model
`fitted_samples(scale = "response")`	Expected response, including every distribution parameter that affects it
`fitted_samples(scale = "linear_predictor")`	First linear predictor, for location
`posterior_samples()`	Responses conditional on each joint coefficient draw
`predicted_samples()`	Responses conditional on the fixed fitted coefficients
`derivative_samples(scale = "response")`	Derivative of the expected response, retaining joint coefficient uncertainty
`derivative_samples(scale = "linear_predictor")`	Derivative of the first predictor
`smooth_samples()`	Selected smooth contribution to its own predictor

The historical .parameter = "location" label in fitted samples remains on both scales. On the response scale, .fitted is E[Y]. In particular, Gumbel, GEV and sinh-arcsinh means generally differ from location; the Gamma LSS mean requires exponentiation and the hurdle Poisson mean uses both predictors. These calculations no longer emit the earlier warning about location needing to equal the mean. A GEV mean does not exist at shape >= 1; requesting that mean or its derivative gives an error, while response draws remain available.

gaulss() represents its second response parameter as inverse SD. Sampling uses its reciprocal and divides SD by the square root of an explicit sampling weight. The native gaulss() likelihood ignores fitting weights; default response weights are ones, preserving the fitted likelihood and existing diagnostic rules. The other LSS response generators currently ignore sampling weights, as their mgcv implementations do. Formula offsets are evaluated at supplied covariates; stored fitting offsets are retained with data = NULL. The installed mgcv version rejects offsets when fitting shash() and drops call-level offsets for general-family fits; gratia cannot reconstruct offsets never applied by the fitter. Prefer supported formula offsets.

Grouped response families

posterior_samples() also accepts fitting-data subsets and new covariates for mgcv::gfam() models. Every row must identify its component family, either through the response-index expression in the supplied data or through family_index. Observed response values are not needed. For example:

posterior_samples(grouped_model, data = evaluation_data,
  family_index = component_index, n_draws = 1000, chains = 4, seed = 281)

The same joint coefficient draw is used for every component. Component-specific links, fitted nuisance parameters and formula offsets are retained. New-data sampling weights default to one; explicit weights override them. Output retains .parameter = "component1", "component2", and so on, together with .row, .draw and any chain columns. Missing predictors, membership or sampling weights produce missing responses. Fitting rows excluded by na.exclude are restored. Smooth samples retain their predictor contributions, and response derivatives use each row’s component mean. For derivatives, membership must be recoverable from data; derivative_samples() does not accept family_index. Ordinal response-scale fitted values and their derivatives remain unsupported; response draws and predictor-scale calculations remain available. Censored grouped components require an observation process for posterior response sampling; use simulate() for latent uncensored draws.

MH coefficient sampling rejects BAM fits with nonzero prescribed smooth contributions, such as a nonzero soap-film boundary. Native fitting and prediction can use different offsets for those fits. Supplied coefficient draws can still be evaluated with the native prediction conventions.

Censoring and survival responses

Censored-family generators return latent uncensored responses; posterior predictive checks need a censoring observation process. Cox survival models lack a response generator for this interface, and expected survival times cannot be obtained from a location parameter alone.