inlabru as a Bayesian GLM with parsnip • inlaops

Scope

bru_reg() is scoped to inlabru used as a fast Bayesian GLM: models with fixed effects and IID random intercepts only. If your model includes any of the following, stop here and use inlabru’s native interface instead:

Spatial random effects: BYM2, SPDE, LGCP
Structured temporal effects: AR1, RW1, RW2
Custom nonlinear predictors with mapper functions

For those models, cross-validation RMSE is not a meaningful quantity and the parsnip abstraction does not hold. See vignette("scoring") for the statistically appropriate evaluation approach.

For simple Bayesian GLMs — Poisson regression for counts, Gaussian regression for continuous outcomes, binomial for proportions — the parsnip layer gives you a consistent interface for comparison with frequentist or ML baselines and clean integration with the rest of the tidymodels ecosystem.

Setup

library(inlabru)
library(inlaops)
library(parsnip)
library(broom)
library(rsample)
library(yardstick)
library(dplyr)

Example: Poisson regression for count data

set.seed(123)
n   <- 300
dat <- data.frame(
  area   = rnorm(n),
  effort = runif(n, 0.5, 2.0)
)
log_mu     <- 0.5 + 0.8 * dat$area + log(dat$effort)
dat$counts <- rpois(n, exp(log_mu))

bru_reg() does not auto-encode factor or character predictors — fitting counts ~ Intercept(1) + habitat on a raw factor or character column errors with a pointer to three workarounds: pre-encode via recipes::step_dummy() + recipes::bake() before calling fit() (recipes/workflows cannot be wrapped around bru_reg because they rebuild formulas without inlabru’s component syntax), convert to integer indicator columns yourself, or use inlabru’s native habitat(habitat, model = "factor_full") component syntax.

Fit with bru_reg()

spec <- bru_reg(mode = "regression", family = "poisson")

fit <- spec |>
  set_engine("inlabru") |>
  fit(counts ~ Intercept(1) + area + effort, data = dat)

fit
#> parsnip model object
#> 
#> inlabru version: 2.15.0 
#> INLA version: 26.08.07 
#> Latent components:
#> Intercept: main = linear(1)
#> area: main = linear(area)
#> effort: main = linear(effort)
#> Observation models:
#>   Model tag: <No tag>
#>     Family: 'poisson'
#>     Data class: 'data.frame'
#>     Response class: 'integer'
#>     Predictor: counts ~ Intercept + area + effort
#>     Additive/Linear/Rowwise: TRUE/TRUE/TRUE
#>     Used components: effect[Intercept, area, effort], latent[] 
#> Time used:
#>     Pre = 0.188, Running = 0.336, Post = 0.0734, Total = 0.598 
#> Fixed effects:
#>             mean    sd 0.025quant 0.5quant 0.975quant   mode kld
#> Intercept -0.271 0.127     -0.519   -0.271     -0.023 -0.271   0
#> area       0.762 0.035      0.694    0.762      0.829  0.762   0
#> effort     0.751 0.086      0.583    0.751      0.920  0.751   0
#> 
#> Deviance Information Criterion (DIC) ...............: 994.69
#> Deviance Information Criterion (DIC, saturated) ....: 292.49
#> Effective number of parameters .....................: 3.00
#> 
#> Watanabe-Akaike information criterion (WAIC) ...: 994.85
#> Effective number of parameters .................: 3.02
#> 
#> Marginal log-Likelihood:  -513.88 
#>  is computed 
#> Posterior summaries for the linear predictor and the fitted values are computed
#> (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

The engine automatically sets control.compute = list(config = TRUE, dic = TRUE, waic = TRUE). The fit is serving-ready for the vetiver bridge without any extra options. For large data sets where WAIC is too expensive, opt out via set_engine("inlabru", dic = FALSE, waic = FALSE); config = TRUE is always set because the vetiver /predict endpoint needs it.

Inspect with broom

The tidy(), glance(), and augment() methods for bru fits come from inlabru (>= 2.15.0).

tidy(fit$fit)
#> # A tibble: 3 × 5
#>   term      estimate std.error conf.low conf.high
#>   <chr>        <dbl>     <dbl>    <dbl>     <dbl>
#> 1 Intercept   -0.271    0.127    -0.519   -0.0227
#> 2 area         0.762    0.0346    0.694    0.829 
#> 3 effort       0.751    0.0860    0.583    0.920

glance(fit$fit)
#> # A tibble: 1 × 5
#>     dic  waic marginal_loglik  nobs elapsed
#>   <dbl> <dbl>           <dbl> <int>   <dbl>
#> 1  995.  995.           -514.   300   0.695

Predict

new_dat <- data.frame(area = c(-1, 0, 1), effort = c(1, 1, 1))
predict(fit, new_data = new_dat)
#> # A tibble: 3 × 1
#>   .pred
#>   <dbl>
#> 1 0.757
#> 2 1.62 
#> 3 3.46

The parsnip predict() returns the posterior mean on the response scale. For predictive uncertainty, use augment() directly:

augment(fit$fit, data = new_dat, pred_formula = ~ exp(Intercept + area + effort))
#> # A tibble: 3 × 6
#>    area effort .fitted .fitted_low .fitted_high .fitted_sd
#>   <dbl>  <dbl>   <dbl>       <dbl>        <dbl>      <dbl>
#> 1    -1      1   0.749       0.643        0.868     0.0583
#> 2     0      1   1.61        1.45         1.77      0.0853
#> 3     1      1   3.44        3.14         3.73      0.157

Comparing against a frequentist baseline

A side-by-side comparison with a frequentist GLM is a one-liner per engine. Note that bru_reg requires inlabru component syntax (Intercept(1) is explicit), while glm() uses the standard implicit-intercept form — the formulas differ on purpose.

# Frequentist baseline
fit_glm <- glm(counts ~ area + effort, data = dat, family = poisson())

# Bayesian with inlabru, via parsnip
bru_spec <- bru_reg(family = "poisson") |> set_engine("inlabru")
fit_bru  <- bru_spec |>
  fit(counts ~ Intercept(1) + area + effort, data = dat)

# Coefficient comparison (response-scale interval for the Bayesian fit)
broom::tidy(fit_glm)
#> # A tibble: 3 × 5
#>   term        estimate std.error statistic   p.value
#>   <chr>          <dbl>     <dbl>     <dbl>     <dbl>
#> 1 (Intercept)   -0.268    0.127      -2.12 3.42e-  2
#> 2 area           0.761    0.0346     22.0  1.93e-107
#> 3 effort         0.751    0.0860      8.73 2.55e- 18
broom::tidy(fit_bru$fit)
#> # A tibble: 3 × 5
#>   term      estimate std.error conf.low conf.high
#>   <chr>        <dbl>     <dbl>    <dbl>     <dbl>
#> 1 Intercept   -0.271    0.127    -0.519   -0.0227
#> 2 area         0.762    0.0346    0.694    0.829 
#> 3 effort       0.751    0.0860    0.583    0.920

recipes and workflows are not compatible with bru_reg because they rebuild formulas in the standard implicit-intercept shape, which bru() rejects (it needs the explicit Intercept(1) component). For preprocessing, manipulate data before the call to fit() rather than wrapping in a recipe.

Cross-validation with bru_cv_predict()

For fixed-effects models, bru_cv_predict() wraps vfold_cv() and returns a tidy holdout prediction tibble:

cv_preds <- bru_cv_predict(
  formula      = counts ~ Intercept(1) + area + effort,
  family       = "poisson",
  data         = dat,
  v            = 5L,
  pred_formula = ~ exp(Intercept + area + effort),
  n_samples    = 500L
)

cv_preds
#> # A tibble: 300 × 5
#>    .truth .pred .pred_low .pred_high .pred_sd
#>     <int> <dbl>     <dbl>      <dbl>    <dbl>
#>  1      2 1.91      1.71       2.14    0.107 
#>  2     11 9.95      8.82      11.1     0.603 
#>  3      3 5.62      5.17       6.06    0.236 
#>  4      6 2.07      1.82       2.36    0.133 
#>  5      2 2.19      1.94       2.48    0.133 
#>  6      6 2.86      2.52       3.26    0.193 
#>  7      2 2.02      1.80       2.27    0.117 
#>  8      1 0.957     0.796      1.14    0.0872
#>  9      0 0.314     0.247      0.394   0.0380
#> 10      5 2.65      2.24       3.11    0.212 
#> # ℹ 290 more rows

# Compute RMSE across folds
yardstick::rmse(cv_preds, truth = .truth, estimate = .pred)
#> # A tibble: 1 × 3
#>   .metric .estimator .estimate
#>   <chr>   <chr>          <dbl>
#> 1 rmse    standard        1.63

For proper scoring rules that use the full posterior distribution rather than just the mean, see vignette("scoring").

Connecting to vetiver

A model fit via bru_reg() is automatically serving-ready:

library(vetiver)
library(pins)

v <- vetiver_model(fit$fit, "poisson_glm")
v$metadata$user <- list(
  pred_formula = "~ exp(Intercept + area + effort)",
  n_samples    = 1000L
)

board <- board_temp()
vetiver_pin_write(board, v)
#> Creating new version '20261016T203236Z-0bd13'
#> Writing to pin 'poisson_glm'
#> 
#> Create a Model Card for your published model
#> • Model Cards provide a framework for transparent, responsible reporting
#> • Use the vetiver `.Rmd` template as a place to start

See vignette("vetiver") for the full serving and monitoring workflow.