--- title: "Getting started with tulpa" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Getting started with tulpa} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 4 ) set.seed(20260529) ``` ```{r cran-gate, include = FALSE} # 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) ``` ```{r load, message = FALSE} 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. ```{r sim-gaussian} 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. ```{r fit-gaussian} fit <- tulpa(y ~ x, data = df, family = "gaussian", mode = "laplace", phi = 0.8) 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, 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. ```{r summary-gaussian} 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. ```{r confint-gaussian} 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 `phi = 0.8` is the matching value. 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. ```{r sim-re} 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) 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. ```{r ranef} 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. ```{r compare} m0 <- tulpa(y ~ 1, data = df, family = "gaussian", mode = "laplace", phi = 0.8) m1 <- tulpa(y ~ x, data = df, family = "gaussian", mode = "laplace", phi = 0.8) m2 <- tulpa(y ~ x + (1 | g), data = df, family = "gaussian", mode = "laplace", sigma_re = 0.6, phi = 0.8) 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. ```{r predict} 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. ```{r predict-plot, eval = has_ggplot && EVAL_FITS, fig.alt = "Predicted response across x with a 95 percent credible band"} 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. ```{r families} # 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. ```{r modes, eval = FALSE} 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. ```{r tiers} 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. ```{r tiers-compare} 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, 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. ```{r tier-fields} c(backend = fit_mala$backend, tier = fit_mala$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. ```{r trace-plot, fig.width = 6, fig.height = 4, fig.alt = "Trace of the MALA chain for the fixed effects"} 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. ```{r spatial-build} 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: the mode at each region's latent column traces the smooth signal the covariate could not explain. ```{r spatial-plot, eval = has_ggplot && EVAL_FITS, fig.alt = "Estimated spatial field effect per region around the ring"} 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. ```{r temporal-fit} 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, 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. ```{r temporal-plot, fig.alt = "Estimated temporal random-walk trend against the simulated 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. ```{r prior-fit} 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. ```{r prior-plot, fig.alt = "Prior predictive densities of the simulated response under the chosen priors"} 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. ```{r tgmrf-build, message = FALSE} 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. ```{r tgmrf-plot, fig.alt = "Sparsity pattern of the user-defined RW1 precision matrix"} 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. ```{r tidy} 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. ```{r glance} 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. ```