Specifying and checking priors

Why priors get their own object

Every Bayesian fit carries priors, whether you write them down or not. tulpa keeps them in one place: a tulpa_priors object that names a prior for each kind of parameter the engine knows about. You build it with tulpa_priors(), fill the slots with prior_*() constructors, and hand the result to the functions that need it. The same object drives prior predictive simulation, so the priors you check are the priors you fit with.

This vignette walks through the prior families, the defaults, how to move a fixed-effect prior and watch the posterior follow, and how to look at what a prior implies about the data before any of it touches the likelihood. The closing section turns that into a short set of working rules.

Three jobs run through the whole vignette, and it helps to name them up front because the right prior depends on which one you are doing. A fixed-effect coefficient lives on the link scale and can be any real number, so it gets a prior over the whole line. A variance or scale parameter is positive by construction, so it gets a prior on the positive half-line and the question is how fast that prior decays toward zero. A correlation or mixing proportion is trapped in the unit interval, so it gets a prior on (0, 1) and the only real decision is whether you lean it toward one end. tulpa gives you one constructor family for each job, and tulpa_priors() keeps one of each in a single object.

Prior families available

A prior in tulpa is a small tagged list built by one of seven constructors. Each one carries a distribution name and its parameters, and prints in a compact form. The seven cover three jobs: location parameters that range over the whole real line, scale parameters that must stay positive, and proportions bounded in the unit interval.

prior_normal(0, 2.5)
prior_normal(0, 1)

prior_half_normal(1)
prior_half_cauchy(2.5)
prior_gamma(2, 0.1)
prior_exponential(1)

c(half_normal_mean = prior_half_normal(1)$sd * sqrt(2 / pi),
  exponential_mean = 1 / prior_exponential(1)$rate,
  gamma_mean       = prior_gamma(2, 0.1)$shape / prior_gamma(2, 0.1)$rate)

prior_pc(U = 1, alpha = 0.01)
prior_pc(U = 0.5, alpha = 0.05)

bin <- tulpa_family("binomial",
  function(eta, params, n_obs, ...) rbinom(n_obs, 1, plogis(eta[[1]])))
gdat <- data.frame(y = 0, x = rnorm(200), g = factor(rep(1:20, each = 10)))

group_spread <- function(pp) {
  vapply(seq_along(pp$y), function(d) {
    p <- plogis(pp$linpred[[d]][[1]])
    sd(tapply(p, gdat$g, mean))
  }, numeric(1))
}
pp_tight <- prior_predict(y ~ x + (1 | g), family = bin, data = gdat,
  n_draws = 200, priors = tulpa_priors(sigma = prior_pc(1, 0.01)), seed = 7)
pp_wide  <- prior_predict(y ~ x + (1 | g), family = bin, data = gdat,
  n_draws = 200, priors = tulpa_priors(sigma = prior_pc(5, 0.01)), seed = 7)
rbind(tight = quantile(group_spread(pp_tight), c(0.5, 0.9, 0.99)),
      wide  = quantile(group_spread(pp_wide),  c(0.5, 0.9, 0.99)))

prior_beta(1, 1)
prior_beta(2, 2)
prior_beta(5, 2)

Three classes of parameter, three kinds of prior

The table above sorts priors by the slot they fill, but the deeper division is by what the parameter does in the model, because that determines how much the prior should be allowed to say.

Fixed effects sit in the linear predictor and the data usually identify them well once you have a few dozen informative rows. The prior on them is regularisation: a backstop for thin data or collinear predictors, not a structural assumption. A wide normal is the right default precisely because you want the likelihood to win whenever it can.

Variance and scale components are different. They are weakly identified by construction, especially when the grouping factor has few levels or the responses are binary, so the prior carries real weight in the posterior and keeps doing so as the data grow until the number of groups grows with them. A PC prior is the recommended default here because it states its strength as a tail probability and shrinks toward the simpler model, which is the safe direction when the data underdetermine the variance. The half-Cauchy is the fallback when you expect a genuinely large component that the PC penalty would fight.

Spatial and temporal hyperparameters are the most weakly identified of all. A single realisation of a spatial field or a time series carries little information about its range or its autocorrelation, so the prior often shapes the posterior more than the data do. The rho_temporal and rho_spatial slots take beta priors on (0, 1): the default prior_beta(2, 2) on temporal correlation is symmetric and pulls gently off the endpoints, where an AR(1) process is hard to estimate, and the default prior_beta(1, 1) on the BYM2 mixing proportion stays flat because there is rarely prior reason to favour structured over unstructured variation. SPDE and Matern range parameters live on the positive half-line and take PC priors of the same kind as a variance component, set from a range you would call large for your study region. The engine assembles these structured priors from a latent specification through prior_from_spec(), which reads a tulpa_temporal or tulpa_spatial object and returns the matching prior, so the structured paths and the manual constructors stay in step.

Default priors

tulpa_priors() with no arguments returns the package defaults. Print it to see all five slots at once.

tulpa_priors()

The defaults read as follows. Fixed effects get prior_normal(0, 2.5), wide enough to let coefficients roam across the plausible link-scale range without being flat. Random-effect standard deviations get a PC prior with P(sigma > 1) = 0.01, which favours smaller variance components and guards against the variance running off when groups are few. The dispersion phi gets a PC prior with P(phi > 10) = 0.01. Temporal correlation defaults to prior_beta(2, 2), symmetric around 0.5, and the spatial mixing proportion to prior_beta(1, 1), a flat unit prior with no lean toward structured or unstructured variation.

priors_default() prints the same defaults with a sentence of interpretation and the customisation entry point for each slot.

priors_default()

Change one slot by naming it; the rest stay at their defaults.

tulpa_priors(
  beta  = prior_normal(0, 1),
  sigma = prior_pc(U = 0.5, alpha = 0.01)
)

Because each slot validates its argument, passing something that is not a tulpa_prior object fails early with a clear message rather than surfacing deep inside a fit.

The five defaults are deliberately conservative. prior_normal(0, 2.5) on fixed effects keeps coefficients roughly in the link-scale band [-5, 5] without going flat. The PC prior P(sigma > 1) = 0.01 on random-effect standard deviations favours small variance components and protects against the variance running away when groups are few. The PC prior P(phi > 10) = 0.01 on the dispersion regularises toward the Poisson limit. The two beta priors stay neutral on the correlation parameters. Those choices are right for a first fit on a healthy dataset, and the cases for overriding them are specific: a strong external estimate to centre a fixed effect on, a variance you expect to be large enough that the default PC prior would over-shrink, a dispersion you can bound tightly from prior counts, or a temporal correlation you have reason to believe is high. Absent one of those, the defaults are the place to start.

priors_default() also takes a family argument. Pass a tulpa_family object from a model package and it prints only the slots that family uses, with a note on what each parameter means for that likelihood, which is the quick way to check whether a family even has a dispersion or a temporal slot before you try to set it.

These five slots cover the direct / conditioning paths and the ModelData samplers. Two backends sit outside them and follow a different convention. The nested-Laplace integrator’s own scale axes (icar, rw1, rw2, ar1’s tau, iid) carry no hyperprior at all: the grid is uniform in log(theta) and the outer weights are a plain softmax of the log marginal, so the effective prior is flat in log-scale, everywhere, by construction of the integration rather than by a tulpa_priors() slot. [tulpa_re_cov_nested()] and [tulpa_eb()] match that convention by default (hyperprior = "flat"); pass hyperprior = "pc_lkj" to opt into the PC + LKJ prior these two functions can also build via [re_cov_pc_lkj_prior()]. [tulpa_re_cov_gibbs()] cannot go fully flat – its Sigma | b step is a conjugate Inverse-Wishart draw, which needs a proper prior to sample from – so it defaults to the weakest proper choice (prior_df = n_coefs + 1) instead.

Priors on fixed effects

The engine’s Laplace path takes a Gaussian prior on the fixed effects through the beta_prior argument of tulpa(), written as list(mean =, sd =). The sd is required; mean defaults to 0. Each may be a scalar applied to every coefficient or a vector with one entry per coefficient. Leaving beta_prior at NULL keeps a weak built-in prior that barely moves the fit.

A prior earns its keep when the data are thin. Simulate a small Gaussian dataset with a slope of 1.2, the kind of sample where the likelihood alone leaves real uncertainty.

n  <- 25
x  <- rnorm(n)
y  <- 0.5 + 1.2 * x + rnorm(n, sd = 1.5)
df <- data.frame(y = y, x = x)

Fit it once with the default weak prior. The slope estimate is driven by the likelihood.

fit_weak <- tulpa(y ~ x, data = df, family = "gaussian",
                  mode = "laplace", phi = 1.5^2)
coef(fit_weak)

Now impose a tight prior that says the slope is near zero: mean 0 and a small standard deviation. The intercept keeps a loose prior, so only the slope is pulled.

fit_tight <- tulpa(y ~ x, data = df, family = "gaussian",
                   mode = "laplace", phi = 1.5^2,
                   beta_prior = list(mean = c(0, 0), sd = c(10, 0.2)))
coef(fit_tight)

The slope under the tight prior sits well below its weakly informed value: the prior at zero and the likelihood at the data strike a compromise, and the small prior standard deviation gives the prior most of the weight. Put the two fits side by side.

data.frame(
  term  = names(coef(fit_weak)),
  weak  = round(coef(fit_weak), 3),
  tight = round(coef(fit_tight), 3)
)

This is regularisation made explicit. Shrinking a coefficient toward a value accepts a little bias in exchange for lower variance, which is worth doing when the data cannot pin the coefficient down on their own. The credible intervals tell the same story from the uncertainty side.

confint(fit_weak)["x", ]
confint(fit_tight)["x", ]

The tight prior narrows the slope’s interval and shifts it toward zero. A vector sd lets you regularise some coefficients hard while leaving others free, which is the usual pattern: pin nuisance terms, let the effects of interest follow the data.

The amount of movement is set by the ratio of prior to likelihood information, which on this Gaussian fit has a closed form worth carrying in your head. For a single coefficient the posterior mean is a precision-weighted average of the prior mean and the least-squares estimate, with weights 1 / sd_prior^2 and 1 / se_mle^2. When the prior standard deviation is much larger than the coefficient’s standard error the likelihood dominates and the prior is nearly invisible; when it is much smaller the prior wins and the estimate sits near the prior mean. The slope’s standard error here is on the order of the residual standard deviation over sqrt(n) times the predictor spread, a few tenths, so a prior standard deviation of 0.2 is comparable to the data’s own precision and the two split the difference, which is what the table shows. Setting a vector sd is then a per-coefficient choice of where on that continuum to sit.

The fixed-effect prior threads through the plain Laplace and sampler paths. Some structured paths (spatial fields, the SPDE integrator) carry their own built-in fixed-effect prior and will tell you so if you pass beta_prior there.

Prior predictive checks

A prior is a statement about parameters, but its consequences land on the data scale, often somewhere surprising. A prior that looks mild on a log link can imply counts in the millions. The way to catch this is a prior predictive check: draw parameters from the priors, push them through the linear predictor, simulate responses, and look at the responses.

prior_predict() does the drawing. It needs a formula, a tulpa_family that supplies a simulator, the data (for its design and dimensions, not its response), and the priors. The family is a thin object built by tulpa_family(): a name plus a simulate_fn(eta, params, n_obs, ...) that turns a linear predictor into a response. Build a Poisson family whose simulator exponentiates the linear predictor.

pois <- tulpa_family(
  name = "poisson",
  simulate_fn = function(eta, params, n_obs, ...) rpois(n_obs, exp(eta[[1]]))
)

Use a small covariate frame; the response column can be a placeholder because the priors generate the response.

dat <- data.frame(y = rep(0, 60), x = rnorm(60))

Draw from a deliberately vague prior on the fixed effects. A standard deviation of 5 on a log link is enormous: it lets the intercept and slope wander far enough to produce astronomically large counts.

pp_vague <- prior_predict(
  y ~ x, family = pois, data = dat, n_draws = 200,
  priors = tulpa_priors(beta = prior_normal(0, 5)), seed = 1
)
pp_vague

The returned object holds one simulated dataset per draw in pp_vague$y. Look at the largest count any draw produced.

max(vapply(pp_vague$y, max, numeric(1)))

That number is far past anything a real count process would produce. The prior is too vague: it places real mass on data the model should consider impossible. Tighten the fixed-effect prior to a standard deviation of 1 and redraw.

pp_ok <- prior_predict(
  y ~ x, family = pois, data = dat, n_draws = 200,
  priors = tulpa_priors(beta = prior_normal(0, 1)), seed = 1
)
max(vapply(pp_ok$y, max, numeric(1)))

The largest simulated count now sits in a range a count model can take seriously. Quantiles across all draws sharpen the contrast.

vague_all <- unlist(pp_vague$y)
ok_all    <- unlist(pp_ok$y)
rbind(
  vague    = quantile(vague_all, c(0.5, 0.9, 0.99)),
  sensible = quantile(ok_all,    c(0.5, 0.9, 0.99))
)

The median is modest under both priors, but the upper tail under the vague prior is orders of magnitude heavier. The vague prior sits at the right centre and fails on its width: the damage hides in the tail, where a glance at the median would miss it entirely.

plot() overlays the simulated datasets so you can see the spread directly. Capping the counts keeps the vague-prior plot readable, since a handful of draws would otherwise stretch the axis past everything else.

pp_capped <- pp_vague
pp_capped$y <- lapply(pp_vague$y, function(yi) pmin(yi, 200))
plot(pp_capped, max_draws = 40)

plot(pp_ok, max_draws = 40)

The sensible-prior draws cluster in a believable band; the vague-prior draws, even capped, spread across the full range. Reading the $y numbers and the overlay together is the check: if the prior predictive covers values the response could never take, the prior is too vague, and the families above let you pull it in. The same machinery extends to any family by swapping the simulate_fn, and to random-effect spread by setting the sigma prior, which prior_predict() draws and feeds into the group effects.

The PC-prior example earlier produced exactly that kind of random-effect spread, and plotting the two group_spread distributions side by side shows the prior’s effect at a glance. The wide PC prior pushes a long right tail onto the between-group spread; the tight default keeps it compact.

library(ggplot2)
sp <- rbind(
  data.frame(prior = "tight (U=1)", spread = group_spread(pp_tight)),
  data.frame(prior = "wide (U=5)",  spread = group_spread(pp_wide)))
ggplot(sp, aes(spread, fill = prior)) +
  geom_density(alpha = 0.5) +
  labs(x = "SD of group-mean probabilities", y = "Prior predictive density") +
  theme(panel.background = element_rect(fill = "transparent"),
        plot.background  = element_rect(fill = "transparent"))

Reading this plot is the variance-component analogue of the count check: the question is whether the spread the prior expects between groups matches what you would believe before seeing the data, and the wide prior’s mass near the right edge of the unit interval is the warning sign.

Choosing priors

A working order of operations for setting priors in tulpa:

• Start from the defaults and change what you have reason to change. tulpa_priors() is weakly informative across the board: a wide normal on fixed effects, PC priors that shrink variance components, symmetric priors on correlations. For a first fit on a healthy dataset, the defaults are a reasonable place to stand.

• Reach for PC priors on variance components. prior_pc(U, alpha) lets you state the prior as a tail probability, which is the parameter you can defend in a methods section. Set U to a standard deviation you would call large on the link scale and alpha to a small probability of exceeding it. When you expect the occasional genuinely large standard deviation, reach for prior_half_cauchy(): its heavier tail tolerates that case where the firm PC penalty would fight it.

• Make fixed-effect priors carry weight only when the data are thin. With hundreds of informative rows the likelihood dominates and beta_prior barely moves the fit, but with few rows, a poorly separated predictor, or a coefficient you want to regularise, a tight beta_prior becomes the lever the worked example demonstrated, and a vector sd lets you pull selected coefficients hard while leaving the rest free.

• Run a prior predictive check before trusting a prior on the link scale. Link-scale intuition is unreliable. A standard deviation that feels mild can imply impossible data once exponentiated, and the failure shows up in the upper tail rather than the centre, so a summary mean will not warn you. prior_predict() plus a one-line tulpa_family() simulator turns the question into numbers you can read off $y and an overlay you can scan; tighten until the simulated data look like data the process could produce.

• Match the prior family to the parameter’s support. Normal for the real line, a positive prior (PC, half-normal, half-Cauchy, gamma, exponential) for standard deviations and dispersions, beta for proportions; the constructors accept any pairing, so a normal on a variance slot will build cleanly and then bite at fit time. The map in the families section is the safe default to copy.

• Keep the prior object you checked. A single tulpa_priors object feeds both prior_predict() and the fit, so a prior that survives the predictive check is the prior you go on to fit with, with no second specification drifting out of sync behind it.

See also

• ?tulpa_priors, ?prior_normal, ?prior_pc and the other prior_* constructor pages.
• ?prior_predict and ?tulpa_family for prior predictive simulation.
• The getting-started vignette for fitting, prediction, and model comparison.
• Simpson, D., Rue, H., Riebler, A., Martins, T. G., & Sorbye, S. H. (2017). Penalising model component complexity. Statistical Science, 32(1), 1-28.