Tutorial 8: Analysis Pipelines

The formula surface compiles a regression specification and a dataframe into a QVR program. This chapter inspects that program, fits hierarchical and Bernoulli regressions with SVI and NUTS, compares models with PSIS-LOO, and runs a posterior-predictive check on synthetic data.

Concepts

• Formula: A brms / lme4-style string like "y ~ x + (1 | g)" that names a response, predictors, and random-effect structure.
• fit: The one-liner that parses the formula against a dataframe, compiles to a QVR MonadicProgram, and dispatches to NUTS / HMC / SVI.
• BayesianFit: The returned record. Carries the parsed Formula, the compiled program, the posterior (either an MCMCResult or a Guide), and the inference-time observations dict.
• Posterior: For NUTS / HMC, an MCMCResult with (chain, draw)-shaped sample tensors per site. For SVI, a fitted Guide from which you draw via rsample.
• DataTree: The ArviZ 1.x container used by the diagnostics layer. to_datatree accepts an MCMCResult; to_datatree_from_svi accepts SVI samples, and to_datatree_any dispatches between supported inputs.

Setup

The analysis-pipeline surface is gated behind the formulas, data, and diagnostics extras:

pip install "quivers[formulas,data,diagnostics]"

import numpy as np
import pandas as pd
import torch
from quivers.formulas import fit, formula_to_qvr
from quivers.diagnostics import to_datatree, compare, posterior_predictive_check

A first fit: Gaussian regression in one line

Generate 200 rows of synthetic data with a known linear relationship:

torch.manual_seed(0)
np.random.seed(0)

N = 200
x = np.random.randn(N)
y = 1.0 + 2.0 * x + 0.5 * np.random.randn(N)
df = pd.DataFrame({"y": y, "x": x})

Fit y ~ x under a Gaussian family. SVI is the fastest inference method (it fits a variational Guide by optimization rather than running an MCMC chain) and is enough for a quick sanity check:

result = fit(
    "y ~ x",
    data=df,
    family="gaussian",
    method="svi",
    num_samples=2000,    # under SVI, num_samples is the number of optimisation steps
    seed=0,
)

result is a BayesianFit: a frozen dx.Model with these load-bearing fields:

Field	What it holds
.formula	The parsed Formula: columns, family, response.
.program	The compiled MonadicProgram.
.posterior	An MCMCResult (NUTS / HMC) or a Guide (SVI).
.observations	The inference-time observations dict (response values + host-data covariates + plate-index codes).

Inspect the QVR program the formula compiled to:

print(result.qvr_source)

object Resp : FinSet 200
program model : Resp -> Resp
    sample intercept <- Normal(0.0, 5.0)
    sample beta_x <- Normal(0.0, 5.0)
    let beta_x_per_row = (beta_x * x)
    sample sigma <- HalfCauchy(2.0)
    let eta = (intercept + beta_x_per_row)
    let mu = eta
    observe y : Resp <- Normal(mu, sigma)
    return y

export model

Per-row predictor values do not appear in the program: the (beta_x * x) term references a free variable x that the inference path resolves from the observations dict at fit time (the host-data channel; see quivers.inference.trace). The intercept latent stays unindexed because the constant-1 design column doesn't need a per-row binding.

You can also emit the source without fitting:

src = formula_to_qvr("y ~ x", data=df)
result.dump_qvr("regression.qvr")          # write the same source to disk

Recover the coefficients

Under SVI the posterior is an AutoNormalGuide. Draw from it and check the means against the true values (intercept=1.0, beta_x=2.0):

draws = [result.posterior.rsample(torch.zeros(1, 1)) for _ in range(200)]
ic = torch.stack([d["intercept"] for d in draws])
bx = torch.stack([d["beta_x"] for d in draws])
print(f"intercept (true 1.0): {ic.mean().item():.3f} ± {ic.std().item():.3f}")
print(f"beta_x    (true 2.0): {bx.mean().item():.3f} ± {bx.std().item():.3f}")

intercept (true 1.0): 0.927 ± 0.036
beta_x    (true 2.0): 2.023 ± 0.030

These printed values record one seeded, short optimization run; they are not a convergence guarantee. Increase the optimization budget and inspect diagnostics, or use a diagnosed MCMC fit, before interpreting coefficient recovery.

Random effects: hierarchical model under SVI

Eight groups with their own intercept offsets, twelve rows per group:

torch.manual_seed(0)
np.random.seed(0)

N_g, N_per = 8, 12
g = np.repeat(np.arange(N_g), N_per)
alpha = np.random.randn(N_g) * 0.6   # group-level offsets
x = np.random.randn(N_g * N_per)
y = 1.5 + 0.7 * x + alpha[g] + 0.4 * np.random.randn(N_g * N_per)

df = pd.DataFrame({"y": y, "x": x, "g": g})

The formula adds a (1 | g) random intercept:

result = fit(
    "y ~ x + (1 | g)",
    data=df,
    family="gaussian",
    method="svi",
    num_samples=4000,
    seed=0,
)
print(result.qvr_source)

object Resp : FinSet 96
object g : FinSet 8
program model : Resp -> Resp
    sample intercept <- Normal(0.0, 5.0)
    sample beta_x <- Normal(0.0, 5.0)
    let beta_x_per_row = (beta_x * x)
    sample sigma_g_Intercept <- HalfNormal(1.0)
    sample z_g_Intercept : g <- Normal(0.0, 1.0)
    let alpha_g_per_row = (sigma_g_Intercept * z_g_Intercept[g_idx])
    sample sigma <- HalfCauchy(2.0)
    let eta = ((intercept + beta_x_per_row) + alpha_g_per_row)
    let mu = eta
    observe y : Resp <- Normal(mu, sigma)
    return y

export model

The random-intercept structure expands to a HalfNormal scale latent sigma_g_Intercept, a standard-Normal plate draw z_g_Intercept of size |g|, and a per-row contribution sigma_g_Intercept * z_g_Intercept[g_idx]. This is the non-centered parameterization: alpha_g[i] = sigma * z[i] is computed deterministically rather than drawn as alpha_g[i] ~ Normal(0, sigma). Non-centering often improves posterior geometry in weakly informed hierarchical models.

Recovery:

def mean_param(res, name, n=200):
    draws = torch.stack(
        [res.posterior.rsample(torch.zeros(1, 1))[name] for _ in range(n)]
    )
    return draws.mean().item(), draws.std().item()

for name, true in [
    ("intercept", 1.5),
    ("beta_x", 0.7),
    ("sigma_g_Intercept", 0.6),
    ("sigma", 0.4),
]:
    m, s = mean_param(result, name)
    print(f"{name:20s} (true {true:.2f}): {m:.3f} ± {s:.3f}")

intercept            (true 1.50): 1.748 ± 0.035
beta_x               (true 0.70): 0.620 ± 0.038
sigma_g_Intercept    (true 0.60): 0.273 ± 0.014
sigma                (true 0.40): 0.375 ± 0.028

This short SVI fit recovers the fixed effects more closely than the group-level scale. Mean-field VI can underestimate posterior uncertainty and distort variance-component estimates. For a richer approximation, pass guide=AutoMultivariateNormalGuide; for sampling-based inference, use method="nuts" and inspect its diagnostics.

Knobs: method, guide, reparameterize

The following choices control the inference path:

• method="nuts" | "hmc" | "svi". The default "nuts" samples the joint posterior with the No-U-Turn extension to HMC and produces an MCMCResult. "svi" fits a Guide by optimization and may be more practical for large models. Either method still requires diagnostics.
• guide=Cls (SVI only). Defaults to AutoNormalGuide (mean-field diagonal Normal). Other choices: AutoMultivariateNormalGuide for a full-rank Cholesky; AutoLowRankMultivariateNormalGuide for O(D·rank) memory on high-dimensional models; AutoLaplaceApproximation for a Gaussian fit at the MAP; AutoIAFGuide for an inverse-autoregressive normalizing flow. Pass the class, not an instance: fit(..., guide=AutoMultivariateNormalGuide).
• reparameterize="noncentered" | "centered". formula_to_qvr uses this option when emitting random-effect terms. The current fit implementation records the option on BayesianFit but does not pass it into the runtime compiler, which thus uses non-centering. Use formula_to_qvr(..., reparameterize="centered") when you need centered source and inspect result.qvr_source to confirm the fitted program.

Model comparison with PSIS-LOO

Two models on the same Bernoulli outcome: with and without the predictor. The smaller N and short MCMC keeps the example fast.

torch.manual_seed(0)
np.random.seed(0)

N = 80
x = np.random.randn(N)
logit = -0.5 + 1.2 * x
y = np.random.binomial(1, 1 / (1 + np.exp(-logit))).astype(float)
df = pd.DataFrame({"y": y, "x": x})

res_with_x = fit(
    "y ~ x", data=df, family="bernoulli", method="nuts",
    num_warmup=100, num_samples=200, num_chains=2, seed=0,
)
res_intercept = fit(
    "y ~ 1", data=df, family="bernoulli", method="nuts",
    num_warmup=100, num_samples=200, num_chains=2, seed=0,
)

Compute per-row log-likelihoods under each model's posterior (the ArviZ comparison needs them):

y_obs = torch.as_tensor(y, dtype=torch.float32)
x_t = torch.as_tensor(x, dtype=torch.float32)

def per_row_log_lik(res, y_obs, x_t):
    samples = res.posterior.samples
    ic = samples["intercept"]                         # (chain, draw)
    bx = samples.get("beta_x", torch.zeros_like(ic))  # zero if intercept-only
    logits = ic[..., None] + bx[..., None] * x_t
    p = torch.sigmoid(logits)
    return y_obs * torch.log(p + 1e-12) + (1 - y_obs) * torch.log(1 - p + 1e-12)

ll_with_x = per_row_log_lik(res_with_x, y_obs, x_t)
ll_intercept = per_row_log_lik(res_intercept, y_obs, x_t)

Wrap each posterior into an ArviZ DataTree with the canonical groups:

idata_with_x = to_datatree(
    res_with_x.posterior,
    observed_data={"y": y_obs},
    log_likelihood={"y": ll_with_x},
)
idata_intercept = to_datatree(
    res_intercept.posterior,
    observed_data={"y": y_obs},
    log_likelihood={"y": ll_intercept},
)

compare delegates to ArviZ's PSIS-LOO ranking and stacking weights:

print(compare({"with_x": idata_with_x, "intercept_only": idata_intercept}))

                rank  elpd    p  elpd_diff  weight   se  dse  warning
with_x             0 -43.0  0.0        0.0    0.98  4.4  0.0    False
intercept_only     1 -50.0  0.0       -6.0    0.02  3.3  3.4    False

In this illustrative output, with_x receives most of the stacking weight. Interpret the comparison from the values produced by your run, including the standard errors and Pareto-\(k\) warnings.

Posterior-predictive checks

posterior_predictive_check computes the posterior-predictive p-value for a chosen test statistic. Generate replicated responses under each posterior draw of (intercept, beta_x) and check whether the observed mean is in the bulk of the replicated distribution:

samples = res_with_x.posterior.samples
ic, bx = samples["intercept"], samples["beta_x"]
p_pp = torch.sigmoid(ic[..., None] + bx[..., None] * x_t)
pp = torch.bernoulli(p_pp).float()         # (chain, draw, N) replicated y's

idata = to_datatree(
    res_with_x.posterior,
    observed_data={"y": y_obs},
    posterior_predictive={"y": pp},
)
ppc = posterior_predictive_check(idata, observed_name="y", statistic="mean")
print(f"observed mean: {float(ppc['observed']):.3f}")
print(f"predictive mean: {ppc['predictive_mean']:.3f}")
print(f"ppp (mean):    {ppc['ppp']:.3f}")

observed mean: 0.312
predictive mean: 0.312
ppp (mean):    0.557

A posterior-predictive p-value near 0.5 places the observed mean near the middle of the replicated distribution for this statistic. Values near 0 or 1 can indicate that the model does not capture the chosen statistic, though the interpretation depends on how the check was constructed.

Built-in statistics: mean, median, sd, var, min, max, q05, q95. The by= argument splits the calculation along a named dim (e.g. by="Verb" for a per-category check); for arbitrary statistics, pass a callable.

What about the original formula?

BayesianFit keeps both the parsed Formula and the compiled program. You can round-trip:

print(result.formula.formula)         # original formula string
print(result.formula.response_name)   # "y"
for col in result.formula.fixed_columns:
    print(f"  {col.qvr_name}: term={col.term!r}, is_intercept={col.is_intercept}")

The lens behind the scenes (FormulaToQVRModule) is a typed dx.Lens whose forward emits (Module, FormulaData). FormulaData is the strict subset of the formula that isn't encoded in the QVR program: the per-row data arrays, the original (pre-_qvr_name) identifier names, presentation labels, the formula string itself. Calling lens.backward(module, complement) decodes the structural skeleton out of the Module and fuses it with the complement to reproduce the original Formula byte-for-byte. See the formulas API reference for the full lens signature.

Summary

You have:

• Fit a Gaussian regression with one line of code and recovered the true coefficients via SVI.
• Inspected the emitted .qvr source and seen how predictors flow in through the host-data channel rather than as latent draws.
• Fit a hierarchical model with (1 | g) random intercepts and seen where SVI's mean-field approximation under-shrinks group-level variance.
• Run NUTS on a Bernoulli regression and compared two models via PSIS-LOO with stacking weights.
• Wrapped MCMC results into an ArviZ DataTree and run a posterior-predictive check.
• Inspected the lens machinery that maps a Formula to a QVR Module.

Next

• The Analysis Pipelines guide is the reference for the whole quivers.formulas / quivers.data / quivers.diagnostics surface, including the full prior-override syntax and the DatasetSchema bridge from dataframes to QVR programs without writing a formula.
• The Variational Inference tutorial shows the lower-level path: building a MonadicProgram by hand, tracing, and running SVI without the formula surface.
• The DSL guide is the right place to go if you want to write or hand-edit .qvr programs directly.