BayesianRegressionModels.jl

Appearance

BayesianRegressionModels.jlA brms-shaped formula DSL for Julia

Same formula grammar, but your columns do not all have to come from one equal-length data frame. Runs directly in Turing or transpiles to Stan through StanBlocks.

Gallery
API
GitHub

BayesianRegressionModels.jl ​

A formula DSL for Bayesian regression. The macro @brm parses brms-style syntax (y ~ 1 + a + (1 | g), log(err) ~ 1 + b, y ~ Normal(loc, err)) into a BRMI intermediate representation, then lowers into a backend-specific executor:

The BRM feature atlas gives every executable example the same build-generated four-way view: BRM authoring, the emitted StanBlocks model, generated Stan, and the selected Turing model.

The Warfarin PK/PD examples render both the faithful public two-stage workflow and a joint one-posterior model where shared latent PK effects feed the PK and PD likelihoods.

The Adaptive centering case studies compare fixed coordinates, post-hoc selection and online adaptation, with explicit gradient costs and common scientific quantities:

  • Motorcycle HSGP: two Gaussian processes for the mean and changing noise level.

  • Eight schools: PosteriorDB priors and manual fully centered selection.

  • Radon: county intercepts and floor slopes, with categorical predictive checks.

  • Pupil: numeric scale predictor: automatic total coefficients, brms sum-to-zero comparisons and both adaptation losses.

  • Pupil: hierarchical residual SD: automatic totals in both the mean and residual-scale predictors.

  • Air pollution: regional intercepts and slopes, exact marginalization and matched scientific quantities.

brm-comparison
BRM regression model
julia
intro_model = (@brm begin
    y ~ Normal(loc, err)
    loc ~ 1 + age + sex + (1 + age | subj)
    err ~ Exponential(1)
end)((;
    age=[21.0, 38.0, 55.0, 29.0, 47.0, 61.0],
    sex=[1, 2, 1, 2, 1, 2],
    subj=[1, 1, 2, 2, 3, 3],
    y=[0.2, 1.1, -0.4, 0.7, 1.4, 1.0],
))
julia
BRMI:
  loc ~ 1 + age + sex + ((1 + age) | subj)
  err ~ Exponential(1)
  y ~ Normal(loc, err)
  age: data (eltype=Float64, n=6)
  sex: data (eltype=Int64, n=6)
  subj: data (eltype=Int64, n=6)
julia
SBBRMI with data keys = [:age, :n_subj, :n_terms_loc_subj, :sex, :sex_idx, :sex_n_levels, :subj, :subj_idx, :y]
emitted @slic body:
begin
    err ~ exponential(1.0 ./ 1)
    X_loc = hcat(rep_vector(1.0, num_elements(age)), age)
    pop_loc ~ popefs(; X = X_loc)
    cat_loc_sex ~ _sb_cat(; x = sex_idx, n_levels = sex_n_levels)
    Z_loc_subj = hcat(rep_vector(1.0, num_elements(age)), age)
    r_loc_subj ~ ranef_correlated(; Z = Z_loc_subj, group_idx = subj_idx, n_groups = n_subj, n_terms = n_terms_loc_subj)
    loc = pop_loc + cat_loc_sex + r_loc_subj
    y ~ normal(loc, err)
end
stan
functions {
matrix hcat(
    vector x,
    vector y
) {
    int n = dims(x)[1];
    if (dims(y)[1] != n) reject("hcat: dim mismatch — `y` dim 1 (= ", dims(y)[1], ") does not match `n` (= ", n, "), inferred from `x` dim 1. `n` sizes: `x` dim 1 (= ", dims(x)[1], "), `y` dim 1 (= ", dims(y)[1], ").");
    return append_col(x, y);
}
vector normal_lpdfs(
    vector obs,
    vector loc,
    real scale
) {
    return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
    vector x1,
    vector x2,
    real x3
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = normal_lpdfs(broadcasted_getindex(x1, i), broadcasted_getindex(x2, i), x3);
    }
    return rv;
}
real normal_lpdfs(
    real args1,
    real args2,
    real args3
) {
    return normal_lpdf(args1 | args2, args3);
}
real broadcasted_getindex(vector x, int i) {
    return x[i];
}
vector normal_vector_rng(
    int anontok__1,
    vector a,
    real b
) {
    int n = anontok__1;
    if((n == 0)) {
        vector[n] rv;
        return rv;
    } else {
        return to_vector(normal_rng(a, b));
    }
}
}
data {
    int age_n;
    vector[age_n] age;
    int sex_n_levels;
    int sex_idx_n;
    array[sex_idx_n] int sex_idx;
    int n_terms_loc_subj;
    int n_subj;
    int subj_idx_n;
    array[subj_idx_n] int subj_idx;
    int y_n;
    vector[y_n] y;
}
transformed data {
    matrix[age_n, 2] X_loc = hcat(rep_vector(1.0, num_elements(age)), age);
    int pop_loc_n_covariates = 2;
    matrix[age_n, 2] Z_loc_subj = hcat(rep_vector(1.0, num_elements(age)), age);
}
parameters {
    real<lower=0.0> err;
    vector[pop_loc_n_covariates] pop_loc_beta_pop;
    vector[(sex_n_levels - 1)] cat_loc_sex_beta;
    cholesky_factor_corr[n_terms_loc_subj] r_loc_subj_L;
    vector<lower=0.0>[n_terms_loc_subj] r_loc_subj_tau;
    vector[(n_terms_loc_subj * n_subj)] r_loc_subj_z_flat;
}
transformed parameters {
    vector[age_n] pop_loc = (X_loc * pop_loc_beta_pop);
    vector[sex_idx_n] cat_loc_sex = append_row(0.0, cat_loc_sex_beta)[sex_idx];
    matrix[n_terms_loc_subj, n_subj] r_loc_subj_z = to_matrix(r_loc_subj_z_flat, n_terms_loc_subj, n_subj);
    matrix[n_subj, n_terms_loc_subj] r_loc_subj_b = ((diag_pre_multiply(r_loc_subj_tau, r_loc_subj_L) * r_loc_subj_z)');
    vector[subj_idx_n] r_loc_subj = rows_dot_product(Z_loc_subj, r_loc_subj_b[subj_idx, :]);
    vector[age_n] loc = (pop_loc + cat_loc_sex + r_loc_subj);
}
model {
    err ~ exponential((1.0 ./ 1));
    pop_loc_beta_pop ~ std_normal();
    cat_loc_sex_beta ~ std_normal();
    r_loc_subj_L ~ lkj_corr_cholesky(1.0);
    r_loc_subj_tau ~ std_normal();
    r_loc_subj_z_flat ~ std_normal();
    y ~ normal(loc, err);
}
generated quantities {
    vector[y_n] y_likelihood = normal_lpdfs(y, loc, err);
    vector[y_n] y_gen = normal_vector_rng(y_n, loc, err);
}
julia
#= line 0 =# Turing.@model(function brm_model(y, X_loc, group_matrix_loc_1, group_indices_loc_1, group_levels_loc_1)
        err ~ Distributions.Exponential(1)
        beta_pop ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 3))
        eta_loc = X_loc * beta_pop
        group_effect_1 = Base.zeros(Base.length(y))
        group_1_1 ~ DynamicPPL.to_submodel(BRM.turing_default_correlated_group(group_matrix_loc_1, group_indices_loc_1, group_levels_loc_1, 1.0))
        group_effect_1 = group_effect_1 + group_1_1.effect
        eta_loc = eta_loc + group_effect_1
        loc = eta_loc
        begin
            for i = Base.eachindex(y)
                y[i] ~ Distributions.Normal(loc[i], err)
            end
        end
        (; loc = loc, err = err, response = y)
    end)

Coming from brms ​

BRM is deliberately brms-shaped, so most of what you know transfers. The differences worth knowing up front are one structural gain and a genuinely shorter feature list.

What carries over ​

Formula grammar, random-effect syntax, and the two features people usually reach brms for:

brms@brm
y ~ 1 + age + sexy ~ 1 + age + sex
(1 + age | subj)(1 + age | subj)
(1 + age |p| subj) — correlated across formulas(1 + age |p| subj)
(1 | gr(subj, by = diagnosis))(1 | gr(subj, by = diagnosis))
bf(y1) + bf(y2)two ~ lines in the same block
Distributional regression — bf(y ~ x, sigma ~ x)every distributional parameter is just another ~ line
Nonlinear terms — bf(y ~ a * exp(-b * x), a ~ 1, b ~ 1, nl = TRUE)the same, with no nl switch

Distributional regression needs no special form: give the parameter its own formula, and apply the link yourself where you want one.

brm-comparison
Distributional Gaussian regression
julia
distributional = (@brm begin
    y        ~ Normal(mu, sigma)
    mu       ~ 1 + age
    log(sigma) ~ 1 + age          # explicit link, addressed as `sigma`
end)((; age=[21.0, 38.0, 55.0, 29.0], y=[0.2, 1.1, -0.4, 0.7]))
julia
BRMI:
  mu ~ 1 + age
  log(sigma) ~ 1 + age
  y ~ Normal(mu, sigma)
  age: data (eltype=Float64, n=4)
julia
SBBRMI with data keys = [:age, :y]
emitted @slic body:
begin
    X_mu = hcat(rep_vector(1.0, num_elements(age)), age)
    pop_mu ~ _popefs_coefs(; X = X_mu)
    X_log_sigma = hcat(rep_vector(1.0, num_elements(age)), age)
    pop_log_sigma ~ popefs(; X = X_log_sigma)
    log_sigma = pop_log_sigma
    sigma = exp(log_sigma)
    y ~ normal_id_glm(X_mu, 0.0, pop_mu, sigma)
    mu = X_mu * pop_mu
end
stan
functions {
matrix hcat(
    vector x,
    vector y
) {
    int n = dims(x)[1];
    if (dims(y)[1] != n) reject("hcat: dim mismatch — `y` dim 1 (= ", dims(y)[1], ") does not match `n` (= ", n, "), inferred from `x` dim 1. `n` sizes: `x` dim 1 (= ", dims(x)[1], "), `y` dim 1 (= ", dims(y)[1], ").");
    return append_col(x, y);
}
vector normal_id_glm_lpdfs(
    vector y,
    matrix X,
    real alpha,
    vector beta,
    vector sigma
) {
    int n = dims(y)[1];
    if (dims(X)[1] != n) reject("normal_id_glm_lpdfs: dim mismatch — `X` dim 1 (= ", dims(X)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `X` dim 1 (= ", dims(X)[1], ").");
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = normal_lpdf(y[i] | (alpha + (X[i, :] * beta)), sigma);
    }
    return rv;
}
vector normal_id_glm_vector_rng(
    int anontok__1,
    matrix X,
    real alpha,
    vector beta,
    vector sigma
) {
    int m = anontok__1;
    if (dims(X)[1] != m) reject("normal_id_glm_rng: dim mismatch — `X` dim 1 (= ", dims(X)[1], ") does not match `m` (= ", m, "), inferred from `anontok__1` dim 1. `m` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `X` dim 1 (= ", dims(X)[1], ").");
    if((m == 0)) {
        vector[m] rv;
        return rv;
    } else {
        return normal_id_glm_rng(X, alpha, beta, sigma);
    }
}
vector normal_id_glm_rng(
    matrix X,
    real alpha,
    vector beta,
    vector sigma
) {
    int m = dims(X)[1];
    if((m == 0)) {
        vector[m] rv;
        return rv;
    } else {
        return to_vector(normal_rng((rep_vector(alpha, m) + (X * beta)), sigma));
    }
}
}
data {
    int age_n;
    vector[age_n] age;
    int y_n;
    vector[y_n] y;
}
transformed data {
    matrix[age_n, 2] X_mu = hcat(rep_vector(1.0, num_elements(age)), age);
    int pop_mu_n_covariates = 2;
    matrix[age_n, 2] X_log_sigma = hcat(rep_vector(1.0, num_elements(age)), age);
    int pop_log_sigma_n_covariates = 2;
}
parameters {
    vector[pop_mu_n_covariates] pop_mu_beta_pop;
    vector[pop_log_sigma_n_covariates] pop_log_sigma_beta_pop;
}
transformed parameters {
    vector[pop_mu_n_covariates] pop_mu = pop_mu_beta_pop;
    vector[age_n] pop_log_sigma = (X_log_sigma * pop_log_sigma_beta_pop);
    vector[age_n] log_sigma = pop_log_sigma;
    vector[age_n] sigma = exp(log_sigma);
}
model {
    pop_mu_beta_pop ~ std_normal();
    pop_log_sigma_beta_pop ~ std_normal();
    y ~ normal_id_glm(X_mu, 0.0, pop_mu, sigma);
}
generated quantities {
    vector[age_n] y_likelihood = normal_id_glm_lpdfs(y, X_mu, 0.0, pop_mu, sigma);
    vector[age_n] y_gen = normal_id_glm_vector_rng(y_n, X_mu, 0.0, pop_mu, sigma);
    vector[age_n] mu = (X_mu * pop_mu);
}
julia
#= line 0 =# Turing.@model(function brm_model(y, X_mu, X_sigma)
        beta_pop ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 2))
        eta_mu = X_mu * beta_pop
        mu = eta_mu
        beta_pop_sigma ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 2))
        eta_sigma = X_sigma * beta_pop_sigma
        sigma = Base.exp.(eta_sigma)
        begin
            for i = Base.eachindex(y)
                y[i] ~ Distributions.Normal(mu[i], sigma[i])
            end
        end
        (; mu = mu, sigma = sigma, response = y)
    end)

Nonlinear models need no nl = TRUE and no nlf(). A declaration is a named value, so composing declarations into an arbitrary Julia expression is the whole feature:

brm-comparison
Nonlinear predictor composition
julia
nonlinear = (@brm begin
    a     ~ 1 + (1 | g)           # ordinary linear predictors …
    b     ~ 1
    y     ~ Normal(a * exp(-b * x), sigma)   # … composed nonlinearly
    sigma ~ Exponential(1)
end)((;
    g=[1, 1, 2, 2], x=[0.0, 0.5, 1.0, 1.5],
    y=[1.0, 0.8, 0.5, 0.3],
))
julia
BRMI:
  g: data (eltype=Int64, n=4)
  a ~ 1 + (1 | g)
  b ~ 1
  x: data (eltype=Float64, n=4)
  sigma ~ Exponential(1)
  y ~ Normal((a * exp((-(b) * x))), sigma)
julia
SBBRMI with data keys = [:g, :total_A_a, :total_group_a, :total_location_a, :total_ng_a, :total_nk_a, :total_np_a, :total_precision_a, :x, :y]
configured submodels:
_brm_total_scales_configured_1 = Base.merge(BayesianRegressionModels._brm_total_scales, quote
            tau ~ (ValueFamily(brm_vector_prior_faeb6f6956d0662a))(0.0, 1.0; n = 1)
        end)
emitted @slic body:
begin
    total_scale_a ~ _brm_total_scales_configured_1(; n = total_nk_a)
    total_a::matrix[total_ng_a, total_nk_a] ~ brm_total(total_scale_a, total_A_a, total_location_a, total_precision_a)
    population_a = brm_total_recover_rng(total_a, total_scale_a, total_A_a, total_location_a, total_precision_a)
    deviation_a = brm_total_deviations(total_a, total_A_a * population_a)
    total_Z_a = hcat(rep_vector(1.0, num_elements(total_group_a)))
    a = rows_dot_product(total_a[total_group_a, :], total_Z_a)
    X_b = hcat(rep_vector(1.0, num_elements(y)))
    pop_b ~ popefs(; X = X_b)
    b = pop_b
    sigma ~ exponential(1.0 ./ 1)
    y ~ normal(a .* (exp)((-)(b) .* x), sigma)
end
stan
functions {
// value UDF brm_vector_prior_faeb6f6956d0662a_lpdf
real brm_vector_prior_faeb6f6956d0662a_lpdf(
    vector x,
    real arg_1,
    real arg_2
) {
    if((x[1] < 0.0)) {
        return negative_infinity();
    }
    return lognormal_lpdf(x[1] | arg_1, arg_2);
}
real brm_total_lpdf(
    matrix total,
    vector tau,
    matrix A,
    vector location,
    vector precision
) {
    int j = dims(total)[1];
    int k = dims(total)[2];
    int p = dims(A)[2];
    if (dims(tau)[1] != k) reject("brm_total_lpdf: dim mismatch — `tau` dim 1 (= ", dims(tau)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
    if (dims(A)[1] != k) reject("brm_total_lpdf: dim mismatch — `A` dim 1 (= ", dims(A)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
    if (dims(location)[1] != p) reject("brm_total_lpdf: dim mismatch — `location` dim 1 (= ", dims(location)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], ").");
    if (dims(precision)[1] != p) reject("brm_total_lpdf: dim mismatch — `precision` dim 1 (= ", dims(precision)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], ").");
    matrix[dims(precision)[1], dims(precision)[1]] Q = brm_total_precision(tau, A, precision, j);
    vector[dims(total)[2]] average = brm_total_mean(total);
    vector[dims(precision)[1]] beta = brm_total_conditional_mean(total, tau, A, location, precision, Q);
    vector[dims(total)[2]] residual = (average - (A * beta));
    real quadratic = 0.0;
    real lp = ((-0.5 * ((j * k) - p) * 1.8378770664093453) - (j * sum(log(tau))));
    for(c in 1:k) {
        quadratic += ((j * square(residual[c])) / square(tau[c]));
        for(g in 1:j) {
            quadratic += (square((total[g, c] - average[c])) / square(tau[c]));
        }
    }
    for(a in 1:p) {
        if((precision[a] > 0.0)) {
            lp += (0.5 * (log(precision[a]) - 1.8378770664093453));
            quadratic += (precision[a] * square((beta[a] - location[a])));
        }
    }
    return (lp - (0.5 * (log_determinant(Q) + quadratic)));
}
matrix brm_total_precision(
    vector tau,
    matrix A,
    vector precision,
    int n_groups
) {
    int k = dims(tau)[1];
    int p = dims(A)[2];
    if (dims(A)[1] != k) reject("brm_total_precision: dim mismatch — `A` dim 1 (= ", dims(A)[1], ") does not match `k` (= ", k, "), inferred from `tau` dim 1. `k` sizes: `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
    if (dims(precision)[1] != p) reject("brm_total_precision: dim mismatch — `precision` dim 1 (= ", dims(precision)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `precision` dim 1 (= ", dims(precision)[1], ").");
    matrix[dims(precision)[1], dims(precision)[1]] out = diag_matrix(precision);
    for(a in 1:p) {
        for(b in 1:p) {
            for(c in 1:k) {
                out[a, b] += ((n_groups * A[c, a] * A[c, b]) / square(tau[c]));
            }
        }
    }
    return out;
}
vector brm_total_mean(
    matrix total
) {
    int j = dims(total)[1];
    int k = dims(total)[2];
    vector[k] out = rep_vector(0.0, k);
    for(c in 1:k) {
        out[c] = (sum(total[:, c]) / j);
    }
    return out;
}
vector brm_total_conditional_mean(
    matrix total,
    vector tau,
    matrix A,
    vector location,
    vector precision,
    matrix Q
) {
    int j = dims(total)[1];
    int k = dims(total)[2];
    int p = dims(A)[2];
    if (dims(tau)[1] != k) reject("brm_total_conditional_mean: dim mismatch — `tau` dim 1 (= ", dims(tau)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
    if (dims(A)[1] != k) reject("brm_total_conditional_mean: dim mismatch — `A` dim 1 (= ", dims(A)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
    if (dims(location)[1] != p) reject("brm_total_conditional_mean: dim mismatch — `location` dim 1 (= ", dims(location)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], "), `Q` dim 1 (= ", dims(Q)[1], "), `Q` dim 2 (= ", dims(Q)[2], ").");
    if (dims(precision)[1] != p) reject("brm_total_conditional_mean: dim mismatch — `precision` dim 1 (= ", dims(precision)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], "), `Q` dim 1 (= ", dims(Q)[1], "), `Q` dim 2 (= ", dims(Q)[2], ").");
    if (dims(Q)[1] != p) reject("brm_total_conditional_mean: dim mismatch — `Q` dim 1 (= ", dims(Q)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], "), `Q` dim 1 (= ", dims(Q)[1], "), `Q` dim 2 (= ", dims(Q)[2], ").");
    if (dims(Q)[2] != p) reject("brm_total_conditional_mean: dim mismatch — `Q` dim 2 (= ", dims(Q)[2], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], "), `Q` dim 1 (= ", dims(Q)[1], "), `Q` dim 2 (= ", dims(Q)[2], ").");
    vector[dims(total)[2]] average = brm_total_mean(total);
    vector[dims(location)[1]] natural = (precision .* location);
    for(a in 1:p) {
        for(c in 1:k) {
            natural[a] += ((j * A[c, a] * average[c]) / square(tau[c]));
        }
    }
    return mdivide_left_spd(Q, natural);
}
vector brm_total_recover_rng(
    matrix total,
    vector tau,
    matrix A,
    vector location,
    vector precision
) {
    int j = dims(total)[1];
    int k = dims(total)[2];
    int p = dims(A)[2];
    if (dims(tau)[1] != k) reject("brm_total_recover_rng: dim mismatch — `tau` dim 1 (= ", dims(tau)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
    if (dims(A)[1] != k) reject("brm_total_recover_rng: dim mismatch — `A` dim 1 (= ", dims(A)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
    if (dims(location)[1] != p) reject("brm_total_recover_rng: dim mismatch — `location` dim 1 (= ", dims(location)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], ").");
    if (dims(precision)[1] != p) reject("brm_total_recover_rng: dim mismatch — `precision` dim 1 (= ", dims(precision)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], ").");
    matrix[dims(precision)[1], dims(precision)[1]] Q = brm_total_precision(tau, A, precision, j);
    vector[dims(precision)[1]] beta = brm_total_conditional_mean(total, tau, A, location, precision, Q);
    return multi_normal_rng(beta, inverse_spd(Q));
}
matrix brm_total_deviations(
    matrix total,
    vector mu
) {
    int j = dims(total)[1];
    int k = dims(total)[2];
    if (dims(mu)[1] != k) reject("brm_total_deviations: dim mismatch — `mu` dim 1 (= ", dims(mu)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `mu` dim 1 (= ", dims(mu)[1], ").");
    matrix[dims(total)[1], dims(total)[2]] out = total;
    for(c in 1:k) {
        out[:, c] = (total[:, c] - rep_vector(mu[c], j));
    }
    return out;
}
matrix hcat(vector x) {
    int n = dims(x)[1];
    return to_matrix(x, n, 1);
}
vector normal_lpdfs(
    vector obs,
    vector loc,
    real scale
) {
    return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
    vector x1,
    vector x2,
    real x3
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = normal_lpdfs(broadcasted_getindex(x1, i), broadcasted_getindex(x2, i), x3);
    }
    return rv;
}
real normal_lpdfs(
    real args1,
    real args2,
    real args3
) {
    return normal_lpdf(args1 | args2, args3);
}
real broadcasted_getindex(vector x, int i) {
    return x[i];
}
vector normal_vector_rng(
    int anontok__1,
    vector a,
    real b
) {
    int n = anontok__1;
    if((n == 0)) {
        vector[n] rv;
        return rv;
    } else {
        return to_vector(normal_rng(a, b));
    }
}
}
data {
    int total_ng_a;
    int total_nk_a;
    int total_A_a_m;
    int total_A_a_n;
    matrix[total_A_a_m, total_A_a_n] total_A_a;
    int total_location_a_n;
    vector[total_location_a_n] total_location_a;
    int total_precision_a_n;
    vector[total_precision_a_n] total_precision_a;
    int total_group_a_n;
    array[total_group_a_n] int total_group_a;
    int y_n;
    vector[y_n] y;
    int x_n;
    vector[x_n] x;
}
transformed data {
    matrix[num_elements(total_group_a), 1] total_Z_a = hcat(rep_vector(1.0, num_elements(total_group_a)));
    matrix[num_elements(y), 1] X_b = hcat(rep_vector(1.0, num_elements(y)));
    int pop_b_n_covariates = 1;
}
parameters {
    vector<lower=0.0>[1] total_scale_a_tau;
    matrix[total_ng_a, total_nk_a] total_a;
    vector[pop_b_n_covariates] pop_b_beta_pop;
    real<lower=0.0> sigma;
}
transformed parameters {
    vector<lower=0.0>[1] total_scale_a = total_scale_a_tau;
    vector[num_elements(total_group_a)] a = rows_dot_product(total_a[total_group_a, :], total_Z_a);
    vector[num_elements(y)] pop_b = (X_b * pop_b_beta_pop);
    vector[num_elements(y)] b = pop_b;
}
model {
    total_scale_a_tau ~ brm_vector_prior_faeb6f6956d0662a(0.0, 1.0);
    total_a ~ brm_total(total_scale_a, total_A_a, total_location_a, total_precision_a);
    pop_b_beta_pop ~ std_normal();
    sigma ~ exponential((1.0 ./ 1));
    y ~ normal((a .* exp(((-b) .* x))), sigma);
}
generated quantities {
    vector[total_precision_a_n] population_a = brm_total_recover_rng(total_a, total_scale_a, total_A_a, total_location_a, total_precision_a);
    matrix[total_ng_a, total_A_a_m] deviation_a = brm_total_deviations(total_a, (total_A_a * population_a));
    vector[y_n] y_likelihood = normal_lpdfs(y, (a .* exp(((-b) .* x))), sigma);
    vector[y_n] y_gen = normal_vector_rng(y_n, (a .* exp(((-b) .* x))), sigma);
}
julia
#= line 0 =# Turing.@model(function brm_model(y, X_a, group_effects_a_1, X_b, x)
        beta_pop ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 1))
        eta_a = X_a * beta_pop
        group_effect_1 = Base.zeros(Base.length(y))
        group_1_1 ~ DynamicPPL.to_submodel(BRM.turing_group_effect(group_effects_a_1, (nothing,), nothing))
        group_effect_1 = group_effect_1 + group_1_1.effect
        eta_a = eta_a + group_effect_1
        a = eta_a
        beta_pop_b ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 1))
        eta_b = X_b * beta_pop_b
        b = eta_b
        sigma ~ Distributions.Exponential(1)
        begin
            for i = Base.eachindex(y)
                y[i] ~ Distributions.Normal(a[i] * Base.exp(-(b[i]) * x[i]), sigma)
            end
        end
        (; a = a, b = b, sigma = sigma, response = y)
    end)

That lowers to y ~ normal(a .* exp(-b .* x), sigma) in the generated Stan, and — as with any BRM response — you also get the pointwise log-likelihood y_likelihood and predictive draws y_gen for free.

Where BRM goes further ​

Your data does not have to be one data frame with equal-length columns. This is the main structural difference. brms takes a single data.frame, so every column shares one row axis. @brm takes any column collection — a NamedTuple is fine — whose columns may live on different row axes, and different linear predictors in one model may be defined on different axes. The multi-axis population PK kernel is a runnable example: its subject columns have one row per person, while its observation columns have one row per sample.

ragged(x, group) attaches a flat secondary frame to the grouping axis, and kernel broadcasts a do-block cell — arbitrary Julia, including a structural or ODE-like time course — over pre-grouped rows, taking the per-group values of ordinary linear predictors as arguments. In brms this class of model is nlf() plus manual data2 bookkeeping, or Stan by hand.

Two smaller gains: VBRMI gives you a pure-Julia LogDensityProblems object with no Stan toolchain in the loop, and brm_descriptor exposes one reflectable description of everything the model emits, so a consumer mounts a fitted model without keeping a parallel registry of parameter names.

Where BRM is behind ​

Not feature-complete against brms. The gaps a brms user is most likely to hit:

  • Residual correlation is Gaussian and complete-row only. Use [y1, y2] ~ MvNormalCholesky([mu1, mu2], L_res) with a declared LKJCovarianceFactor; partially missing response vectors and non-Gaussian residual copulas are not implemented.

  • No fixed / known covariance groups. by= is the only gr option — no cov=, so no phylogenetic or pedigree random effects.

  • Autocorrelation is first-order only. ar(time; p=1) emits an ordinary AR(1) noise column; dar(time; p=1) emits a direct differenced-AR(1) trajectory. There is no MA, ARMA, compound symmetry, unstructured, CAR or SAR surface.

  • Fewer families. Likelihoods is the complete list, and it is considerably shorter than brms'.

Formula terms has the full catalogue of what is supported.

Configuring priors ​

Population coefficients use independent standard-normal priors by default. Override selected coefficients with separate effect(...) statements; the coefficient names are exactly those returned by popcoefnames:

brm-comparison
Population effect priors
julia
pk = (@brm begin
    log_ka ~ 1 + weight + (1 | pk | subject)
    effect(log_ka, Intercept) ~ Normal(log(1 / 8), 0.8)
    effect(log_ka, weight) ~ Normal(0, 0.1)
    y ~ Normal(log_ka, 1.0)
end)((;
    weight=[55.0, 65.0, 75.0, 85.0],
    subject=[1, 1, 2, 2], y=[-2.1, -1.9, -1.8, -1.7],
))
julia
BRMI:
  weight: data (eltype=Float64, n=4)
  subject: data (eltype=Int64, n=4)
  log_ka ~ 1 + weight + (1 | pk | subject)
  effect(log_ka, Intercept) ~ Normal(log(/(1, 8)), 0.8)
  effect(log_ka, weight) ~ Normal(0, 0.1)
  y ~ Normal(log_ka, 1.0)
julia
SBBRMI with data keys = [:subject, :total_A_log_ka, :total_group_log_ka, :total_location_log_ka, :total_ng_log_ka, :total_nk_log_ka, :total_np_log_ka, :total_precision_log_ka, :weight, :y]
configured submodels:
_brm_total_scales_configured_1 = Base.merge(BayesianRegressionModels._brm_total_scales, quote
            tau ~ (ValueFamily(brm_vector_prior_ca4b8a1c1bc116d6))(0.0, 1.0; n = 1)
        end)
emitted @slic body:
begin
    total_scale_log_ka ~ _brm_total_scales_configured_1(; n = total_nk_log_ka)
    total_log_ka::matrix[total_ng_log_ka, total_nk_log_ka] ~ brm_total(total_scale_log_ka, total_A_log_ka, total_location_log_ka, total_precision_log_ka)
    population_log_ka = brm_total_recover_rng(total_log_ka, total_scale_log_ka, total_A_log_ka, total_location_log_ka, total_precision_log_ka)
    deviation_log_ka = brm_total_deviations(total_log_ka, total_A_log_ka * population_log_ka)
    total_Z_log_ka = hcat(rep_vector(1.0, num_elements(total_group_log_ka)))
    X_log_ka = hcat(weight)
    pop_log_ka ~ _popefs_normal(; X = X_log_ka, beta_loc = [0], beta_scale = [0.1])
    log_ka = rows_dot_product(total_log_ka[total_group_log_ka, :], total_Z_log_ka) + pop_log_ka
    y ~ normal(log_ka, 1.0)
end
stan
functions {
// value UDF brm_vector_prior_ca4b8a1c1bc116d6_lpdf
real brm_vector_prior_ca4b8a1c1bc116d6_lpdf(
    vector x,
    real arg_1,
    real arg_2
) {
    if((x[1] < 0.0)) {
        return negative_infinity();
    }
    return normal_lpdf(x[1] | arg_1, arg_2);
}
real brm_total_lpdf(
    matrix total,
    vector tau,
    matrix A,
    vector location,
    vector precision
) {
    int j = dims(total)[1];
    int k = dims(total)[2];
    int p = dims(A)[2];
    if (dims(tau)[1] != k) reject("brm_total_lpdf: dim mismatch — `tau` dim 1 (= ", dims(tau)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
    if (dims(A)[1] != k) reject("brm_total_lpdf: dim mismatch — `A` dim 1 (= ", dims(A)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
    if (dims(location)[1] != p) reject("brm_total_lpdf: dim mismatch — `location` dim 1 (= ", dims(location)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], ").");
    if (dims(precision)[1] != p) reject("brm_total_lpdf: dim mismatch — `precision` dim 1 (= ", dims(precision)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], ").");
    matrix[dims(precision)[1], dims(precision)[1]] Q = brm_total_precision(tau, A, precision, j);
    vector[dims(total)[2]] average = brm_total_mean(total);
    vector[dims(precision)[1]] beta = brm_total_conditional_mean(total, tau, A, location, precision, Q);
    vector[dims(total)[2]] residual = (average - (A * beta));
    real quadratic = 0.0;
    real lp = ((-0.5 * ((j * k) - p) * 1.8378770664093453) - (j * sum(log(tau))));
    for(c in 1:k) {
        quadratic += ((j * square(residual[c])) / square(tau[c]));
        for(g in 1:j) {
            quadratic += (square((total[g, c] - average[c])) / square(tau[c]));
        }
    }
    for(a in 1:p) {
        if((precision[a] > 0.0)) {
            lp += (0.5 * (log(precision[a]) - 1.8378770664093453));
            quadratic += (precision[a] * square((beta[a] - location[a])));
        }
    }
    return (lp - (0.5 * (log_determinant(Q) + quadratic)));
}
matrix brm_total_precision(
    vector tau,
    matrix A,
    vector precision,
    int n_groups
) {
    int k = dims(tau)[1];
    int p = dims(A)[2];
    if (dims(A)[1] != k) reject("brm_total_precision: dim mismatch — `A` dim 1 (= ", dims(A)[1], ") does not match `k` (= ", k, "), inferred from `tau` dim 1. `k` sizes: `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
    if (dims(precision)[1] != p) reject("brm_total_precision: dim mismatch — `precision` dim 1 (= ", dims(precision)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `precision` dim 1 (= ", dims(precision)[1], ").");
    matrix[dims(precision)[1], dims(precision)[1]] out = diag_matrix(precision);
    for(a in 1:p) {
        for(b in 1:p) {
            for(c in 1:k) {
                out[a, b] += ((n_groups * A[c, a] * A[c, b]) / square(tau[c]));
            }
        }
    }
    return out;
}
vector brm_total_mean(
    matrix total
) {
    int j = dims(total)[1];
    int k = dims(total)[2];
    vector[k] out = rep_vector(0.0, k);
    for(c in 1:k) {
        out[c] = (sum(total[:, c]) / j);
    }
    return out;
}
vector brm_total_conditional_mean(
    matrix total,
    vector tau,
    matrix A,
    vector location,
    vector precision,
    matrix Q
) {
    int j = dims(total)[1];
    int k = dims(total)[2];
    int p = dims(A)[2];
    if (dims(tau)[1] != k) reject("brm_total_conditional_mean: dim mismatch — `tau` dim 1 (= ", dims(tau)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
    if (dims(A)[1] != k) reject("brm_total_conditional_mean: dim mismatch — `A` dim 1 (= ", dims(A)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
    if (dims(location)[1] != p) reject("brm_total_conditional_mean: dim mismatch — `location` dim 1 (= ", dims(location)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], "), `Q` dim 1 (= ", dims(Q)[1], "), `Q` dim 2 (= ", dims(Q)[2], ").");
    if (dims(precision)[1] != p) reject("brm_total_conditional_mean: dim mismatch — `precision` dim 1 (= ", dims(precision)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], "), `Q` dim 1 (= ", dims(Q)[1], "), `Q` dim 2 (= ", dims(Q)[2], ").");
    if (dims(Q)[1] != p) reject("brm_total_conditional_mean: dim mismatch — `Q` dim 1 (= ", dims(Q)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], "), `Q` dim 1 (= ", dims(Q)[1], "), `Q` dim 2 (= ", dims(Q)[2], ").");
    if (dims(Q)[2] != p) reject("brm_total_conditional_mean: dim mismatch — `Q` dim 2 (= ", dims(Q)[2], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], "), `Q` dim 1 (= ", dims(Q)[1], "), `Q` dim 2 (= ", dims(Q)[2], ").");
    vector[dims(total)[2]] average = brm_total_mean(total);
    vector[dims(location)[1]] natural = (precision .* location);
    for(a in 1:p) {
        for(c in 1:k) {
            natural[a] += ((j * A[c, a] * average[c]) / square(tau[c]));
        }
    }
    return mdivide_left_spd(Q, natural);
}
vector brm_total_recover_rng(
    matrix total,
    vector tau,
    matrix A,
    vector location,
    vector precision
) {
    int j = dims(total)[1];
    int k = dims(total)[2];
    int p = dims(A)[2];
    if (dims(tau)[1] != k) reject("brm_total_recover_rng: dim mismatch — `tau` dim 1 (= ", dims(tau)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
    if (dims(A)[1] != k) reject("brm_total_recover_rng: dim mismatch — `A` dim 1 (= ", dims(A)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
    if (dims(location)[1] != p) reject("brm_total_recover_rng: dim mismatch — `location` dim 1 (= ", dims(location)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], ").");
    if (dims(precision)[1] != p) reject("brm_total_recover_rng: dim mismatch — `precision` dim 1 (= ", dims(precision)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], ").");
    matrix[dims(precision)[1], dims(precision)[1]] Q = brm_total_precision(tau, A, precision, j);
    vector[dims(precision)[1]] beta = brm_total_conditional_mean(total, tau, A, location, precision, Q);
    return multi_normal_rng(beta, inverse_spd(Q));
}
matrix brm_total_deviations(
    matrix total,
    vector mu
) {
    int j = dims(total)[1];
    int k = dims(total)[2];
    if (dims(mu)[1] != k) reject("brm_total_deviations: dim mismatch — `mu` dim 1 (= ", dims(mu)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `mu` dim 1 (= ", dims(mu)[1], ").");
    matrix[dims(total)[1], dims(total)[2]] out = total;
    for(c in 1:k) {
        out[:, c] = (total[:, c] - rep_vector(mu[c], j));
    }
    return out;
}
matrix hcat(vector x) {
    int n = dims(x)[1];
    return to_matrix(x, n, 1);
}
vector normal_lpdfs(
    vector obs,
    vector loc,
    real scale
) {
    return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
    vector x1,
    vector x2,
    real x3
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = normal_lpdfs(broadcasted_getindex(x1, i), broadcasted_getindex(x2, i), x3);
    }
    return rv;
}
real normal_lpdfs(
    real args1,
    real args2,
    real args3
) {
    return normal_lpdf(args1 | args2, args3);
}
real broadcasted_getindex(vector x, int i) {
    return x[i];
}
vector normal_vector_rng(
    int anontok__1,
    vector a,
    real b
) {
    int n = anontok__1;
    if((n == 0)) {
        vector[n] rv;
        return rv;
    } else {
        return to_vector(normal_rng(a, b));
    }
}
}
data {
    int total_ng_log_ka;
    int total_nk_log_ka;
    int total_A_log_ka_m;
    int total_A_log_ka_n;
    matrix[total_A_log_ka_m, total_A_log_ka_n] total_A_log_ka;
    int total_location_log_ka_n;
    vector[total_location_log_ka_n] total_location_log_ka;
    int total_precision_log_ka_n;
    vector[total_precision_log_ka_n] total_precision_log_ka;
    int total_group_log_ka_n;
    array[total_group_log_ka_n] int total_group_log_ka;
    int weight_n;
    vector[weight_n] weight;
    int y_n;
    vector[y_n] y;
}
transformed data {
    matrix[num_elements(total_group_log_ka), 1] total_Z_log_ka = hcat(rep_vector(1.0, num_elements(total_group_log_ka)));
    matrix[weight_n, 1] X_log_ka = hcat(weight);
    int pop_log_ka_n_covariates = 1;
}
parameters {
    vector<lower=0.0>[1] total_scale_log_ka_tau;
    matrix[total_ng_log_ka, total_nk_log_ka] total_log_ka;
    vector[pop_log_ka_n_covariates] pop_log_ka_beta_pop;
}
transformed parameters {
    vector<lower=0.0>[1] total_scale_log_ka = total_scale_log_ka_tau;
    vector[weight_n] pop_log_ka = (X_log_ka * pop_log_ka_beta_pop);
    vector[num_elements(total_group_log_ka)] log_ka = (rows_dot_product(total_log_ka[total_group_log_ka, :], total_Z_log_ka) + pop_log_ka);
}
model {
    total_scale_log_ka_tau ~ brm_vector_prior_ca4b8a1c1bc116d6(0.0, 1.0);
    total_log_ka ~ brm_total(total_scale_log_ka, total_A_log_ka, total_location_log_ka, total_precision_log_ka);
    pop_log_ka_beta_pop ~ normal([0]', [0.1]');
    y ~ normal(log_ka, 1.0);
}
generated quantities {
    vector[total_precision_log_ka_n] population_log_ka = brm_total_recover_rng(
        total_log_ka,
        total_scale_log_ka,
        total_A_log_ka,
        total_location_log_ka,
        total_precision_log_ka
    );
    matrix[total_ng_log_ka, total_A_log_ka_m] deviation_log_ka = brm_total_deviations(total_log_ka, (total_A_log_ka * population_log_ka));
    vector[y_n] y_likelihood = normal_lpdfs(y, log_ka, 1.0);
    vector[y_n] y_gen = normal_vector_rng(y_n, log_ka, 1.0);
}
julia
#= line 0 =# Turing.@model(function brm_model(y, X_log_ka, group_effects_log_ka_1)
        beta_pop ~ Distributions.product_distribution([Distributions.Normal(Base.log(1 / 8), 0.8), Distributions.Normal(0, 0.1)])
        eta_log_ka = X_log_ka * beta_pop
        group_effect_1 = Base.zeros(Base.length(y))
        group_1_1 ~ DynamicPPL.to_submodel(BRM.turing_group_effect(group_effects_log_ka_1, (nothing,), nothing))
        group_effect_1 = group_effect_1 + group_1_1.effect
        eta_log_ka = eta_log_ka + group_effect_1
        log_ka = eta_log_ka
        begin
            for i = Base.eachindex(y)
                y[i] ~ Distributions.Normal(log_ka[i], 1.0)
            end
        end
        (; log_ka = log_ka, response = y)
    end)

Both slots are mandatory, and : is the wildcard for either. effect(:, weight) ~ Normal(0, 0.1) is the default layer for :weight — it reaches that column in every predictor owning it, and a more specific address such as effect(log_ka, weight) overrides it. Two addresses of equal specificity reaching one parameter, and unknown addresses, error. The first shipped lowering supports Normal(location, scale) overrides and retains the existing pop_<predictor>_beta_pop vector parameter, its popcoefnames labels, and descriptor provenance. Inspect the captured formula statements with effect_priors(brmi). This surface belongs to SBBRMI; VBRMI does not implement it.

Categorical contrasts ​

A categorical predictor — a bare integer/CategoricalArray column, or one wrapped in factor(...) — is not a beta_pop column, so popcoefnames deliberately never lists it: it owns a separate cat_<predictor>_<column>_beta vector holding its K−1 treatment contrasts, with the reference level pinned at 0. Those contrasts also default to std_normal(), and the same effect(...) address changes them — keyed by the column name, not the emitted predictor-qualified cat_<predictor>_<column> parameter name:

brm-comparison
Categorical effect prior
julia
categorical_prior = (@brm begin
    sigma ~ Exponential(1)
    mu ~ 1 + factor(g) + x
    effect(mu, g) ~ Normal(0.0, 0.5)   # ⇒ cat_mu_g_beta ~ normal(0.0, 0.5);
    y ~ Normal(mu, sigma)
end)((;
    g=[1, 2, 3, 1, 2, 3], x=[-1.0, -0.5, 0.0, 0.5, 1.0, 1.5],
    y=[-2.4, -2.2, -2.0, -1.8, -1.7, -1.5],
))
julia
BRMI:
  sigma ~ Exponential(1)
  g: data (eltype=Int64, n=6)
  x: data (eltype=Float64, n=6)
  mu ~ 1 + factor(g) + x
  effect(mu, g) ~ Normal(0.0, 0.5)
  y ~ Normal(mu, sigma)
julia
SBBRMI with data keys = [:g, :g_idx, :g_n_levels, :x, :y]
emitted @slic body:
begin
    sigma ~ exponential(1.0 ./ 1)
    X_mu = hcat(rep_vector(1.0, num_elements(x)), x)
    pop_mu ~ popefs(; X = X_mu)
    cat_mu_g ~ _sb_cat_normal(; x = g_idx, n_levels = g_n_levels, beta_loc = 0.0, beta_scale = 0.5)
    mu = pop_mu + cat_mu_g
    y ~ normal(mu, sigma)
end
stan
functions {
matrix hcat(
    vector x,
    vector y
) {
    int n = dims(x)[1];
    if (dims(y)[1] != n) reject("hcat: dim mismatch — `y` dim 1 (= ", dims(y)[1], ") does not match `n` (= ", n, "), inferred from `x` dim 1. `n` sizes: `x` dim 1 (= ", dims(x)[1], "), `y` dim 1 (= ", dims(y)[1], ").");
    return append_col(x, y);
}
vector normal_lpdfs(
    vector obs,
    vector loc,
    real scale
) {
    return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
    vector x1,
    vector x2,
    real x3
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = normal_lpdfs(broadcasted_getindex(x1, i), broadcasted_getindex(x2, i), x3);
    }
    return rv;
}
real normal_lpdfs(
    real args1,
    real args2,
    real args3
) {
    return normal_lpdf(args1 | args2, args3);
}
real broadcasted_getindex(vector x, int i) {
    return x[i];
}
vector normal_vector_rng(
    int anontok__1,
    vector a,
    real b
) {
    int n = anontok__1;
    if((n == 0)) {
        vector[n] rv;
        return rv;
    } else {
        return to_vector(normal_rng(a, b));
    }
}
}
data {
    int x_n;
    vector[x_n] x;
    int g_n_levels;
    int g_idx_n;
    array[g_idx_n] int g_idx;
    int y_n;
    vector[y_n] y;
}
transformed data {
    matrix[x_n, 2] X_mu = hcat(rep_vector(1.0, num_elements(x)), x);
    int pop_mu_n_covariates = 2;
}
parameters {
    real<lower=0.0> sigma;
    vector[pop_mu_n_covariates] pop_mu_beta_pop;
    vector[(g_n_levels - 1)] cat_mu_g_beta;
}
transformed parameters {
    vector[x_n] pop_mu = (X_mu * pop_mu_beta_pop);
    vector[g_idx_n] cat_mu_g = append_row(0.0, cat_mu_g_beta)[g_idx];
    vector[x_n] mu = (pop_mu + cat_mu_g);
}
model {
    sigma ~ exponential((1.0 ./ 1));
    pop_mu_beta_pop ~ std_normal();
    cat_mu_g_beta ~ normal(0.0, 0.5);
    y ~ normal(mu, sigma);
}
generated quantities {
    vector[y_n] y_likelihood = normal_lpdfs(y, mu, sigma);
    vector[y_n] y_gen = normal_vector_rng(y_n, mu, sigma);
}
julia
#= line 0 =# Turing.@model(function brm_model(y, X_mu)
        sigma ~ Distributions.Exponential(1)
        beta_pop ~ Distributions.product_distribution([Distributions.Normal(0, 1), Distributions.Normal(0.0, 0.5), Distributions.Normal(0.0, 0.5), Distributions.Normal(0, 1)])
        eta_mu = X_mu * beta_pop
        mu = eta_mu
        begin
            for i = Base.eachindex(y)
                y[i] ~ Distributions.Normal(mu[i], sigma)
            end
        end
        (; mu = mu, sigma = sigma, response = y)
    end)

One shared (location, scale) covers every contrast in the block; a treatment contrast has no per-level address (a cell mean does — see below). The :-predictor form effect(:, g) reaches the corresponding block in every predictor owning it, and the statement composes with population overrides on the same predictor (effect(mu, x) ~ Normal(0, 0.25)) — each addresses its own parameter. A non-default reference level emits cat_<predictor>_<column>__ref_<k>_beta, which the plain column name still addresses whenever that is unambiguous; when two factor(g; ref=…) blocks of one column would both claim it, the bare address is refused and each block is addressed by its reference-qualified column name. Models with no such statement keep the same std_normal() contrast prior under the predictor-qualified name.

Cell means: a categorical predictor without an intercept ​

Treatment contrasts measure each level against a reference, so they need an intercept to measure from. A predictor without one codes its first categorical term by cell means instead — one coefficient per level, no reference level, exactly what 0 + factor(g) means in brms and in R's model.matrix. Each cell mean is addressable on its own, as <column>_lvl_<k> with k the level's position in the fitted level order, so a per-group location can carry a per-group prior:

brm-comparison
Cell means with a per-level prior
julia
cell_means = (@brm begin
    sigma ~ Exponential(1)
    mu ~ 0 + factor(site)
    effect(mu, site) ~ Normal(0.0, 2.0)          # every site
    effect(mu, site_lvl_3) ~ Normal(4.0, 0.5)    # ... except the third
    y ~ Normal(mu, sigma)
end)((;
    site=[1, 2, 3, 1, 2, 3],
    y=[-0.4, 0.2, 4.1, -0.1, 0.5, 3.8],
))
julia
BRMI:
  sigma ~ Exponential(1)
  site: data (eltype=Int64, n=6)
  mu ~ 0 + factor(site)
  effect(mu, site) ~ Normal(0.0, 2.0)
  effect(mu, site_lvl_3) ~ Normal(4.0, 0.5)
  y ~ Normal(mu, sigma)
julia
SBBRMI with data keys = [:site, :site_idx, :site_n_levels, :y]
emitted @slic body:
begin
    sigma ~ exponential(1.0 ./ 1)
    cat_mu_site ~ _sb_cat_cells_normal(; x = site_idx, n_levels = site_n_levels, beta_loc = [0.0, 0.0, 4.0], beta_scale = [2.0, 2.0, 0.5])
    mu = cat_mu_site
    y ~ normal(mu, sigma)
end
stan
functions {
vector normal_lpdfs(
    vector obs,
    vector loc,
    real scale
) {
    return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
    vector x1,
    vector x2,
    real x3
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = normal_lpdfs(broadcasted_getindex(x1, i), broadcasted_getindex(x2, i), x3);
    }
    return rv;
}
real normal_lpdfs(
    real args1,
    real args2,
    real args3
) {
    return normal_lpdf(args1 | args2, args3);
}
real broadcasted_getindex(vector x, int i) {
    return x[i];
}
vector normal_vector_rng(
    int anontok__1,
    vector a,
    real b
) {
    int n = anontok__1;
    if((n == 0)) {
        vector[n] rv;
        return rv;
    } else {
        return to_vector(normal_rng(a, b));
    }
}
}
data {
    int site_n_levels;
    int site_idx_n;
    array[site_idx_n] int site_idx;
    int y_n;
    vector[y_n] y;
}
transformed data {
}
parameters {
    real<lower=0.0> sigma;
    vector[site_n_levels] cat_mu_site_beta;
}
transformed parameters {
    vector[site_idx_n] cat_mu_site = cat_mu_site_beta[site_idx];
    vector[site_idx_n] mu = cat_mu_site;
}
model {
    sigma ~ exponential((1.0 ./ 1));
    cat_mu_site_beta ~ normal([0.0, 0.0, 4.0]', [2.0, 2.0, 0.5]');
    y ~ normal(mu, sigma);
}
generated quantities {
    vector[y_n] y_likelihood = normal_lpdfs(y, mu, sigma);
    vector[y_n] y_gen = normal_vector_rng(y_n, mu, sigma);
}
julia
#= line 0 =# Turing.@model(function brm_model(y, X_mu)
        sigma ~ Distributions.Exponential(1)
        beta_pop ~ Distributions.product_distribution([Distributions.Normal(0.0, 2.0), Distributions.Normal(0.0, 2.0), Distributions.Normal(4.0, 0.5)])
        eta_mu = X_mu * beta_pop
        mu = eta_mu
        begin
            for i = Base.eachindex(y)
                y[i] ~ Distributions.Normal(mu[i], sigma)
            end
        end
        (; mu = mu, sigma = sigma, response = y)
    end)

The emitted carrier is the same cat_<predictor>_<column>_beta vector, of length K rather than K−1 and with no pinned zero. The column address sets one prior over all K cell means; a level address is more specific and overrides it for that level, under the same most-specific-wins rule as every other effect(...) statement.

The rule, in full:

  • BRM has no implicit intercept and 0 is only a marker, so "without an intercept" means the predictor has no 1 term: mu ~ site, mu ~ 0 + site and mu ~ 0 + factor(site) are the same formula.

  • Only the first categorical term is cell-mean coded. A second one's full indicator set would be collinear with the first's, so later categorical terms — and every & interaction — stay treatment-coded.

  • factor(site; cmc=false) — brms' switch, "cell-mean coding" — keeps K−1 treatment contrasts in a predictor without an intercept, and the cell means pass to the next categorical term. A ref= alone does not opt out: as in R, a releveled factor without an intercept is still cell-mean coded, in its releveled order.

  • An ordinal model's estimated thresholds are its location predictor's intercept, so the eta of Ordinal(...) / OrderedLogistic(...) keeps treatment contrasts even though it is written eta ~ 0 + ....

  • A random intercept (1 | g) is not a population intercept.

  • The same rule holds inside a random-effect block, decided on that block's own left-hand side: (0 + c | g) gives every level of c its own group-level effect (margins c_dummy_1 … c_dummy_K, addressable by sd(lp, ID, c_dummy_k) on a shared |ID| block), (1 + c | g) keeps the random intercept plus K−1 dummies, and (0 + factor(c; cmc=false) | g) opts out. BRM merges (1 | g) + (0 + c | g) into one block, so the sibling intercept keeps c treatment-coded there; (0 + c || g) is decided on its original left-hand side before it is split into uncorrelated terms.

brm_population_effect_coordinates reports which coding a block has (coding === :cellmeans or :treatment); a cell-mean block has no reference_level, and cells pairs each level with its posterior coordinate. An r2d2(...) decomposition allocates its shares over treatment contrasts and refuses a cell-mean block.

Term-internal parameters ​

Some terms own parameters no coefficient address can reach — s(x)'s smoothing scale, mo(c)'s Dirichlet increments, me(x, sd)'s latent true covariate, a Gaussian process's length scale and amplitude. They are addressed by naming the term itself in the target slot, under the same head-position grammar:

brm-comparison
Term-internal priors
julia
term_priors_example = (@brm begin
    y ~ Normal(mu, 1.)
    mu ~ 1 + s(age) + mo(dose) + me(w_obs, 0.3) + hsgp(conc; k=5)

    sd(:, s(age))               ~ Exponential(2)     # smoothing SD
    simplex(mu, mo(dose))       ~ Dirichlet(2)       # monotonic increments
    latent(:, me(w_obs))        ~ Normal(0, 5)       # latent true covariate
    length_scale(:, hsgp(conc)) ~ Uniform(0.84, 2)   # GP length scale
    sd(:, hsgp(conc))           ~ Normal(0, 0.5)     # GP marginal amplitude
end)((;
    age=collect(20.0:29.0), dose=repeat(1:5; inner=2),
    w_obs=collect(1.0:10.0), conc=collect(range(-2, 2; length=10)),
    y=collect(range(-1, 1; length=10)),
))
julia
BRMI:
  mu ~ 1 + s(age) + mo(dose) + me(w_obs, 0.3) + hsgp(conc; k=5)
  y ~ Normal(mu, 1.0)
  age: data (eltype=Float64, n=10)
  dose: data (eltype=Int64, n=10)
  w_obs: data (eltype=Float64, n=10)
  conc: data (eltype=Float64, n=10)
  effect(term_sd, s(age), :) ~ Exponential(2)
  effect(term_simplex, mo(dose), mu) ~ Dirichlet(2)
  effect(term_latent, me(w_obs), :) ~ Normal(0, 5)
  effect(term_length_scale, hsgp(conc), :) ~ Uniform(0.84, 2)
  effect(term_sd, hsgp(conc), :) ~ Normal(0, 0.5)
julia
SBBRMI with data keys = [:PHI_hsgp_conc, :Xnull_age, :Zpen_age, :age, :conc, :dose, :dose_idx, :omega2_hsgp_conc, :rho_lower_hsgp_conc, :sd_w_obs, :w_obs, :y]
configured submodels:
_sb_s_generic_configured_1 = Base.merge(BayesianRegressionModels._sb_s_generic, quote
            sd_pen ~ (ValueFamily(brm_vector_prior_655a23ff68a31dcc))(0.5; n = 1)
        end)
_sb_hsgp_configured_1 = Base.merge(BayesianRegressionModels._sb_hsgp, quote
            rho_iso ~ uniform(0.84, 2; lower = 0.84, upper = 2.0)
            sigma ~ normal(0.0, 0.5; lower = 0.0)
        end)
emitted @slic body:
begin
    mo_dose ~ _sb_mo(; x = dose_idx, alpha = rep_vector(2.0, 4))
    me_w_obs ~ _sb_me(; x_obs = w_obs, sd_x = sd_w_obs, x_true_loc = 0, x_true_scale = 5)
    X_mu = hcat(rep_vector(1.0, num_elements(dose)), mo_dose, me_w_obs)
    pop_mu ~ popefs(; X = X_mu)
    s_age ~ _sb_s_generic_configured_1(; Xnull = Xnull_age, Zpen = Zpen_age)
    hsgp_conc ~ _sb_hsgp_configured_1(; PHI = PHI_hsgp_conc, omega2 = omega2_hsgp_conc, rho_lower = rho_lower_hsgp_conc)
    mu = pop_mu + s_age + hsgp_conc
    y ~ normal(mu, 1.0)
end
stan
functions {
vector normal_lpdfs(
    vector obs,
    vector loc,
    real scale
) {
    return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
    vector x1,
    vector x2,
    real x3
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = normal_lpdfs(broadcasted_getindex(x1, i), broadcasted_getindex(x2, i), x3);
    }
    return rv;
}
real normal_lpdfs(
    real args1,
    real args2,
    real args3
) {
    return normal_lpdf(args1 | args2, args3);
}
real broadcasted_getindex(vector x, int i) {
    return x[i];
}
vector normal_vector_rng(
    int anontok__1,
    vector a,
    real b
) {
    int n = anontok__1;
    if((n == 0)) {
        vector[n] rv;
        return rv;
    } else {
        return to_vector(normal_rng(a, b));
    }
}
matrix hcat(vector x, vector y, vector z) {
    return hcat(hcat(x, y), z);
}
matrix hcat(
    matrix x,
    vector y
) {
    int m = dims(x)[1];
    int n = dims(x)[2];
    if (dims(y)[1] != m) reject("hcat: dim mismatch — `y` dim 1 (= ", dims(y)[1], ") does not match `m` (= ", m, "), inferred from `x` dim 1. `m` sizes: `x` dim 1 (= ", dims(x)[1], "), `y` dim 1 (= ", dims(y)[1], ").");
    return append_col(x, y);
}
matrix hcat(
    vector x,
    vector y
) {
    int n = dims(x)[1];
    if (dims(y)[1] != n) reject("hcat: dim mismatch — `y` dim 1 (= ", dims(y)[1], ") does not match `n` (= ", n, "), inferred from `x` dim 1. `n` sizes: `x` dim 1 (= ", dims(x)[1], "), `y` dim 1 (= ", dims(y)[1], ").");
    return append_col(x, y);
}
// value UDF brm_vector_prior_655a23ff68a31dcc_lpdf
real brm_vector_prior_655a23ff68a31dcc_lpdf(
    vector x,
    real arg_1
) {
    if((x[1] < 0.0)) {
        return negative_infinity();
    }
    return exponential_lpdf(x[1] | arg_1);
}
vector brm_hsgp_sqrt_spd(
    matrix omega2,
    real sigma,
    vector rho
) {
    int m = dims(omega2)[1];
    int d = dims(omega2)[2];
    if (dims(rho)[1] != d) reject("brm_hsgp_sqrt_spd: dim mismatch — `rho` dim 1 (= ", dims(rho)[1], ") does not match `d` (= ", d, "), inferred from `omega2` dim 2. `d` sizes: `omega2` dim 2 (= ", dims(omega2)[2], "), `rho` dim 1 (= ", dims(rho)[1], ").");
    vector[m] rv;
    real scale = sigma;
    for(axis in 1:d) {
        scale *= sqrt((rho[axis] * 2.5066282746310002));
    }
    for(b in 1:m) {
        real exponent = 0.0;
        for(axis in 1:d) {
            exponent += (rho[axis] * rho[axis] * omega2[b, axis]);
        }
        rv[b] = (scale * exp((-0.25 * exponent)));
    }
    return rv;
}
}
data {
    int dose_idx_n;
    array[dose_idx_n] int dose_idx;
    int w_obs_n;
    vector[w_obs_n] w_obs;
    real sd_w_obs;
    int dose_n;
    array[dose_n] int dose;
    int Zpen_age_n;
    int Xnull_age_m;
    int Xnull_age_n;
    matrix[Xnull_age_m, Xnull_age_n] Xnull_age;
    int Zpen_age_m;
    matrix[Zpen_age_m, Zpen_age_n] Zpen_age;
    int omega2_hsgp_conc_m;
    int omega2_hsgp_conc_n;
    matrix[omega2_hsgp_conc_m, omega2_hsgp_conc_n] omega2_hsgp_conc;
    int PHI_hsgp_conc_m;
    int PHI_hsgp_conc_n;
    matrix[PHI_hsgp_conc_m, PHI_hsgp_conc_n] PHI_hsgp_conc;
    int y_n;
    vector[y_n] y;
}
transformed data {
    int pop_mu_n_covariates = (2 + 1);
    int s_age_n_pen = Zpen_age_n;
    int hsgp_conc_n_basis = omega2_hsgp_conc_m;
    int hsgp_conc_n_axes = omega2_hsgp_conc_n;
}
parameters {
    simplex[4] mo_dose_simplex_incr;
    vector[num_elements(w_obs)] me_w_obs_x_true;
    vector[pop_mu_n_covariates] pop_mu_beta_pop;
    vector[2] s_age_b_fixed;
    vector<lower=0.0>[1] s_age_sd_pen;
    vector[s_age_n_pen] s_age_b_pen_raw;
    real<lower=0.84, upper=2.0> hsgp_conc_rho_iso;
    real<lower=0.0> hsgp_conc_sigma;
    vector[hsgp_conc_n_basis] hsgp_conc_beta_raw;
}
transformed parameters {
    vector[dose_idx_n] mo_dose = cumulative_sum(append_row(0.0, mo_dose_simplex_incr))[dose_idx];
    vector[num_elements(w_obs)] me_w_obs = me_w_obs_x_true;
    matrix[num_elements(w_obs), (2 + 1)] X_mu = hcat(rep_vector(1.0, num_elements(dose)), mo_dose, me_w_obs);
    vector[num_elements(w_obs)] pop_mu = (X_mu * pop_mu_beta_pop);
    vector[s_age_n_pen] s_age_b_pen = (s_age_sd_pen[1] * s_age_b_pen_raw);
    vector[Xnull_age_m] s_age = ((Xnull_age * s_age_b_fixed) + (Zpen_age * s_age_b_pen));
    vector[hsgp_conc_n_axes] hsgp_conc_rho = rep_vector(hsgp_conc_rho_iso, hsgp_conc_n_axes);
    vector[omega2_hsgp_conc_m] hsgp_conc_sqrt_spd = brm_hsgp_sqrt_spd(omega2_hsgp_conc, hsgp_conc_sigma, hsgp_conc_rho);
    vector[PHI_hsgp_conc_m] hsgp_conc = (PHI_hsgp_conc * (hsgp_conc_sqrt_spd .* hsgp_conc_beta_raw));
    vector[num_elements(w_obs)] mu = (pop_mu + s_age + hsgp_conc);
}
model {
    mo_dose_simplex_incr ~ dirichlet(rep_vector(2.0, 4));
    me_w_obs_x_true ~ normal(0, 5);
    w_obs ~ normal(me_w_obs_x_true, sd_w_obs);
    pop_mu_beta_pop ~ std_normal();
    s_age_sd_pen ~ brm_vector_prior_655a23ff68a31dcc(0.5);
    s_age_b_pen_raw ~ std_normal();
    hsgp_conc_rho_iso ~ uniform(0.84, 2);
    hsgp_conc_sigma ~ normal(0.0, 0.5);
    hsgp_conc_beta_raw ~ std_normal();
    y ~ normal(mu, 1.0);
}
generated quantities {
    vector[w_obs_n] w_obs_likelihood = normal_lpdfs(w_obs, me_w_obs_x_true, sd_w_obs);
    vector[w_obs_n] w_obs_gen = normal_vector_rng(w_obs_n, me_w_obs_x_true, sd_w_obs);
    vector[y_n] y_likelihood = normal_lpdfs(y, mu, 1.0);
    vector[y_n] y_gen = normal_vector_rng(y_n, mu, 1.0);
}
julia
#= line 0 =# Turing.@model(function brm_model(y, X_mu, terms_mu_1, age, terms_mu_2, dose, terms_mu_3, w_obs, terms_mu_4, conc)
        beta_pop ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 1))
        eta_mu = X_mu * beta_pop
        term_mu_1 ~ DynamicPPL.to_submodel(BRM.turing_term_model(terms_mu_1, Base.length(y), NamedTuple{$(QuoteNode((:sd,)))}((Distributions.Exponential(2),)), NamedTuple{$(QuoteNode((:age,)))}((age,))))
        eta_mu = eta_mu + term_mu_1.effect
        term_mu_2 ~ DynamicPPL.to_submodel(BRM.turing_term_model(terms_mu_2, Base.length(y), NamedTuple{$(QuoteNode((:simplex,)))}((Distributions.Dirichlet(Base.vect(2.0, 2.0, 2.0, 2.0)),)), NamedTuple{$(QuoteNode((:dose,)))}((dose,))))
        eta_mu = eta_mu + term_mu_2.effect
        term_mu_3 ~ DynamicPPL.to_submodel(BRM.turing_term_model(terms_mu_3, Base.length(y), NamedTuple{$(QuoteNode((:latent,)))}((Distributions.Normal(0, 5),)), NamedTuple{$(QuoteNode((:w_obs,)))}((w_obs,))))
        eta_mu = eta_mu + term_mu_3.effect
        term_mu_4 ~ DynamicPPL.to_submodel(BayesianRegressionModelsTuringExt._brm_turing_hsgp_ncp_model(terms_mu_4, Base.length(y), NamedTuple{$(QuoteNode((:rho, :sigma)))}((Distributions.Uniform(0.84, 2), Distributions.Normal(0, 0.5))), NamedTuple{$(QuoteNode((:conc,)))}((conc,))))
        eta_mu = eta_mu + term_mu_4.effect
        mu = eta_mu
        begin
            for i = Base.eachindex(y)
                y[i] ~ Distributions.Normal(mu[i], 1.0)
            end
        end
        (; mu = mu, response = y)
    end)

The term is spelled the way the formula spells it, minus numeric and keyword arguments — me(w_obs, 0.3) is addressed as me(w_obs). term_priors(brmi) returns the captured statements. See Term-internal priors for the full table, the t2 component slot, why the standardized raw innovations are deliberately not configurable, and the approximation-validity floor an hsgp term puts on its length scale by default.

See Formula terms and Likelihoods for the supported syntax and backend-specific contracts. The Gallery provides live, interactive examples — input formula, the SLIC submodel body, the transpiled Stan source, and the auto-generated posterior-predictive check, all in one card. The API page lists every public binding.