Nothing
Getting started with tulpa
In tulpa: Templated Unified Library for Posterior Approximation in Bayesian Hierarchical Models
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.width = 6, fig.height = 4
)
set.seed(20260529)
# The model fits below are evaluated when the article is built locally, in CI
# and for the pkgdown site (all of which set NOT_CRAN). CRAN's check farm gives
# the whole check a ten-minute budget, which these fits do not fit inside, so
# there the code is shown without being run.
EVAL_FITS <- identical(Sys.getenv("NOT_CRAN"), "true")
knitr::opts_chunk$set(eval = EVAL_FITS)
# Read here, not in a gated chunk: chunk options are evaluated whether or not
# the chunk runs, so a flag a chunk option reads has to exist either way.
has_ggplot <- requireNamespace("ggplot2", quietly = TRUE)
Installation
# Current release, built for Windows and macOS
install.packages("tulpa",
repos = c("https://gcol33.r-universe.dev",
"https://cloud.r-project.org"))
# Development version from GitHub
install.packages("pak")
pak::pak("gcol33/tulpa")
# Pin a release
pak::pak("gcol33/tulpa@v0.6.0")
pak resolves the dependency tree, including tulpaMesh (on CRAN, used for SPDE mesh construction). tulpa compiles its C++ backend on first install, so a C++17 toolchain is required: Rtools on Windows, Xcode CLI tools on macOS, r-base-dev on Linux.
library(tulpa)
What tulpa is
tulpa is a Bayesian hierarchical modelling engine. You write a formula such
as y ~ x + (1 | g), choose a family, and tulpa() parses the formula,
builds the design, picks an inference backend, and returns a fitted object
you can summarise, predict from, plot, and compare. The call shape is the
one you already know from lm() and glmer(); the machinery behind it is a
nested-approximation Bayesian engine rather than a frequentist optimiser or
a single MCMC sampler.
The engine carries the structural pieces that hierarchical models reuse:
fixed effects, random intercepts and slopes, spatial fields (SPDE, NNGP,
HSGP, ICAR, BYM2), temporal random walks, spatially and temporally varying
coefficients, and user-defined latent blocks. It also carries a tier system
for how the posterior is computed: a fast Laplace approximation, a
gradient sampler, an independence-MH correction, a Polya-Gamma Gibbs
sampler, or a Pathfinder variational fit, each with a stated epistemic
guarantee rather than a silent default. The mode that ran is recorded on the
fit, so the uncertainty you report can always be traced to the method that
produced it.
The design philosophy is nested approximation plus debias. A hierarchical
posterior splits into blocks. Some blocks are well behaved and a cheap
deterministic approximation handles them accurately: a Gaussian latent field
under enough data is close to Gaussian near its mode, and a Laplace step
nails it. Other directions are skewed or heavy-tailed, and there a short
burst of exact MCMC corrects the residual bias the approximation leaves
behind. The engine is built to compose those two ideas rather than choosing
between them once and for all. That framing places tulpa between two
familiar tools. INLA does the nested approximation but stops there, so it
inherits the approximation's bias on non-Gaussian residuals. Stan runs exact
MCMC on everything, which is correct but pays the full sampling price on
every block, including the ones a Laplace step would have settled
instantly. tulpa is the synthesis: approximate the easy blocks, sample only
the hard directions, and never hide which one ran where.
The observation likelihood is the one piece tulpa keeps deliberately small.
The engine ships five built-in families (gaussian, binomial, poisson,
neg_binomial_2, beta), which cover the common regression cases. Model
packages built on the engine plug their own likelihoods into the same
inference machinery through a LikelihoodSpec and inherit every tier, every
latent block, and every accessor for free. Occupancy and detection models
live in tulpaObs, ratio and rate models in tulpaRatio, generalized linear
mixed models in tulpaGlmm. This vignette stays on the engine's own surface:
simulate data with a known truth, fit it with tulpa(), pull the estimates
out, predict, compare specifications, and take a short tour of the latent
blocks. Each tour stop links to the dedicated vignette that goes deep.
One habit runs through everything below: simulate against a known truth.
Real data never reveal the values that generated them, so a fit that merely
looks reasonable proves nothing about whether the method works. Every worked
example here builds its data from coefficients you can see in the code, then
asks the fit to find them again. When the recovered intercept, slope, field,
or trend lands on the value that went in, the call did what it claimed. When
it does not, the gap is a bug you can catch before the same call meets data
where no truth is available to check against. The numbers in the printed
output are there to be compared against the numbers in the simulation, not
admired on their own.
A first model
The honest way to learn a fitting interface is to feed it data whose answer
you already know, then check that the fit finds it. Start with a Gaussian
response built from a fixed intercept of 0.5 and a slope of 1.2, with
residual noise at a standard deviation of 0.8. Four hundred rows is enough
for the estimates to settle close to the truth without being so many that
the fit is trivially easy.
n <- 400
x <- rnorm(n)
y <- 0.5 + 1.2 * x + rnorm(n, sd = 0.8)
df <- data.frame(y = y, x = x)
Fit with tulpa(). The formula and family mirror glm(); the new
argument is mode, which selects the inference backend. mode = "laplace"
is the Tier-2 Gaussian approximation: it finds the posterior mode and the
curvature there. For a Gaussian likelihood that curvature is exact, so the
fit returns the analytic posterior immediately, with no iteration and no
Monte Carlo error.
fit <- tulpa(y ~ x, data = df, family = "gaussian",
mode = "laplace", phi = 0.8^2)
coef(fit)
The coefficients land near the truth. coef() reports the posterior mean of
each fixed effect, the same point estimate you would read off a frequentist
fit, but here it is a posterior summary rather than a maximum-likelihood
point. The slope sits close to 1.2 and the intercept close to 0.5, which is
the first thing to check on any simulation: did the estimates recover the
values you put in.
A point estimate alone is half a fit. summary() adds the posterior
standard error and the 2.5% and 97.5% credible bounds for each fixed effect,
alongside how the fit was produced and how reliable it is, which is what turns "the slope is about 1.2" into "the slope is 1.2 with a
95% credible interval that excludes zero by a wide margin". Read the standard
error as the posterior spread of the coefficient and the bounds as the
interval the coefficient falls in with 95% posterior probability.
summary(fit)
The credible bounds bracket the true values, which is the second recovery
check. confint() returns the same bounds as a matrix, handy when you want
to program against them rather than read them, and vcov() returns the full
fixed-effect covariance, including the off-diagonal covariance between the
intercept and the slope. Both read the fixed-effect block of the Hessian the
Laplace fit already computed, so neither triggers an extra sampling step;
they are deterministic functions of a fit that has already run.
confint(fit)
The phi argument is the family's dispersion parameter, and its meaning
depends on the family. For gaussian it is the residual variance, for
neg_binomial_2 the size parameter that sets overdispersion, for beta the
precision. The Laplace and sampler paths both condition on the value you
pass: they treat phi as known and fit everything else given it. Here the
true residual standard deviation is 0.8, so the matching value is its
square, phi = 0.8^2. When you do not know phi in advance, the nested-Laplace and EM
layers integrate it rather than fixing it, which the inference vignettes
walk through; for a first fit, passing a plausible value is enough to get
estimates and intervals you can sanity-check.
That is the whole basic loop. Simulate or load data, call tulpa() with a
formula and a family, read coef() and summary() to check the estimates
and their uncertainty, and lean on confint() and vcov() when you need the
bounds or the covariance in a usable shape. Everything that follows builds on
this loop: random effects add a term to the formula, other families change
one argument, prediction and comparison are accessors on the object the loop
already returned. The interface stays flat as the models get richer.
It is worth pausing on what the returned object carries, because the rest of
the vignette pulls from it. The fit is a list with named fields. The fixed
effects and their covariance are there, the latent modes are there for fits
that have a latent block, the backend and tier that ran are there, and the
log marginal likelihood is there. The accessors are thin readers over those
fields: coef() reads the fixed effects, vcov() inverts the stored Hessian
block, logLik() reads the marginal likelihood. Nothing recomputes the fit,
so calling five accessors costs five field reads rather than five fits. That
matters when you compare a dozen specifications or bootstrap a prediction
band, because the expensive step ran once when you called tulpa() and every
summary after it is cheap.
Random effects
Grouped data break the independence a plain regression assumes. Measurements
taken on the same subject, plot, school, or site share something the
covariates do not capture, and ignoring that shared signal makes the
standard errors too small. A random intercept is the standard fix: one
offset per group, tied together by a common variance so groups borrow
strength from each other instead of each estimating its own level in
isolation. It enters the formula through the (1 | g) term, exactly as in
lme4.
Simulate twelve groups, each with its own offset drawn from a normal with
standard deviation 0.6, then add that offset to the linear predictor. The
fixed structure is the same intercept and slope as before, so the only new
ingredient is the group effect.
g <- factor(sample(1:12, n, replace = TRUE))
u <- rnorm(12, sd = 0.6)
df$y <- 0.5 + 1.2 * df$x + u[g] + rnorm(n, sd = 0.8)
df$g <- g
fit_re <- tulpa(y ~ x + (1 | g), data = df, family = "gaussian",
mode = "laplace", sigma_re = 0.6, phi = 0.8^2)
coef(fit_re)
coef() reports the fixed effects, which are unchanged in meaning: the
intercept and slope read against a background that now accounts for the group
structure. The group-level offsets themselves come from ranef(), which
returns one row per group carrying the posterior mean of that group's effect,
the quantity a mixed-model package calls the BLUP (best linear unbiased
predictor). These are the estimated departures of each group from the
population mean, shrunk toward zero by the random-effect variance.
head(ranef(fit_re), 4)
This fit conditions on sigma_re rather than estimating it: you pass the
random-effect standard deviation and the engine fits everything else given
that value. Passing one number applies it to the term; with several
random-effect terms you pass one value per term. Conditioning is a deliberate
design choice on the fast path, and it has a clean escape hatch. When you do
not have a value for sigma_re, fit at a couple of plausible ones (say 0.3
and 1.0) and watch whether the fixed effects move. If they hold steady, the
exact value barely matters and you report the fixed effects with confidence.
If they shift, the variance is informative and you should integrate it rather
than guess, which the nested-Laplace path does automatically.
The BLUPs that ranef() returns are shrunk estimates, and the shrinkage is
the point of pooling. A group with many observations gets a BLUP close to its
own raw mean, because the data in that group outweigh the prior pull toward
zero. A group with few observations gets a BLUP pulled hard toward the
population mean, because there is little group-specific signal to resist the
pull. That is the borrowing of strength a random effect buys over a fixed
effect per group: thin groups lean on the variance the whole sample informs,
rather than reporting a noisy mean from a handful of rows. The twelve groups
here each carry roughly thirty rows, so the shrinkage is mild and the BLUPs
sit close to the offsets that generated them, which is the recovery check for
the group level.
Random slopes are different in kind. A term such as (1 + x | g) lets the
slope itself vary by group, and the intercept and slope variances plus their
correlation form a covariance matrix rather than a single scalar. There is no
one sigma_re to condition on, so when a slope term is present the Laplace
path redirects on its own to a backend that integrates the whole covariance
matrix over a grid, with a PC-and-LKJ hyperprior and a Pareto-k-hat accuracy
diagnostic on the result. You reach that path not by naming a backend but by
writing the slope term and fitting at mode = "laplace". The inference-modes
vignette walks through a worked random-slope fit and the covariance it
returns.
Comparing models
A hierarchical model is rarely the only candidate. You usually have a ladder
of specifications, from an intercept-only baseline up through covariates and
structural terms, and the question is which rungs earn their place.
compare_models() ranks fits by a shared criterion and lays them out in one
table. With criterion = "loglik" it reports logLik() for each model,
which on a Laplace fit is the approximate log marginal likelihood: the
probability of the data with the latent quantities integrated out, the model
evidence. That is the right currency for comparing nested specifications,
because it already pays for any extra flexibility a richer model brings.
Build a three-rung ladder on the grouped data: an intercept-only model, a
model adding the slope, and a model adding the random intercept on top of the
slope. Each is a single tulpa() call differing only in the formula.
m0 <- tulpa(y ~ 1, data = df, family = "gaussian", mode = "laplace", phi = 0.8^2)
m1 <- tulpa(y ~ x, data = df, family = "gaussian", mode = "laplace", phi = 0.8^2)
m2 <- tulpa(y ~ x + (1 | g), data = df, family = "gaussian",
mode = "laplace", sigma_re = 0.6, phi = 0.8^2)
compare_models(intercept = m0, slope = m1, slope_re = m2, criterion = "loglik")
The table reports each model's name (taken straight from the argument
labels, so label them as a ladder), its fixed-effect parameter count, and its
log marginal likelihood. Adding the slope lifts the evidence sharply, because
the slope is real and large. Adding the random intercept lifts it further,
matching how the data were generated: the group offsets were genuinely there,
so the model that accounts for them is more probable. Evidence lives on the
log scale, so a difference of a few units is already a strong preference and a
difference in the tens is decisive.
Read the parameter count with one caveat. It counts fixed effects only, so
the random-intercept model shows the same count as the slope-only model even
though it carries the group structure. The random intercept is conditioned on
a fixed sigma_re rather than estimated as a fixed parameter, so the
bookkeeping does not increment the column, but the evidence still rewards the
structure because it integrates over the group offsets. Lean on the logLik
differences and treat the count as a fixed-effect tally rather than a model
complexity score. WAIC and LOO are available for fits that carry a pointwise
log-likelihood, which the model packages provide on top of their full
observation models; the engine's own comparison runs on the marginal
likelihood and the joint log-density, which is exactly what the ladder above
uses. The dedicated model-comparison vignette covers the WAIC and LOO paths
and how to read each criterion.
Prediction
A fitted model earns its keep when it predicts. predict() gives the
fixed-effect prediction at new covariate values you supply through
newdata. Random effects are held at their population mean of zero, so the
result is the marginal trend: what the covariate does on average across
groups, rather than what happens in any one group. se.fit = TRUE adds the
link-scale standard error and credible bounds, drawn from the same
fixed-effect covariance the summary reports, so the band on the prediction
inherits the uncertainty in the coefficients.
nd <- data.frame(x = seq(-2, 2, length.out = 50))
pr <- predict(fit_re, newdata = nd, se.fit = TRUE)
head(pr, 3)
The return value carries the prediction in fit and the bounds in lower
and upper, one row per row of newdata. Plotting the trend with its
credible band turns the table into the marginal effect you would put in a
paper: the line is the expected response across the covariate, and the
ribbon is the 95% credible region around it.
library(ggplot2)
ggplot(data.frame(x = nd$x, fit = pr$fit, lo = pr$lower, hi = pr$upper),
aes(x, fit)) +
geom_ribbon(aes(ymin = lo, ymax = hi), fill = "steelblue", alpha = 0.25) +
geom_line(linewidth = 1) +
labs(x = "x", y = "predicted y") +
theme(panel.background = element_rect(fill = "transparent"),
plot.background = element_rect(fill = "transparent"))
The band is narrow in the middle of the covariate range and flares at the
edges, which is the signature of a linear fit: the prediction is most certain
near the bulk of the data and least certain where you extrapolate. For a
non-Gaussian family, type = "response" maps the prediction through the
inverse link, and the credible bounds stay inside the response range, so a
binomial prediction reports probabilities in [0, 1] and a Poisson
prediction reports positive means. The default type = "link" keeps
everything on the linear-predictor scale, which is the natural scale for a
symmetric credible band; map to the response scale when you want the answer
in the units a reader expects.
Other families
The five built-in families share one call shape. Only family changes, plus
any dispersion phi the family needs. The same coef(), summary(),
predict(), and comparison machinery applies, because the engine treats the
likelihood as a pluggable component and everything above it is family-blind.
# Poisson counts
dp <- data.frame(y = rpois(n, exp(0.2 + 0.6 * x)), x = x)
coef(tulpa(y ~ x, data = dp, family = "poisson", mode = "laplace"))
# Beta proportions in (0, 1), precision phi
mu <- plogis(0.1 + 0.8 * x)
db <- data.frame(y = rbeta(n, mu * 8, (1 - mu) * 8), x = x)
coef(tulpa(y ~ x, data = db, family = "beta", mode = "laplace", phi = 8))
The Poisson slope recovers the 0.6 it was built from and the beta slope
recovers its 0.8, each on its own link scale: log for the Poisson rate, logit
for the beta mean. Counts go through poisson or, when they are
overdispersed, through neg_binomial_2 with phi setting the size. Binary
and binomial outcomes go through binomial, which the next sections use
heavily. Proportions strictly inside the unit interval go through beta with
phi as the precision. When your response does not fit one of these five,
the right move is a model package rather than a workaround in the formula:
zero-inflation, hurdles, detection and occupancy, and ratio responses all
have a dedicated likelihood in a package that plugs into this same engine.
Inference modes
Every fit so far ran at mode = "laplace", the Tier-2 Gaussian
approximation. The same model can be fit at a different tier, and the tier
matters because it pins down what the credible intervals mean. tulpa sorts
every backend into one of three tiers, and inference_mode_info() prints the
full map with the guarantee attached to each.
inference_mode_info()
Tier 1 (Exact) draws from the posterior, so its intervals are posterior
uncertainty up to Monte Carlo error that shrinks with more draws. MALA, the
independence-MH samplers, and the Polya-Gamma Gibbs sampler live here. Tier 2
(Structured) is accurate conditional on a stated assumption, a Gaussian
posterior shape for Laplace and Pathfinder, which holds for latent Gaussian
models with enough data and fails predictably when it does not. Tier 3
(Optimized) gives a point and a covariance from an optimisation objective with
no general guarantee on the spread, and tulpa reaches it only when you ask for
it by name. mode = "auto" chooses between Tier 1 and Tier 2 by a recorded
rule and never silently drops to Tier 3.
Fit the binomial grouped model two ways, once with the deterministic Laplace
approximation and once with MALA, a gradient sampler that draws from the
posterior directly.
df_b <- df
df_b$y <- rbinom(n, 1, plogis(-0.3 + 1.0 * df$x + u[df$g]))
fit_lap <- tulpa(y ~ x + (1 | g), data = df_b, family = "binomial",
mode = "laplace", sigma_re = 0.6)
fit_mala <- tulpa(y ~ x + (1 | g), data = df_b, family = "binomial",
mode = "mala", sigma_re = 0.6,
control = list(n_iter = 450, warmup = 150))
The two tiers agree on the fixed effects to within their uncertainty, which
is the expected outcome when the posterior is near-Gaussian: the posterior
mean is the easy quantity and every method finds it. The Laplace fit is
deterministic and immediate; the sampler carries Monte Carlo error and stays
agnostic about the posterior shape, which is the guarantee you pay for when
the posterior is not Gaussian.
data.frame(
term = names(coef(fit_lap)),
laplace = round(coef(fit_lap), 3),
mala = round(coef(fit_mala), 3)
)
The agreement here is itself a signal that Laplace was safe for this model.
Every fit records which backend ran and at which tier, and why mode =
"auto" chose what it did. print() names the backend and tier in its first
line, and summary() carries all three in its inference part, so you are
never left guessing which engine produced a number.
unlist(summary(fit_mala)$inference[c("backend", "inference_tier")])
The tiers encode an epistemic guarantee that goes beyond the speed
trade-off. Tier 1 intervals are posterior uncertainty you can quote without
an asterisk. Tier 2
intervals are correct conditional on the Gaussian shape, and the cheapest
check on that condition is mode = "imh_laplace": a high acceptance rate
says the posterior is close to the Laplace Gaussian and the intervals are
safe, a low one says it is far and Laplace was biased. The R-callable
backends, with the tier each belongs to and the situation each suits:
| Mode | Tier | What it does | Reach for it when |
|---|---|---|---|
| "laplace" | 2 | mode + curvature | a fast, deterministic fit; the default sanity check |
| "pathfinder" | 2 | variational draws along an optimisation path | you want draws but not full MCMC |
| "mala" | 1 | gradient-based Metropolis sampler | exact posterior moments, moderate dimension |
| "imh_laplace" | 1 | Laplace proposal with independence-MH correction | debiasing a Laplace fit at low parameter dimension |
| "gibbs" | 1 | Polya-Gamma sampler (binomial / negbin) | a conjugate-style exact fit for those families |
The trace plot of a sampler fit is one window onto whether the chain mixed.
plot(fit, type = "trace") draws the post-warmup draws for the leading fixed
effects; a healthy chain looks like a fuzzy horizontal band with no drift or
sticking.
plot(fit_mala, type = "trace")
The agreement between the tiers here is the typical case for a latent
Gaussian model with enough data, and it is also the case where the cheap
method is the right one. The interesting models are the ones where the
agreement breaks: a few dozen rows, a sparse binomial with most outcomes
zero, a variance component the data barely identify. There the posterior
turns skewed or heavy-tailed, the Gaussian approximation symmetrises it, and
the Tier-2 standard error drifts away from the Tier-1 one. The samplers track
the true spread through that drift; the Laplace fit keeps reporting the
symmetric interval its assumption forces. The moment to pay for a sampler is
the moment that drift would matter to your conclusion.
inference_mode_info() lists every backend, including the ones reachable only
from model packages through the C interface. The inference-modes vignette
fits one model at Laplace, MALA, and Pathfinder side by side, times each, and
shows exactly where the cost of the exact guarantee turns into a reason to
pay it.
Spatial fields
Areal data carry a neighbourhood structure: regions that sit next to each
other tend to look alike, and a regression that treats them as independent
draws gets the standard errors wrong and lets a spatially patterned covariate
absorb the regional signal. A spatial field is the place to put that signal.
You add one random effect per region, tie neighbouring effects together so
the field stays smooth, and read the covariate against a background that
already accounts for location.
Pass the adjacency through spatial= and name the per-observation unit with
a spatial(col) term in the formula. type = "icar" puts an intrinsic
conditional-autoregressive field on the units. Here twenty regions sit on a
ring, each adjacent to its two neighbours, with a smooth wave of regional
effect around the loop.
K <- 20
W <- matrix(0, K, K)
for (i in 1:K) { j <- if (i < K) i + 1 else 1; W[i, j] <- W[j, i] <- 1 }
region <- factor(sample(1:K, n, replace = TRUE))
field <- as.numeric(scale(sin(2 * pi * (1:K) / K)))[region]
ds <- data.frame(y = rbinom(n, 1, plogis(-0.2 + 0.7 * x + field)),
x = x, region = region)
fit_sp <- tulpa(y ~ x + spatial(region), data = ds, family = "binomial",
spatial = list(type = "icar", adjacency = W),
mode = "laplace")
coef(fit_sp)
The slope recovers near its true 0.7 while the ICAR field absorbs the
regional structure. Plot the per-region field effect around the ring to see
the wave the fit pulled out of binary data: field_effects() returns the
field at each region, here its posterior mode, and traces the smooth signal
the covariate could not explain.
library(ggplot2)
fe <- field_effects(fit_sp)
ggplot(data.frame(region = fe$unit, effect = fe$estimate), aes(region, effect)) +
geom_line(linewidth = 1, colour = "steelblue") +
geom_point() +
labs(x = "region", y = "field effect") +
theme(panel.background = element_rect(fill = "transparent"),
plot.background = element_rect(fill = "transparent"))
The field absorbs the regional signal so the slope does not have to. Without
the field, a covariate that happens to vary across the ring the same way the
latent wave does would soak up part of that wave and report a biased slope.
With the field carrying the spatial structure, the covariate reports its own
effect, and the standard error widens to its honest size because the field,
not the slope, now carries the autocorrelation. The intrinsic prior fixes the
field only up to an additive constant, so the intercept floats against a
near-flat prior and its standard error runs wide by design; read the slopes,
which are identified, and treat the intercept as a nuisance level.
Continuous fields use a coordinate spec instead of an adjacency:
spatial_gp(~ lon + lat) for a nearest-neighbour GP,
spatial_gp(~ lon + lat, approx = "hsgp") for a Hilbert-space GP,
spatial_spde(~ lon + lat, data) for a Matern
SPDE solved on a mesh through fit_spde(). The spatial vignette walks
through the areal path end to end, recovers a known field, integrates the
smoothing precision rather than fixing it, and shows how the field's
uncertainty flows into the fixed-effect standard errors.
Temporal fields
Datasets that carry a clock carry the temporal twin of spatial confounding.
A smooth background that drifts over time will be read as covariate signal if
the covariate trends the same way, and the slope comes out too steep. A
flexible function of time soaks up the background and leaves the slope to
explain only what wiggles faster than the trend. tulpa offers that function
as a first-order random walk, declared with temporal_rw1() and passed
through the temporal= argument.
Simulate counts over thirty time points with a smooth sinusoidal trend, a
real covariate effect, and Poisson noise, then fit with the random-walk trend
in place.
Tt <- 30
time <- sample(seq_len(Tt), n, replace = TRUE)
trend <- 1.2 * sin(2 * pi * seq_len(Tt) / Tt)
xt <- rnorm(n)
dt <- data.frame(y = rpois(n, exp(0.4 + 0.5 * xt + trend[time])),
x = xt, time = time)
fit_t <- tulpa(y ~ x, data = dt, family = "poisson",
temporal = temporal_rw1("time"), mode = "auto")
coef(fit_t)
The slope lands near its true 0.5 with the trend absorbed into the field
rather than the coefficient. mode = "auto" routes the temporal field to the
nested-Laplace backend, and tulpa_outer_grid() returns its grid: the
per-grid latent modes in modes (one row per grid point, columns ordered
[fixed effects, field]) and the grid weights. Weight-average the field
column for each time point to integrate over the smoothing precision, then
plot the recovered walk against the truth.
g <- tulpa_outer_grid(fit_t)
te <- colSums(g$weights * g$modes)[g$n_fixed + seq_len(Tt)]
ggplot(data.frame(time = seq_len(Tt),
est = te - mean(te),
truth = trend - mean(trend)),
aes(time)) +
geom_line(aes(y = truth), linewidth = 1, colour = "grey50") +
geom_line(aes(y = est), linewidth = 1, colour = "steelblue") +
labs(x = "time", y = "trend (centred)") +
theme(panel.background = element_rect(fill = "transparent"),
plot.background = element_rect(fill = "transparent"))
The estimated walk follows the wave, centred to remove the level the
sum-to-zero handling leaves free. The front door carries temporal_rw1()
today; the second-order walk, the AR(1) process, and per-group panel trends
have direct fitters the temporal vignette covers, which walks the full path
from model to recovery to honest uncertainty.
Priors
Every Bayesian fit carries priors whether you write them down or not. tulpa
keeps them in one tulpa_priors object that names a prior for each kind of
parameter the engine knows about, built with tulpa_priors() and filled with
prior_*() constructors. The same object drives prior predictive simulation,
so the priors you check are the priors you fit with. Three jobs run through
the prior system: a fixed effect lives on the link scale and gets a prior over
the whole real line, a variance is positive and gets a prior on the
half-line, a proportion is trapped in the unit interval.
Look at what a fixed-effect prior implies before any data touch it.
prior_predict() draws coefficients from their priors, pushes them through
the linear predictor, and simulates responses, so you can see whether the
prior expects data the response could never produce.
fam <- tulpa_family(
"poisson",
function(eta, params, n_obs, ...) rpois(n_obs, exp(eta[[1]]))
)
pp <- prior_predict(y ~ x, family = fam, data = dp, n_draws = 50,
priors = tulpa_priors(beta = prior_normal(0, 1)))
pp
pp$y holds one simulated dataset per prior draw. Overlay their densities to
see the range the prior expects, then tighten the prior if it covers ranges
the response could never take.
sims <- data.frame(
value = unlist(pp$y),
draw = factor(rep(seq_along(pp$y), lengths(pp$y)))
)
ggplot(sims, aes(value, group = draw)) +
geom_density(colour = "steelblue", alpha = 0.3) +
labs(x = "simulated y", y = "density") +
theme(panel.background = element_rect(fill = "transparent"),
plot.background = element_rect(fill = "transparent"))
A wide prior predictive that spills past the plausible range says the priors
are too vague; tightening prior_normal(0, 2.5) toward prior_normal(0, 1)
pulls it in. The priors vignette covers all seven constructors, the
penalised-complexity prior prior_pc() for scale parameters, the defaults,
and how to move a prior and watch the posterior follow.
Custom latent blocks with tgmrf()
When the latent structure you need is not a built-in field, tgmrf() lets a
script define its own Gaussian Markov random field and plug it in as a
latent(...) term that reaches every inference tier. You supply two closures
and an init vector: Q(theta) returns the sparse precision matrix at a given
hyperparameter, prior(theta) returns its log-density, and init sets the
starting hyperparameter. No autodiff, no DSL, no codegen; the engine calls
the closures numerically and owns everything from the Newton system up.
Define a first-order random-walk precision on a chain of twenty nodes, the
same structure the temporal path builds internally, but written by hand to
show the contract.
library(Matrix)
m <- 20
D2 <- diff(diag(m), differences = 1)
R <- crossprod(D2) # RW1 structure matrix
blk <- tgmrf(
Q = function(theta) as(theta[1] * (R + 1e-4 * diag(m)), "dgCMatrix"),
prior = function(theta) dgamma(theta[1], 2, 1, log = TRUE),
init = c(tau = 1)
)
blk
The constructor evaluates Q(init) once to catch errors early, infer the
field length, and capture the sparsity pattern the inner solver reuses at
every grid point. The returned block carries the precision family, the
hyperparameter names, and the registry key the nested-Laplace dispatcher
recognises, and prints a compact summary of what it found.
Q1 <- as.matrix(blk$Q(c(tau = 1)))
ix <- which(Q1 != 0, arr.ind = TRUE)
ggplot(data.frame(row = ix[, 1], col = ix[, 2]), aes(col, row)) +
geom_tile(fill = "steelblue") +
scale_y_reverse() +
labs(x = "column", y = "row") +
theme(panel.background = element_rect(fill = "transparent"),
plot.background = element_rect(fill = "transparent"))
The tridiagonal pattern is the Markov property written as a matrix: each node
is conditionally independent of every non-neighbour given its neighbours, and
that independence is exactly a zero off the three central diagonals. The
tgmrf() vignette builds a periodic AR(1) block end to end, fits it through
the full tier ladder from nested Laplace to exact sampling, and shows how one
block written once reaches all four tiers without modification.
S3 methods
The fitted object answers the standard accessors, and they work the same way
on a Laplace fit and a sampler fit. The Laplace tier reports the Gaussian
approximation it computed; the sampler tier reports empirical posterior
summaries from its draws. Reaching for the accessor you already know from
base R or broom is the intended workflow.
| Accessor | Returns |
|---|---|
| coef() | named fixed-effect estimates |
| summary() | the full report: backend and tier, fixed effects with SE and credible bounds, variance components, hyperparameters, reliability; coef(summary(fit)) is the fixed-effect table |
| confint() | credible intervals (matrix) |
| vcov() | fixed-effect covariance |
| ranef() | group-level random effects |
| tidy() | one-row-per-term data frame (broom style) |
| glance() | one-row model summary |
| diagnostics() | convergence or approximation reliability, per parameter |
| tulpa_outer_grid() | a nested-Laplace fit's hyperparameter grid: nodes, weights, per-cell modes |
| predict(), fitted() | fixed-effect prediction / in-sample mean |
| logLik() | log marginal likelihood (Laplace) or mean log-density |
| plot() | fixed-effect densities (or trace for sampler fits) |
tidy() returns the fixed-effect table in the broom convention: one row per
term, with the estimate, standard error, and credible bounds in named
columns, ready to feed a table or a coefficient plot.
tidy(fit_re)
glance() returns a single-row model summary carrying the fixed-effect
count, the number of draws, the log-likelihood, the divergence count, the
mean acceptance rate, and the convergence flag. On a Laplace fit the
sampler-only fields read NA, which is the honest report; on a sampler fit
they carry the diagnostics that say whether the chain mixed.
glance(fit_mala)[c("n_samples", "logLik", "mean_accept", "n_divergent")]
plot() defaults to the fixed-effect densities, a quick visual of where each
coefficient sits and how wide its posterior is. On a sampler fit type =
"trace" shows the chains instead, as in the inference section above, and
type = "pairs" shows the joint posterior of the leading parameters.
coef(), confint(), vcov(), and ranef() round out the set, each
reading off the fit that has already run rather than triggering new
computation. Because the model packages inherit class = c("model_fit",
"tulpa_fit"), every accessor here works unchanged on an occupancy fit, a
ratio fit, or a GLMM fit, with the package adding only its own model-specific
methods on top.
Practical guidance
A few rules of thumb for everyday use:
- Start at Tier 2.
mode = "laplace" is deterministic and returns in
well under a second on tens of thousands of rows. Use it for the first
look and for model comparison by marginal likelihood. Move to a sampler
(mode = "mala") when you need exact posterior moments or the posterior
is visibly non-Gaussian, and check first with mode = "imh_laplace",
whose acceptance rate is a direct verdict on the Laplace approximation.
sigma_re is conditioned, not estimated, on the fixed-sigma_re
paths. With a random intercept the fit assumes the value you pass. If
you do not have one, fit at a couple of values (say 0.3 and 1.0) and
check whether the fixed effects move; if they do, integrate the variance
through the nested-Laplace path rather than guessing. Random slopes route
there on their own.- Sample size. A random-effect term needs enough groups to inform its
variance: with fewer than about 5 to 8 groups, the group-level estimates
are dominated by the prior, and a fixed effect for the grouping is often
the honest choice. Spatial fields want dozens to thousands of regions,
many of them thin, to earn their machinery.
- Watch convergence on the sampler tiers.
glance() reports
mean_accept and the divergence count; an acceptance rate far from
0.4 to 0.8, or any divergences, means the chain needs more warmup or a
smaller step size (control$epsilon). A Tier-1 fit is only nominally
exact until the chain has been shown to mix.
- Reach for a model package when the likelihood is not built in. The
five families here cover the common cases. Zero-inflation, hurdles,
detection and occupancy, ratio responses, and posterior-predictive
checking live in the packages built on this engine, not in
tulpa()
itself, and they inherit every tier and accessor shown here.
Where to go next
- Inference details and the tier guarantees: the inference-modes and
algorithm vignettes.
- Spatial and temporal fields, spatially varying coefficients: the
spatial, temporal, and SVC vignettes.
- Priors and prior predictive checks: the priors vignette and
prior_predict().
-
User-defined latent structure: the tgmrf() vignette.
-
Model comparison by marginal likelihood, WAIC, and LOO: the
model-comparison vignette.
- Posterior-predictive checks, WAIC, and LOO on a full observation model:
the model packages (tulpaObs, tulpaRatio, tulpaGlmm) that build on this
engine.
```
Try the tulpa package in your browser
Any scripts or data that you put into this service are public.
tulpa documentation built on Oct. 19, 2026, 5:09 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 4 ) set.seed(20260529)
# The model fits below are evaluated when the article is built locally, in CI # and for the pkgdown site (all of which set NOT_CRAN). CRAN's check farm gives # the whole check a ten-minute budget, which these fits do not fit inside, so # there the code is shown without being run. EVAL_FITS <- identical(Sys.getenv("NOT_CRAN"), "true") knitr::opts_chunk$set(eval = EVAL_FITS) # Read here, not in a gated chunk: chunk options are evaluated whether or not # the chunk runs, so a flag a chunk option reads has to exist either way. has_ggplot <- requireNamespace("ggplot2", quietly = TRUE)
Installation
# Current release, built for Windows and macOS install.packages("tulpa", repos = c("https://gcol33.r-universe.dev", "https://cloud.r-project.org")) # Development version from GitHub install.packages("pak") pak::pak("gcol33/tulpa") # Pin a release pak::pak("gcol33/tulpa@v0.6.0")
pak resolves the dependency tree, including tulpaMesh (on CRAN, used for SPDE mesh construction). tulpa compiles its C++ backend on first install, so a C++17 toolchain is required: Rtools on Windows, Xcode CLI tools on macOS, r-base-dev on Linux.
library(tulpa)
What tulpa is
tulpa is a Bayesian hierarchical modelling engine. You write a formula such
as y ~ x + (1 | g), choose a family, and tulpa() parses the formula,
builds the design, picks an inference backend, and returns a fitted object
you can summarise, predict from, plot, and compare. The call shape is the
one you already know from lm() and glmer(); the machinery behind it is a
nested-approximation Bayesian engine rather than a frequentist optimiser or
a single MCMC sampler.
The engine carries the structural pieces that hierarchical models reuse: fixed effects, random intercepts and slopes, spatial fields (SPDE, NNGP, HSGP, ICAR, BYM2), temporal random walks, spatially and temporally varying coefficients, and user-defined latent blocks. It also carries a tier system for how the posterior is computed: a fast Laplace approximation, a gradient sampler, an independence-MH correction, a Polya-Gamma Gibbs sampler, or a Pathfinder variational fit, each with a stated epistemic guarantee rather than a silent default. The mode that ran is recorded on the fit, so the uncertainty you report can always be traced to the method that produced it.
The design philosophy is nested approximation plus debias. A hierarchical posterior splits into blocks. Some blocks are well behaved and a cheap deterministic approximation handles them accurately: a Gaussian latent field under enough data is close to Gaussian near its mode, and a Laplace step nails it. Other directions are skewed or heavy-tailed, and there a short burst of exact MCMC corrects the residual bias the approximation leaves behind. The engine is built to compose those two ideas rather than choosing between them once and for all. That framing places tulpa between two familiar tools. INLA does the nested approximation but stops there, so it inherits the approximation's bias on non-Gaussian residuals. Stan runs exact MCMC on everything, which is correct but pays the full sampling price on every block, including the ones a Laplace step would have settled instantly. tulpa is the synthesis: approximate the easy blocks, sample only the hard directions, and never hide which one ran where.
The observation likelihood is the one piece tulpa keeps deliberately small.
The engine ships five built-in families (gaussian, binomial, poisson,
neg_binomial_2, beta), which cover the common regression cases. Model
packages built on the engine plug their own likelihoods into the same
inference machinery through a LikelihoodSpec and inherit every tier, every
latent block, and every accessor for free. Occupancy and detection models
live in tulpaObs, ratio and rate models in tulpaRatio, generalized linear
mixed models in tulpaGlmm. This vignette stays on the engine's own surface:
simulate data with a known truth, fit it with tulpa(), pull the estimates
out, predict, compare specifications, and take a short tour of the latent
blocks. Each tour stop links to the dedicated vignette that goes deep.
One habit runs through everything below: simulate against a known truth. Real data never reveal the values that generated them, so a fit that merely looks reasonable proves nothing about whether the method works. Every worked example here builds its data from coefficients you can see in the code, then asks the fit to find them again. When the recovered intercept, slope, field, or trend lands on the value that went in, the call did what it claimed. When it does not, the gap is a bug you can catch before the same call meets data where no truth is available to check against. The numbers in the printed output are there to be compared against the numbers in the simulation, not admired on their own.
A first model
The honest way to learn a fitting interface is to feed it data whose answer you already know, then check that the fit finds it. Start with a Gaussian response built from a fixed intercept of 0.5 and a slope of 1.2, with residual noise at a standard deviation of 0.8. Four hundred rows is enough for the estimates to settle close to the truth without being so many that the fit is trivially easy.
n <- 400 x <- rnorm(n) y <- 0.5 + 1.2 * x + rnorm(n, sd = 0.8) df <- data.frame(y = y, x = x)
Fit with tulpa(). The formula and family mirror glm(); the new
argument is mode, which selects the inference backend. mode = "laplace"
is the Tier-2 Gaussian approximation: it finds the posterior mode and the
curvature there. For a Gaussian likelihood that curvature is exact, so the
fit returns the analytic posterior immediately, with no iteration and no
Monte Carlo error.
fit <- tulpa(y ~ x, data = df, family = "gaussian", mode = "laplace", phi = 0.8^2) coef(fit)
The coefficients land near the truth. coef() reports the posterior mean of
each fixed effect, the same point estimate you would read off a frequentist
fit, but here it is a posterior summary rather than a maximum-likelihood
point. The slope sits close to 1.2 and the intercept close to 0.5, which is
the first thing to check on any simulation: did the estimates recover the
values you put in.
A point estimate alone is half a fit. summary() adds the posterior
standard error and the 2.5% and 97.5% credible bounds for each fixed effect,
alongside how the fit was produced and how reliable it is, which is what turns "the slope is about 1.2" into "the slope is 1.2 with a
95% credible interval that excludes zero by a wide margin". Read the standard
error as the posterior spread of the coefficient and the bounds as the
interval the coefficient falls in with 95% posterior probability.
summary(fit)
The credible bounds bracket the true values, which is the second recovery
check. confint() returns the same bounds as a matrix, handy when you want
to program against them rather than read them, and vcov() returns the full
fixed-effect covariance, including the off-diagonal covariance between the
intercept and the slope. Both read the fixed-effect block of the Hessian the
Laplace fit already computed, so neither triggers an extra sampling step;
they are deterministic functions of a fit that has already run.
confint(fit)
The phi argument is the family's dispersion parameter, and its meaning
depends on the family. For gaussian it is the residual variance, for
neg_binomial_2 the size parameter that sets overdispersion, for beta the
precision. The Laplace and sampler paths both condition on the value you
pass: they treat phi as known and fit everything else given it. Here the
true residual standard deviation is 0.8, so the matching value is its
square, phi = 0.8^2. When you do not know phi in advance, the nested-Laplace and EM
layers integrate it rather than fixing it, which the inference vignettes
walk through; for a first fit, passing a plausible value is enough to get
estimates and intervals you can sanity-check.
That is the whole basic loop. Simulate or load data, call tulpa() with a
formula and a family, read coef() and summary() to check the estimates
and their uncertainty, and lean on confint() and vcov() when you need the
bounds or the covariance in a usable shape. Everything that follows builds on
this loop: random effects add a term to the formula, other families change
one argument, prediction and comparison are accessors on the object the loop
already returned. The interface stays flat as the models get richer.
It is worth pausing on what the returned object carries, because the rest of
the vignette pulls from it. The fit is a list with named fields. The fixed
effects and their covariance are there, the latent modes are there for fits
that have a latent block, the backend and tier that ran are there, and the
log marginal likelihood is there. The accessors are thin readers over those
fields: coef() reads the fixed effects, vcov() inverts the stored Hessian
block, logLik() reads the marginal likelihood. Nothing recomputes the fit,
so calling five accessors costs five field reads rather than five fits. That
matters when you compare a dozen specifications or bootstrap a prediction
band, because the expensive step ran once when you called tulpa() and every
summary after it is cheap.
Random effects
Grouped data break the independence a plain regression assumes. Measurements
taken on the same subject, plot, school, or site share something the
covariates do not capture, and ignoring that shared signal makes the
standard errors too small. A random intercept is the standard fix: one
offset per group, tied together by a common variance so groups borrow
strength from each other instead of each estimating its own level in
isolation. It enters the formula through the (1 | g) term, exactly as in
lme4.
Simulate twelve groups, each with its own offset drawn from a normal with standard deviation 0.6, then add that offset to the linear predictor. The fixed structure is the same intercept and slope as before, so the only new ingredient is the group effect.
g <- factor(sample(1:12, n, replace = TRUE)) u <- rnorm(12, sd = 0.6) df$y <- 0.5 + 1.2 * df$x + u[g] + rnorm(n, sd = 0.8) df$g <- g fit_re <- tulpa(y ~ x + (1 | g), data = df, family = "gaussian", mode = "laplace", sigma_re = 0.6, phi = 0.8^2) coef(fit_re)
coef() reports the fixed effects, which are unchanged in meaning: the
intercept and slope read against a background that now accounts for the group
structure. The group-level offsets themselves come from ranef(), which
returns one row per group carrying the posterior mean of that group's effect,
the quantity a mixed-model package calls the BLUP (best linear unbiased
predictor). These are the estimated departures of each group from the
population mean, shrunk toward zero by the random-effect variance.
head(ranef(fit_re), 4)
This fit conditions on sigma_re rather than estimating it: you pass the
random-effect standard deviation and the engine fits everything else given
that value. Passing one number applies it to the term; with several
random-effect terms you pass one value per term. Conditioning is a deliberate
design choice on the fast path, and it has a clean escape hatch. When you do
not have a value for sigma_re, fit at a couple of plausible ones (say 0.3
and 1.0) and watch whether the fixed effects move. If they hold steady, the
exact value barely matters and you report the fixed effects with confidence.
If they shift, the variance is informative and you should integrate it rather
than guess, which the nested-Laplace path does automatically.
The BLUPs that ranef() returns are shrunk estimates, and the shrinkage is
the point of pooling. A group with many observations gets a BLUP close to its
own raw mean, because the data in that group outweigh the prior pull toward
zero. A group with few observations gets a BLUP pulled hard toward the
population mean, because there is little group-specific signal to resist the
pull. That is the borrowing of strength a random effect buys over a fixed
effect per group: thin groups lean on the variance the whole sample informs,
rather than reporting a noisy mean from a handful of rows. The twelve groups
here each carry roughly thirty rows, so the shrinkage is mild and the BLUPs
sit close to the offsets that generated them, which is the recovery check for
the group level.
Random slopes are different in kind. A term such as (1 + x | g) lets the
slope itself vary by group, and the intercept and slope variances plus their
correlation form a covariance matrix rather than a single scalar. There is no
one sigma_re to condition on, so when a slope term is present the Laplace
path redirects on its own to a backend that integrates the whole covariance
matrix over a grid, with a PC-and-LKJ hyperprior and a Pareto-k-hat accuracy
diagnostic on the result. You reach that path not by naming a backend but by
writing the slope term and fitting at mode = "laplace". The inference-modes
vignette walks through a worked random-slope fit and the covariance it
returns.
Comparing models
A hierarchical model is rarely the only candidate. You usually have a ladder
of specifications, from an intercept-only baseline up through covariates and
structural terms, and the question is which rungs earn their place.
compare_models() ranks fits by a shared criterion and lays them out in one
table. With criterion = "loglik" it reports logLik() for each model,
which on a Laplace fit is the approximate log marginal likelihood: the
probability of the data with the latent quantities integrated out, the model
evidence. That is the right currency for comparing nested specifications,
because it already pays for any extra flexibility a richer model brings.
Build a three-rung ladder on the grouped data: an intercept-only model, a
model adding the slope, and a model adding the random intercept on top of the
slope. Each is a single tulpa() call differing only in the formula.
m0 <- tulpa(y ~ 1, data = df, family = "gaussian", mode = "laplace", phi = 0.8^2) m1 <- tulpa(y ~ x, data = df, family = "gaussian", mode = "laplace", phi = 0.8^2) m2 <- tulpa(y ~ x + (1 | g), data = df, family = "gaussian", mode = "laplace", sigma_re = 0.6, phi = 0.8^2) compare_models(intercept = m0, slope = m1, slope_re = m2, criterion = "loglik")
The table reports each model's name (taken straight from the argument labels, so label them as a ladder), its fixed-effect parameter count, and its log marginal likelihood. Adding the slope lifts the evidence sharply, because the slope is real and large. Adding the random intercept lifts it further, matching how the data were generated: the group offsets were genuinely there, so the model that accounts for them is more probable. Evidence lives on the log scale, so a difference of a few units is already a strong preference and a difference in the tens is decisive.
Read the parameter count with one caveat. It counts fixed effects only, so
the random-intercept model shows the same count as the slope-only model even
though it carries the group structure. The random intercept is conditioned on
a fixed sigma_re rather than estimated as a fixed parameter, so the
bookkeeping does not increment the column, but the evidence still rewards the
structure because it integrates over the group offsets. Lean on the logLik
differences and treat the count as a fixed-effect tally rather than a model
complexity score. WAIC and LOO are available for fits that carry a pointwise
log-likelihood, which the model packages provide on top of their full
observation models; the engine's own comparison runs on the marginal
likelihood and the joint log-density, which is exactly what the ladder above
uses. The dedicated model-comparison vignette covers the WAIC and LOO paths
and how to read each criterion.
Prediction
A fitted model earns its keep when it predicts. predict() gives the
fixed-effect prediction at new covariate values you supply through
newdata. Random effects are held at their population mean of zero, so the
result is the marginal trend: what the covariate does on average across
groups, rather than what happens in any one group. se.fit = TRUE adds the
link-scale standard error and credible bounds, drawn from the same
fixed-effect covariance the summary reports, so the band on the prediction
inherits the uncertainty in the coefficients.
nd <- data.frame(x = seq(-2, 2, length.out = 50)) pr <- predict(fit_re, newdata = nd, se.fit = TRUE) head(pr, 3)
The return value carries the prediction in fit and the bounds in lower
and upper, one row per row of newdata. Plotting the trend with its
credible band turns the table into the marginal effect you would put in a
paper: the line is the expected response across the covariate, and the
ribbon is the 95% credible region around it.
library(ggplot2) ggplot(data.frame(x = nd$x, fit = pr$fit, lo = pr$lower, hi = pr$upper), aes(x, fit)) + geom_ribbon(aes(ymin = lo, ymax = hi), fill = "steelblue", alpha = 0.25) + geom_line(linewidth = 1) + labs(x = "x", y = "predicted y") + theme(panel.background = element_rect(fill = "transparent"), plot.background = element_rect(fill = "transparent"))
The band is narrow in the middle of the covariate range and flares at the
edges, which is the signature of a linear fit: the prediction is most certain
near the bulk of the data and least certain where you extrapolate. For a
non-Gaussian family, type = "response" maps the prediction through the
inverse link, and the credible bounds stay inside the response range, so a
binomial prediction reports probabilities in [0, 1] and a Poisson
prediction reports positive means. The default type = "link" keeps
everything on the linear-predictor scale, which is the natural scale for a
symmetric credible band; map to the response scale when you want the answer
in the units a reader expects.
Other families
The five built-in families share one call shape. Only family changes, plus
any dispersion phi the family needs. The same coef(), summary(),
predict(), and comparison machinery applies, because the engine treats the
likelihood as a pluggable component and everything above it is family-blind.
# Poisson counts dp <- data.frame(y = rpois(n, exp(0.2 + 0.6 * x)), x = x) coef(tulpa(y ~ x, data = dp, family = "poisson", mode = "laplace")) # Beta proportions in (0, 1), precision phi mu <- plogis(0.1 + 0.8 * x) db <- data.frame(y = rbeta(n, mu * 8, (1 - mu) * 8), x = x) coef(tulpa(y ~ x, data = db, family = "beta", mode = "laplace", phi = 8))
The Poisson slope recovers the 0.6 it was built from and the beta slope
recovers its 0.8, each on its own link scale: log for the Poisson rate, logit
for the beta mean. Counts go through poisson or, when they are
overdispersed, through neg_binomial_2 with phi setting the size. Binary
and binomial outcomes go through binomial, which the next sections use
heavily. Proportions strictly inside the unit interval go through beta with
phi as the precision. When your response does not fit one of these five,
the right move is a model package rather than a workaround in the formula:
zero-inflation, hurdles, detection and occupancy, and ratio responses all
have a dedicated likelihood in a package that plugs into this same engine.
Inference modes
Every fit so far ran at mode = "laplace", the Tier-2 Gaussian
approximation. The same model can be fit at a different tier, and the tier
matters because it pins down what the credible intervals mean. tulpa sorts
every backend into one of three tiers, and inference_mode_info() prints the
full map with the guarantee attached to each.
inference_mode_info()
Tier 1 (Exact) draws from the posterior, so its intervals are posterior
uncertainty up to Monte Carlo error that shrinks with more draws. MALA, the
independence-MH samplers, and the Polya-Gamma Gibbs sampler live here. Tier 2
(Structured) is accurate conditional on a stated assumption, a Gaussian
posterior shape for Laplace and Pathfinder, which holds for latent Gaussian
models with enough data and fails predictably when it does not. Tier 3
(Optimized) gives a point and a covariance from an optimisation objective with
no general guarantee on the spread, and tulpa reaches it only when you ask for
it by name. mode = "auto" chooses between Tier 1 and Tier 2 by a recorded
rule and never silently drops to Tier 3.
Fit the binomial grouped model two ways, once with the deterministic Laplace approximation and once with MALA, a gradient sampler that draws from the posterior directly.
df_b <- df df_b$y <- rbinom(n, 1, plogis(-0.3 + 1.0 * df$x + u[df$g])) fit_lap <- tulpa(y ~ x + (1 | g), data = df_b, family = "binomial", mode = "laplace", sigma_re = 0.6) fit_mala <- tulpa(y ~ x + (1 | g), data = df_b, family = "binomial", mode = "mala", sigma_re = 0.6, control = list(n_iter = 450, warmup = 150))
The two tiers agree on the fixed effects to within their uncertainty, which is the expected outcome when the posterior is near-Gaussian: the posterior mean is the easy quantity and every method finds it. The Laplace fit is deterministic and immediate; the sampler carries Monte Carlo error and stays agnostic about the posterior shape, which is the guarantee you pay for when the posterior is not Gaussian.
data.frame( term = names(coef(fit_lap)), laplace = round(coef(fit_lap), 3), mala = round(coef(fit_mala), 3) )
The agreement here is itself a signal that Laplace was safe for this model.
Every fit records which backend ran and at which tier, and why mode =
"auto" chose what it did. print() names the backend and tier in its first
line, and summary() carries all three in its inference part, so you are
never left guessing which engine produced a number.
unlist(summary(fit_mala)$inference[c("backend", "inference_tier")])
The tiers encode an epistemic guarantee that goes beyond the speed
trade-off. Tier 1 intervals are posterior uncertainty you can quote without
an asterisk. Tier 2
intervals are correct conditional on the Gaussian shape, and the cheapest
check on that condition is mode = "imh_laplace": a high acceptance rate
says the posterior is close to the Laplace Gaussian and the intervals are
safe, a low one says it is far and Laplace was biased. The R-callable
backends, with the tier each belongs to and the situation each suits:
| Mode | Tier | What it does | Reach for it when |
|---|---|---|---|
| "laplace" | 2 | mode + curvature | a fast, deterministic fit; the default sanity check |
| "pathfinder" | 2 | variational draws along an optimisation path | you want draws but not full MCMC |
| "mala" | 1 | gradient-based Metropolis sampler | exact posterior moments, moderate dimension |
| "imh_laplace" | 1 | Laplace proposal with independence-MH correction | debiasing a Laplace fit at low parameter dimension |
| "gibbs" | 1 | Polya-Gamma sampler (binomial / negbin) | a conjugate-style exact fit for those families |
The trace plot of a sampler fit is one window onto whether the chain mixed.
plot(fit, type = "trace") draws the post-warmup draws for the leading fixed
effects; a healthy chain looks like a fuzzy horizontal band with no drift or
sticking.
plot(fit_mala, type = "trace")
The agreement between the tiers here is the typical case for a latent Gaussian model with enough data, and it is also the case where the cheap method is the right one. The interesting models are the ones where the agreement breaks: a few dozen rows, a sparse binomial with most outcomes zero, a variance component the data barely identify. There the posterior turns skewed or heavy-tailed, the Gaussian approximation symmetrises it, and the Tier-2 standard error drifts away from the Tier-1 one. The samplers track the true spread through that drift; the Laplace fit keeps reporting the symmetric interval its assumption forces. The moment to pay for a sampler is the moment that drift would matter to your conclusion.
inference_mode_info() lists every backend, including the ones reachable only
from model packages through the C interface. The inference-modes vignette
fits one model at Laplace, MALA, and Pathfinder side by side, times each, and
shows exactly where the cost of the exact guarantee turns into a reason to
pay it.
Spatial fields
Areal data carry a neighbourhood structure: regions that sit next to each other tend to look alike, and a regression that treats them as independent draws gets the standard errors wrong and lets a spatially patterned covariate absorb the regional signal. A spatial field is the place to put that signal. You add one random effect per region, tie neighbouring effects together so the field stays smooth, and read the covariate against a background that already accounts for location.
Pass the adjacency through spatial= and name the per-observation unit with
a spatial(col) term in the formula. type = "icar" puts an intrinsic
conditional-autoregressive field on the units. Here twenty regions sit on a
ring, each adjacent to its two neighbours, with a smooth wave of regional
effect around the loop.
K <- 20 W <- matrix(0, K, K) for (i in 1:K) { j <- if (i < K) i + 1 else 1; W[i, j] <- W[j, i] <- 1 } region <- factor(sample(1:K, n, replace = TRUE)) field <- as.numeric(scale(sin(2 * pi * (1:K) / K)))[region] ds <- data.frame(y = rbinom(n, 1, plogis(-0.2 + 0.7 * x + field)), x = x, region = region) fit_sp <- tulpa(y ~ x + spatial(region), data = ds, family = "binomial", spatial = list(type = "icar", adjacency = W), mode = "laplace") coef(fit_sp)
The slope recovers near its true 0.7 while the ICAR field absorbs the
regional structure. Plot the per-region field effect around the ring to see
the wave the fit pulled out of binary data: field_effects() returns the
field at each region, here its posterior mode, and traces the smooth signal
the covariate could not explain.
library(ggplot2) fe <- field_effects(fit_sp) ggplot(data.frame(region = fe$unit, effect = fe$estimate), aes(region, effect)) + geom_line(linewidth = 1, colour = "steelblue") + geom_point() + labs(x = "region", y = "field effect") + theme(panel.background = element_rect(fill = "transparent"), plot.background = element_rect(fill = "transparent"))
The field absorbs the regional signal so the slope does not have to. Without the field, a covariate that happens to vary across the ring the same way the latent wave does would soak up part of that wave and report a biased slope. With the field carrying the spatial structure, the covariate reports its own effect, and the standard error widens to its honest size because the field, not the slope, now carries the autocorrelation. The intrinsic prior fixes the field only up to an additive constant, so the intercept floats against a near-flat prior and its standard error runs wide by design; read the slopes, which are identified, and treat the intercept as a nuisance level.
Continuous fields use a coordinate spec instead of an adjacency:
spatial_gp(~ lon + lat) for a nearest-neighbour GP,
spatial_gp(~ lon + lat, approx = "hsgp") for a Hilbert-space GP,
spatial_spde(~ lon + lat, data) for a Matern
SPDE solved on a mesh through fit_spde(). The spatial vignette walks
through the areal path end to end, recovers a known field, integrates the
smoothing precision rather than fixing it, and shows how the field's
uncertainty flows into the fixed-effect standard errors.
Temporal fields
Datasets that carry a clock carry the temporal twin of spatial confounding.
A smooth background that drifts over time will be read as covariate signal if
the covariate trends the same way, and the slope comes out too steep. A
flexible function of time soaks up the background and leaves the slope to
explain only what wiggles faster than the trend. tulpa offers that function
as a first-order random walk, declared with temporal_rw1() and passed
through the temporal= argument.
Simulate counts over thirty time points with a smooth sinusoidal trend, a real covariate effect, and Poisson noise, then fit with the random-walk trend in place.
Tt <- 30 time <- sample(seq_len(Tt), n, replace = TRUE) trend <- 1.2 * sin(2 * pi * seq_len(Tt) / Tt) xt <- rnorm(n) dt <- data.frame(y = rpois(n, exp(0.4 + 0.5 * xt + trend[time])), x = xt, time = time) fit_t <- tulpa(y ~ x, data = dt, family = "poisson", temporal = temporal_rw1("time"), mode = "auto") coef(fit_t)
The slope lands near its true 0.5 with the trend absorbed into the field
rather than the coefficient. mode = "auto" routes the temporal field to the
nested-Laplace backend, and tulpa_outer_grid() returns its grid: the
per-grid latent modes in modes (one row per grid point, columns ordered
[fixed effects, field]) and the grid weights. Weight-average the field
column for each time point to integrate over the smoothing precision, then
plot the recovered walk against the truth.
g <- tulpa_outer_grid(fit_t) te <- colSums(g$weights * g$modes)[g$n_fixed + seq_len(Tt)] ggplot(data.frame(time = seq_len(Tt), est = te - mean(te), truth = trend - mean(trend)), aes(time)) + geom_line(aes(y = truth), linewidth = 1, colour = "grey50") + geom_line(aes(y = est), linewidth = 1, colour = "steelblue") + labs(x = "time", y = "trend (centred)") + theme(panel.background = element_rect(fill = "transparent"), plot.background = element_rect(fill = "transparent"))
The estimated walk follows the wave, centred to remove the level the
sum-to-zero handling leaves free. The front door carries temporal_rw1()
today; the second-order walk, the AR(1) process, and per-group panel trends
have direct fitters the temporal vignette covers, which walks the full path
from model to recovery to honest uncertainty.
Priors
Every Bayesian fit carries priors whether you write them down or not. tulpa
keeps them in one tulpa_priors object that names a prior for each kind of
parameter the engine knows about, built with tulpa_priors() and filled with
prior_*() constructors. The same object drives prior predictive simulation,
so the priors you check are the priors you fit with. Three jobs run through
the prior system: a fixed effect lives on the link scale and gets a prior over
the whole real line, a variance is positive and gets a prior on the
half-line, a proportion is trapped in the unit interval.
Look at what a fixed-effect prior implies before any data touch it.
prior_predict() draws coefficients from their priors, pushes them through
the linear predictor, and simulates responses, so you can see whether the
prior expects data the response could never produce.
fam <- tulpa_family( "poisson", function(eta, params, n_obs, ...) rpois(n_obs, exp(eta[[1]])) ) pp <- prior_predict(y ~ x, family = fam, data = dp, n_draws = 50, priors = tulpa_priors(beta = prior_normal(0, 1))) pp
pp$y holds one simulated dataset per prior draw. Overlay their densities to
see the range the prior expects, then tighten the prior if it covers ranges
the response could never take.
sims <- data.frame( value = unlist(pp$y), draw = factor(rep(seq_along(pp$y), lengths(pp$y))) ) ggplot(sims, aes(value, group = draw)) + geom_density(colour = "steelblue", alpha = 0.3) + labs(x = "simulated y", y = "density") + theme(panel.background = element_rect(fill = "transparent"), plot.background = element_rect(fill = "transparent"))
A wide prior predictive that spills past the plausible range says the priors
are too vague; tightening prior_normal(0, 2.5) toward prior_normal(0, 1)
pulls it in. The priors vignette covers all seven constructors, the
penalised-complexity prior prior_pc() for scale parameters, the defaults,
and how to move a prior and watch the posterior follow.
Custom latent blocks with tgmrf()
When the latent structure you need is not a built-in field, tgmrf() lets a
script define its own Gaussian Markov random field and plug it in as a
latent(...) term that reaches every inference tier. You supply two closures
and an init vector: Q(theta) returns the sparse precision matrix at a given
hyperparameter, prior(theta) returns its log-density, and init sets the
starting hyperparameter. No autodiff, no DSL, no codegen; the engine calls
the closures numerically and owns everything from the Newton system up.
Define a first-order random-walk precision on a chain of twenty nodes, the same structure the temporal path builds internally, but written by hand to show the contract.
library(Matrix) m <- 20 D2 <- diff(diag(m), differences = 1) R <- crossprod(D2) # RW1 structure matrix blk <- tgmrf( Q = function(theta) as(theta[1] * (R + 1e-4 * diag(m)), "dgCMatrix"), prior = function(theta) dgamma(theta[1], 2, 1, log = TRUE), init = c(tau = 1) ) blk
The constructor evaluates Q(init) once to catch errors early, infer the
field length, and capture the sparsity pattern the inner solver reuses at
every grid point. The returned block carries the precision family, the
hyperparameter names, and the registry key the nested-Laplace dispatcher
recognises, and prints a compact summary of what it found.
Q1 <- as.matrix(blk$Q(c(tau = 1))) ix <- which(Q1 != 0, arr.ind = TRUE) ggplot(data.frame(row = ix[, 1], col = ix[, 2]), aes(col, row)) + geom_tile(fill = "steelblue") + scale_y_reverse() + labs(x = "column", y = "row") + theme(panel.background = element_rect(fill = "transparent"), plot.background = element_rect(fill = "transparent"))
The tridiagonal pattern is the Markov property written as a matrix: each node
is conditionally independent of every non-neighbour given its neighbours, and
that independence is exactly a zero off the three central diagonals. The
tgmrf() vignette builds a periodic AR(1) block end to end, fits it through
the full tier ladder from nested Laplace to exact sampling, and shows how one
block written once reaches all four tiers without modification.
S3 methods
The fitted object answers the standard accessors, and they work the same way on a Laplace fit and a sampler fit. The Laplace tier reports the Gaussian approximation it computed; the sampler tier reports empirical posterior summaries from its draws. Reaching for the accessor you already know from base R or broom is the intended workflow.
| Accessor | Returns |
|---|---|
| coef() | named fixed-effect estimates |
| summary() | the full report: backend and tier, fixed effects with SE and credible bounds, variance components, hyperparameters, reliability; coef(summary(fit)) is the fixed-effect table |
| confint() | credible intervals (matrix) |
| vcov() | fixed-effect covariance |
| ranef() | group-level random effects |
| tidy() | one-row-per-term data frame (broom style) |
| glance() | one-row model summary |
| diagnostics() | convergence or approximation reliability, per parameter |
| tulpa_outer_grid() | a nested-Laplace fit's hyperparameter grid: nodes, weights, per-cell modes |
| predict(), fitted() | fixed-effect prediction / in-sample mean |
| logLik() | log marginal likelihood (Laplace) or mean log-density |
| plot() | fixed-effect densities (or trace for sampler fits) |
tidy() returns the fixed-effect table in the broom convention: one row per
term, with the estimate, standard error, and credible bounds in named
columns, ready to feed a table or a coefficient plot.
tidy(fit_re)
glance() returns a single-row model summary carrying the fixed-effect
count, the number of draws, the log-likelihood, the divergence count, the
mean acceptance rate, and the convergence flag. On a Laplace fit the
sampler-only fields read NA, which is the honest report; on a sampler fit
they carry the diagnostics that say whether the chain mixed.
glance(fit_mala)[c("n_samples", "logLik", "mean_accept", "n_divergent")]
plot() defaults to the fixed-effect densities, a quick visual of where each
coefficient sits and how wide its posterior is. On a sampler fit type =
"trace" shows the chains instead, as in the inference section above, and
type = "pairs" shows the joint posterior of the leading parameters.
coef(), confint(), vcov(), and ranef() round out the set, each
reading off the fit that has already run rather than triggering new
computation. Because the model packages inherit class = c("model_fit",
"tulpa_fit"), every accessor here works unchanged on an occupancy fit, a
ratio fit, or a GLMM fit, with the package adding only its own model-specific
methods on top.
Practical guidance
A few rules of thumb for everyday use:
- Start at Tier 2.
mode = "laplace"is deterministic and returns in well under a second on tens of thousands of rows. Use it for the first look and for model comparison by marginal likelihood. Move to a sampler (mode = "mala") when you need exact posterior moments or the posterior is visibly non-Gaussian, and check first withmode = "imh_laplace", whose acceptance rate is a direct verdict on the Laplace approximation. sigma_reis conditioned, not estimated, on the fixed-sigma_repaths. With a random intercept the fit assumes the value you pass. If you do not have one, fit at a couple of values (say 0.3 and 1.0) and check whether the fixed effects move; if they do, integrate the variance through the nested-Laplace path rather than guessing. Random slopes route there on their own.- Sample size. A random-effect term needs enough groups to inform its variance: with fewer than about 5 to 8 groups, the group-level estimates are dominated by the prior, and a fixed effect for the grouping is often the honest choice. Spatial fields want dozens to thousands of regions, many of them thin, to earn their machinery.
- Watch convergence on the sampler tiers.
glance()reportsmean_acceptand the divergence count; an acceptance rate far from 0.4 to 0.8, or any divergences, means the chain needs more warmup or a smaller step size (control$epsilon). A Tier-1 fit is only nominally exact until the chain has been shown to mix. - Reach for a model package when the likelihood is not built in. The
five families here cover the common cases. Zero-inflation, hurdles,
detection and occupancy, ratio responses, and posterior-predictive
checking live in the packages built on this engine, not in
tulpa()itself, and they inherit every tier and accessor shown here.
Where to go next
- Inference details and the tier guarantees: the inference-modes and algorithm vignettes.
- Spatial and temporal fields, spatially varying coefficients: the spatial, temporal, and SVC vignettes.
- Priors and prior predictive checks: the priors vignette and
prior_predict(). -
User-defined latent structure: the
tgmrf()vignette. -
Model comparison by marginal likelihood, WAIC, and LOO: the model-comparison vignette.
- Posterior-predictive checks, WAIC, and LOO on a full observation model: the model packages (tulpaObs, tulpaRatio, tulpaGlmm) that build on this engine. ```
Try the tulpa package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.