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)
#> (Intercept) x
#> 0.4956826 1.1717667The 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, 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)
#> estimate std.error 2.5% 97.5%
#> (Intercept) 0.4956826 0.04000155 0.417281 0.5740842
#> x 1.1717667 0.04087976 1.091644 1.2518896The 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)
#> 2.5% 97.5%
#> (Intercept) 0.417281 0.5740842
#> x 1.091644 1.2518896The 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)
#> (Intercept) x
#> 0.3719796 1.1491057coef() 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)
#> term estimate sd conf.low conf.high source
#> 1 g[1] -1.1846742 NA NA NA mode
#> 2 g[2] -0.6793481 NA NA NA mode
#> 3 g[3] -0.8067578 NA NA NA mode
#> 4 g[4] -0.4506177 NA NA NA modeThis 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")
#> model n_params logLik
#> 1 intercept 1 -1059.0421
#> 2 slope 2 -672.4025
#> 3 slope_re 2 -496.2809The 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)
#> fit se.fit lower upper
#> 1 -1.926232 0.1962305 -2.310836 -1.541627
#> 2 -1.832427 0.1948400 -2.214307 -1.450548
#> 3 -1.738623 0.1934976 -2.117871 -1.359374The 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"))
#> (Intercept) x
#> 0.2170233 0.5695186
# 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))
#> (Intercept) x
#> 0.1752357 0.8561175The 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.
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)
)
#> term laplace mala
#> (Intercept) (Intercept) -0.472 -0.613
#> x x 1.089 1.091The agreement here is itself a signal that Laplace was safe for this
model. Every fit records which backend ran and at which tier, through
the backend and inference_tier fields, and
selection_reason records why mode = "auto"
chose what it did, so you are never left guessing which engine produced
a number.
c(backend = fit_mala$backend, tier = fit_mala$inference_tier)
#> backend tier
#> "mala" "1"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)
#> (Intercept) x
#> -0.008787351 0.751041095The 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: the mode at each region’s latent column traces the smooth signal the covariate could not explain.
library(ggplot2)
fe <- tail(fit_sp$mode, K)
ggplot(data.frame(region = seq_len(K), effect = fe), 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)
#> (Intercept) x
#> 0.3301573 0.5502579The 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, which carries 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.
w <- fit_t$weights / sum(fit_t$weights)
nf <- fit_t$n_fixed
te <- vapply(seq_len(Tt),
function(u) sum(w * fit_t$modes[, nf + u]), numeric(1))
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
#> tulpa prior predictive draws
#> ============================
#> Family: poisson
#> Processes: y
#> Draws: 50
#> Obs: 400pp$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
#> <tgmrf>
#> n_latent : 20
#> theta : 1 dim (tau)
#> init : 1
#> Q nnz : 58 (14.5% of dense)
#> backend : rThe 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() |
fixed effects: estimate, SE, credible bounds |
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 |
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)
#> term estimate std.error conf.low conf.high
#> 1 (Intercept) 0.3719796 0.17794627 0.0232113 0.7207478
#> 2 x 1.1491057 0.04110041 1.0685504 1.2296610glance() 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")]
#> n_samples logLik mean_accept n_divergent
#> 1 375 -233.4807 0.5733333 NAplot() 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. ```