brmspy

Python-first access to R's brms with proper parameter names, ArviZ support, and cmdstanr performance. The easiest way to run brms models from Python.

Github repo and issues

Installation

R Configuration

R>=4 is required before installing brmspy.

On Linux and macOS brmspy will usually auto-detect R_HOME, and the session layer attempts to prepend $R_HOME/lib to LD_LIBRARY_PATH when needed for rpy2 ABI mode.

If you run into errors like “cannot find libR” (or similar dynamic loader issues), set these explicitly:

# Set R_HOME and add lib directory to LD_LIBRARY_PATH (Unix)
export R_HOME=$(R RHOME)
export LD_LIBRARY_PATH="${R_HOME}/lib:${LD_LIBRARY_PATH}"

# Recommended for stability
export RPY2_CFFI_MODE=ABI

Python

pip install brmspy

First-time setup (installs brms, cmdstanr, and CmdStan in R):

from brmspy import brms

with brms.manage() as ctx: # requires R to be installed already
    ctx.install_brms()

Prebuilt Runtimes (Optional)

For faster installation (~20-60 seconds vs 20-30 minutes), use prebuilt runtime bundles:

from brmspy import brms

with brms.manage() as ctx:
    ctx.install_brms(use_prebuilt=True)

Windows RTools

In case you don't have RTools installed, you can use the flag install_rtools = True. This is disabled by default, because the flag runs the full rtools installer and modifies system path. Use with caution!

from brmspy import brms

with brms.manage() as ctx:
    ctx.install_brms(
        use_prebuilt=True,
        install_rtools=True,  # works for both prebuilt and compiled installs
    )

System Requirements

Linux (x86_64):

• glibc >= 2.27 (Ubuntu 18.04+, Debian 10+, RHEL 8+)
• g++ >= 9.0
• R >= 4.0

macOS (Intel & Apple Silicon):

• Xcode Command Line Tools: xcode-select --install
• clang >= 11.0
• R >= 4.0

Windows (x86_64):

• Rtools
• R >= 4.0

Download Rtools from: https://cran.r-project.org/bin/windows/Rtools/

Key Features

• Proper parameter names: Returns b_Intercept, b_zAge, sd_patient__Intercept instead of generic names like b_dim_0
• ArviZ integration: Returns arviz.InferenceData by default for Python workflow
• brms formula syntax: Full support for brms formula interface including random effects
• Dual access: Results include .idata (ArviZ) plus a lightweight .r handle that can be passed back to brmspy to reference the underlying R object (the R object itself stays in the worker process)
• No reimplementation: Delegates all modeling logic to real brms. No Python-side reimplementation, no divergence from native behavior
• Prebuilt Binaries: Fast installation with precompiled runtimes (50x faster, ~25 seconds on Google Colab)
• Stays true to brms: Function names, parameters, and returned objects are designed to be as close as possible to brms
• Composable formula DSL: Build multivariate, non-linear, and distributional formulas by simply adding components together, identical to brms

Examples

1. Quick Start

Basic Bayesian regression with ArviZ diagnostics:

from brmspy import brms
from arviz_stats import summary

# Fit Poisson model with random effects
epilepsy = brms.get_brms_data("epilepsy")
model = brms.brm("count ~ zAge + (1|patient)", data=epilepsy, family="poisson")

# Proper parameter names automatically!
print(summary(model.idata))
#                  mean     sd  hdi_3%  hdi_97%  ...  r_hat
# b_Intercept     1.234  0.123   1.012    1.456  ...   1.00
# b_zAge          0.567  0.089   0.398    0.732  ...   1.00
# sd_patient__... 0.345  0.067   0.223    0.467  ...   1.00

2. Multivariate Models (Python vs R)

Model multiple responses simultaneously with seamless ArviZ integration:

<table> <tr><th>Python (brmspy)</th><th>R (brms)</th></tr> <tr> <td>

from brmspy import brms
from brmspy.brms import bf, set_rescor
from arviz_stats import loo
from arviz_plots import plot_ppc_dist

# Fit multivariate model
mv = brms.brm(
    bf("mvbind(tarsus, back) ~ sex + (1|p|fosternest)")
    + set_rescor(True),
    data=btdata
)

# ArviZ just works!
loo(mv.idata, var_names=["tarsus"])
loo(mv.idata, var_names=["back"])
plot_ppc_dist(mv.idata, var_names=["tarsus"])

</td> <td>

library(brms)
library(loo)

# Fit multivariate model
fit <- brm(
  bf(mvbind(tarsus, back) ~ sex + (1|p|fosternest))
  + set_rescor(TRUE),
  data = BTdata
)

# Separate LOO for each response
loo_tarsus <- loo(fit, resp = "tarsus")
loo_back <- loo(fit, resp = "back")

</td> </tr> </table>

3. Distributional Regression

Model heteroscedasticity (variance depends on predictors):

from brmspy import brms
from brmspy.brms import bf

# Model both mean AND variance
model = brms.brm(
    bf("reaction ~ days", sigma = "~ days"),  # sigma varies with days!
    data=sleep_data,
    family="gaussian"
)

# Extract distributional parameters
print(model.idata.posterior.data_vars)
# b_Intercept, b_days, b_sigma_Intercept, b_sigma_days, ...

4. Complete Diagnostic Workflow with ArviZ

Full model checking in ~10 lines:

from brmspy import brms
from arviz_stats import rhat, ess, compare
from arviz_plots import plot_ppc_dist

model = brms.brm("count ~ zAge * Trt + (1|patient)", data=epilepsy, family="poisson")

# Check convergence
assert rhat(model.idata).to_array().values.max() < 1.01, "Convergence issues!"
assert ess(model.idata).to_array().values.min() > 400, "Low effective sample size!"

# Posterior predictive check
plot_ppc_dist(model.idata, num_samples=100)

# Model comparison
model2 = brms.brm("count ~ zAge + Trt + (1|patient)", data=epilepsy, family="poisson")
comparison = compare({"interaction": model.idata, "additive": model2.idata})
print(comparison)
#              rank  loo    p_loo  d_loo  weight
# interaction     0 -456.2   12.3    0.0    0.89
# additive        1 -461.5   10.8    5.3    0.11

5. Advanced Formulas: Splines & Non-linear Effects

Smooth non-linear relationships with splines:

from brmspy import brms

# Generalized additive model (GAM) with spline
model = brms.brm(
    "y ~ s(x, bs='cr', k=10) + (1 + x | group)",
    data=data,
    family="gaussian"
)

# Polynomial regression
poly_model = brms.brm(
    "y ~ poly(x, 3) + (1|group)",
    data=data
)

# Extract and visualize smooth effects
conditional_effects = brms.call("conditional_effects", model, "x")

Additional Features

Custom Priors:

from brmspy.brms import prior

model = brms.brm(
    "count ~ zAge + (1|patient)",
    data=epilepsy,
    priors=[
        prior("normal(0, 0.5)", class_="b"),
        prior("exponential(1)", class_="sd", group="patient")
    ],
    family="poisson"
)

Predictions:

import pandas as pd
from arviz_plots import plot_dist

new_data = pd.DataFrame({"zAge": [-1, 0, 1], "patient": [999, 999, 999]})

# Expected value (without observation noise)
epred = brms.posterior_epred(model, newdata=new_data)

# Posterior predictive (with noise)
ypred = brms.posterior_predict(model, newdata=new_data)

# Access as InferenceData for ArviZ
plot_dist(epred.idata)

6. Maximalist Example: Kitchen Sink

Everything at once - multivariate responses, different families, distributional parameters, splines, and complete diagnostics:

from brmspy.brms import bf, lf, set_rescor, skew_normal, gaussian
from brmspy import brms
import arviz as az

# Load data
btdata = brms.get_data("BTdata", package="MCMCglmm")

bf_tarsus = (
    bf("tarsus ~ sex + (1|p|fosternest) + (1|q|dam)") +
    lf("sigma ~ 0 + sex") +
    skew_normal()
)

bf_back = (
    bf("back ~ s(hatchdate) + (1|p|fosternest) + (1|q|dam)") +
    gaussian()
)

model = brms.brm(
    bf_tarsus + bf_back + set_rescor(False),
    data=btdata,
    chains=2,
    control={"adapt_delta": 0.95}
)

# ArviZ diagnostics work seamlessly
for response in ["tarsus", "back"]:
    print(f"\n=== {response.upper()} ===")
    
    # Model comparison
    loo = az.loo(model.idata, var_name=response)
    print(f"LOO: {loo.loo:.1f} ± {loo.loo_se:.1f}")
    
    # Posterior predictive check
    az.plot_ppc(model.idata, var_names=[response])
    
    # Parameter summaries
    print(az.summary(
        model.idata,
        var_names=[f"b_{response}"],
        filter_vars="like"
    ))

# Visualize non-linear effect
conditional = brms.call("conditional_effects", model, "hatchdate", resp="back")
# Returns proper pandas DataFrame ready for plotting!

Output shows:

• Proper parameter naming: b_tarsus_Intercept, b_tarsus_sex, b_sigma_sex, sd_fosternest__tarsus_Intercept, etc.
• Separate posterior predictive for each response
• Per-response LOO for model comparison
• All parameters accessible via ArviZ