Formula terms
Catalogue
Every term below lowers on the StanBlocks backend (SBBRMI); the VBRMI column marks the ones the pure-Julia backend also implements. The pages that follow document the terms whose behaviour is not obvious from the signature — the rest are covered by their docstrings on the API page.
Predictor terms
| Term | What it contributes | VBRMI |
|---|---|---|
offset(x) | x added with coefficient one — no beta_pop column | ✓ |
s(x) | rank-10 penalized thin-plate regression spline | — |
t2(x, z) | two-margin tensor-product smooth | — |
gp(x…; cov=:exp_quad, iso=true, jitter=1e-9) | exact latent Gaussian process, noncentered Cholesky draw | — |
gp(x; cov=:periodic, period, jitter=1e-9) | exact latent GP under Stan's periodic kernel, one axis | — |
hsgp(x…; k=20, c=1.5, iso=true, by=nothing, domain=nothing, orthogonal_to=nothing, centeredness=0) | Hilbert-space GP approximation over prod(k) basis functions; optional fixed partial centering | — |
hsgp(x; k=20, cov=:periodic, period) | periodic Hilbert-space basis: k harmonics, 2k cosine/sine functions | — |
ar(time; p=1) | AR(p) noise process ordered by time; only p=1 is emitted | — |
dar(time; p=1) | direct differenced-AR(1) trajectory with bounded persistence and scaled innovations | — |
rw(time) | direct random-walk trajectory: dar with the persistence fixed at zero | — |
cdar(step; by=group, cor=C) | per-group deviations following a damped walk over step with innovations correlated across groups by C | — |
mo(c), mo1(c) | monotonic effect of an ordered factor via Dirichlet increments | — |
me(x, sd) | measurement-error covariate — x is observed with known sd | — |
interval_censored(x; upper=lloq, lower=0) | quantified/BLOQ covariate with bounded latent values on BLOQ rows | — |
factor(c), a bare integer / CategoricalArray column | K−1 treatment contrasts under an intercept; K cell means as the first categorical term of a predictor without one ("Cell means" on the overview page) | ✓ |
factor(c; ref=k) | reference level k for the term's treatment contrasts | ✓ |
factor(c; cmc=false) | keep K−1 treatment contrasts even in a predictor without an intercept (brms' cmc, "cell-mean coding") | ✓ |
protect(x) | materialize a raw data expression as one literal column | ✓ |
A plain RHS expression in raw data columns (log(exposure), x^2) is treated as an implicit protect(...) and materialized the same way.
Differenced-AR(1) trajectories
dar(time; p=1) is not an alias for ar(time; p=1). It emits the path
x[1] = 0
d[t] = beta * d[t-1] + sigma * z[t] # d[0] = 0
x[t+1] = x[t] + d[t]as a direct predictor summand. In log_r_week ~ 1 + dar(week), the population intercept is therefore the initial level x0; the trajectory receives no second beta_pop multiplier. The term samples beta on [0, 1], a positive innovation scale sigma, and length(week)-1 standardized innovations z. Its defaults are the CDC wastewater priors beta ~ Normal(0.5, 0.2) truncated to [0, 1] and sigma ~ Normal(0, 0.2) truncated at zero.
The time column must be nonempty, finite, strictly increasing, and unique. Its values establish order; spacing does not rescale the recurrence. Expand a weekly path onto a daily axis explicitly with a deterministic index operation such as log_r = weekly_expand(log_r_week, week_idx). reprocess may replace the ordered grid, including its length, without changing the Stan source.
Configure the two model-scale priors by addressing the term:
ar(:, dar(week)) ~ Normal(0.4, 0.1) # bounded persistence beta
sd(:, dar(week)) ~ Normal(0.0, 0.3) # positive innovation sigmaar accepts Normal, Beta, or an in-bounds Uniform; sd accepts the positive-scale family set documented below. The standardized z innovations are inspectable through the descriptor but deliberately have no prior override.
Random-walk trajectories
rw(time) is dar(time) with the increments' persistence fixed at zero — the path
x[1] = 0
x[t+1] = x[t] + sigma * z[t]over the sorted distinct values of time, as a direct predictor summand: in log_R ~ 1 + rw(time) the population intercept is the initial level and the term samples a positive innovation scale sigma (default Normal(0, 0.2) truncated at zero) and one standardized innovation z per step after the first. Each row reads the point of its own time value, so a long frame whose rows share times (several groups per day) gets one shared walk — with unique times it is the plain path. Address the scale as sd(:, rw(time)); there is no persistence to address, so ar(:, rw(time)) is refused with a message saying so. The time column must be nonempty and finite; on replay the grid may gain new times (a forecast extends the walk with prior innovations). This is the log-reproduction-number walk of a renewal model (Epidemic renewal models), stated as one formula line.
Random-walk trajectoryrandom_walk_term = (@brm begin
mu ~ 1 + rw(t)
effect(mu, Intercept) ~ Normal(0.0, 1.0)
sd(:, rw(t)) ~ Normal(0.0, 0.05)
y ~ Normal(mu, 0.3)
end)((;
t=collect(1.0:8),
y=[0.4, 1.1, 0.9, 1.6, 1.2, 2.0, 1.8, 2.5],
))BRMI:
t: data (eltype=Float64, n=8)
mu ~ 1 + rw(t)
effect(mu, Intercept) ~ Normal(0.0, 1.0)
effect(term_sd, rw(t), :) ~ Normal(0.0, 0.05)
y ~ Normal(mu, 0.3)SBBRMI with data keys = [:rw_mu_t_n_steps, :rw_mu_t_time_idx, :t, :y]
configured submodels:
_sb_rw1_configured_1 = Base.merge(BayesianRegressionModels._sb_rw1, quote
sigma ~ normal(0.0, 0.05; lower = 0.0)
end)
emitted @slic body:
begin
X_mu = hcat(rep_vector(1.0, num_elements(t)))
pop_mu ~ _popefs_normal(; X = X_mu, beta_loc = [0.0], beta_scale = [1.0])
rw_mu_t ~ _sb_rw1_configured_1(; n_steps = rw_mu_t_n_steps, time_idx = rw_mu_t_time_idx)
mu = pop_mu + rw_mu_t
y ~ normal(mu, 0.3)
endfunctions {
matrix hcat(vector x) {
int n = dims(x)[1];
return to_matrix(x, n, 1);
}
vector random_walk_rows(
real sigma,
vector z,
array[] int idx
) {
int N = dims(idx)[1];
vector[(dims(z)[1] + 1)] x = random_walk_path(sigma, z);
vector[N] out;
for(i in 1:N) {
out[i] = x[idx[i]];
}
return out;
}
vector random_walk_path(
real sigma,
vector z
) {
int n = dims(z)[1];
return append_row(0.0, (sigma * cumulative_sum(z)));
}
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 t_n;
vector[t_n] t;
int rw_mu_t_n_steps;
int rw_mu_t_time_idx_n;
array[rw_mu_t_time_idx_n] int rw_mu_t_time_idx;
int y_n;
vector[y_n] y;
}
transformed data {
matrix[num_elements(t), 1] X_mu = hcat(rep_vector(1.0, num_elements(t)));
int pop_mu_n_covariates = 1;
}
parameters {
vector[pop_mu_n_covariates] pop_mu_beta_pop;
real<lower=0.0> rw_mu_t_sigma;
vector[(rw_mu_t_n_steps - 1)] rw_mu_t_z;
}
transformed parameters {
vector[num_elements(t)] pop_mu = (X_mu * pop_mu_beta_pop);
vector[rw_mu_t_time_idx_n] rw_mu_t = random_walk_rows(rw_mu_t_sigma, rw_mu_t_z, rw_mu_t_time_idx);
vector[num_elements(t)] mu = (pop_mu + rw_mu_t);
}
model {
pop_mu_beta_pop ~ normal([0.0]', [1.0]');
rw_mu_t_sigma ~ normal(0.0, 0.05);
rw_mu_t_z ~ std_normal();
y ~ normal(mu, 0.3);
}
generated quantities {
vector[y_n] y_likelihood = normal_lpdfs(y, mu, 0.3);
vector[y_n] y_gen = normal_vector_rng(y_n, mu, 0.3);
}#= line 0 =# Turing.@model(function brm_model(y, X_mu, terms_mu_1, t)
beta_pop ~ Distributions.product_distribution([Distributions.Normal(0.0, 1.0)])
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.Normal(0.0, 0.05),)), NamedTuple{$(QuoteNode((:t,)))}((t,))))
eta_mu = eta_mu + term_mu_1.effect
mu = eta_mu
begin
for i = Base.eachindex(y)
y[i] ~ Distributions.Normal(mu[i], 0.3)
end
end
(; mu = mu, response = y)
end)Grouped correlated damped walks
cdar(step; by=group, cor=C) gives every level of group its own deviation path over the sorted distinct values of step, damped with a shared persistence rho and driven by innovations that are correlated across groups through the Cholesky factor L of C (L L' = C):
delta[:, 1] = sigma * L * eta[:, 1]
delta[:, w] = rho * delta[:, w-1] + sigma * sqrt(1 - rho^2) * L * eta[:, w]Each row contributes delta[group(row), step(row)] as a direct summand, so log_R ~ 1 + rw(time) + cdar(week; by=patch, cor=C) is a shared random walk plus spatially correlated weekly patch deviations — the six-patch model of Epidemic renewal models, stated on the formula surface. C is a P × P symmetric positive-definite matrix, P the number of group levels, supplied as a data field (a matrix-valued field is accepted) or a literal; it is a fixed hyperparameter, not a sampled covariance. The term samples sigma > 0 (default Normal(0, 0.2) truncated at zero), rho on [0, 1] (default Normal(0.5, 0.2) truncated), and P·W standardized innovations eta; address them as sd(:, cdar(step)) and ar(:, cdar(step)). On replay the group levels and L are frozen from the fit while the step grid may grow — a forecast extends the walk with prior innovations — and an unseen group level is refused.
Grouped correlated damped walkgrouped_walk = (@brm begin
mu ~ 1 + cdar(week; by=patch, cor=C)
sd(:, cdar(week)) ~ Normal(0.0, 0.2)
ar(:, cdar(week)) ~ Normal(0.8, 0.1)
y ~ Normal(mu, 0.3)
end)((;
week=[1, 1, 2, 2, 3, 3],
patch=["a", "b", "a", "b", "a", "b"],
y=[0.2, -0.1, 0.4, 0.0, 0.5, 0.1],
C=[1.0 0.6; 0.6 1.0],
))BRMI:
patch: data (eltype=String, n=6)
C: data (eltype=Float64, n=4)
week: data (eltype=Int64, n=6)
mu ~ 1 + cdar(week; by=patch, cor=C)
effect(term_sd, cdar(week), :) ~ Normal(0.0, 0.2)
effect(term_ar, cdar(week), :) ~ Normal(0.8, 0.1)
y ~ Normal(mu, 0.3)SBBRMI with data keys = [:C, :cdar_mu_week_L, :cdar_mu_week_group_idx, :cdar_mu_week_n_groups, :cdar_mu_week_n_steps, :cdar_mu_week_step_idx, :week, :y]
configured submodels:
_sb_cdar_configured_1 = Base.merge(BayesianRegressionModels._sb_cdar, quote
rho ~ normal(0.8, 0.1; lower = 0.0, upper = 1.0)
end)
emitted @slic body:
begin
X_mu = hcat(rep_vector(1.0, num_elements(week)))
pop_mu ~ popefs(; X = X_mu)
cdar_mu_week ~ _sb_cdar_configured_1(; n_groups = cdar_mu_week_n_groups, n_steps = cdar_mu_week_n_steps, L = cdar_mu_week_L, group_idx = cdar_mu_week_group_idx, step_idx = cdar_mu_week_step_idx)
mu = pop_mu + cdar_mu_week
y ~ normal(mu, 0.3)
endfunctions {
matrix hcat(vector x) {
int n = dims(x)[1];
return to_matrix(x, n, 1);
}
vector correlated_damped_walk(
real sigma,
real rho,
vector eta,
matrix L,
int W,
array[] int group_idx,
array[] int step_idx
) {
int P = dims(L)[1];
int N = dims(group_idx)[1];
if (dims(L)[2] != P) reject("correlated_damped_walk: dim mismatch — `L` dim 2 (= ", dims(L)[2], ") does not match `P` (= ", P, "), inferred from `L` dim 1. `P` sizes: `L` dim 1 (= ", dims(L)[1], "), `L` dim 2 (= ", dims(L)[2], ").");
if (dims(step_idx)[1] != N) reject("correlated_damped_walk: dim mismatch — `step_idx` dim 1 (= ", dims(step_idx)[1], ") does not match `N` (= ", N, "), inferred from `group_idx` dim 1. `N` sizes: `group_idx` dim 1 (= ", dims(group_idx)[1], "), `step_idx` dim 1 (= ", dims(step_idx)[1], ").");
matrix[P, W] E = to_matrix(eta, P, W);
matrix[P, W] delta;
vector[dims(L)[1]] prev = (sigma * (L * col(E, 1)));
delta[:, 1] = prev;
real scale = (sigma * sqrt((1.0 - (rho * rho))));
for(w in 2:W) {
prev = ((rho * prev) + (scale * (L * col(E, w))));
delta[:, w] = prev;
}
vector[N] out;
for(i in 1:N) {
out[i] = delta[group_idx[i], step_idx[i]];
}
return out;
}
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 week_n;
array[week_n] int week;
int cdar_mu_week_n_groups;
int cdar_mu_week_n_steps;
int cdar_mu_week_step_idx_n;
int cdar_mu_week_L_m;
int cdar_mu_week_L_n;
matrix[cdar_mu_week_L_m, cdar_mu_week_L_n] cdar_mu_week_L;
int cdar_mu_week_group_idx_n;
array[cdar_mu_week_group_idx_n] int cdar_mu_week_group_idx;
array[cdar_mu_week_step_idx_n] int cdar_mu_week_step_idx;
int y_n;
vector[y_n] y;
}
transformed data {
matrix[num_elements(week), 1] X_mu = hcat(rep_vector(1.0, num_elements(week)));
int pop_mu_n_covariates = 1;
}
parameters {
vector[pop_mu_n_covariates] pop_mu_beta_pop;
real<lower=0.0> cdar_mu_week_sigma;
real<lower=0.0, upper=1.0> cdar_mu_week_rho;
vector[(cdar_mu_week_n_groups * cdar_mu_week_n_steps)] cdar_mu_week_eta;
}
transformed parameters {
vector[num_elements(week)] pop_mu = (X_mu * pop_mu_beta_pop);
vector[cdar_mu_week_step_idx_n] cdar_mu_week = correlated_damped_walk(
cdar_mu_week_sigma,
cdar_mu_week_rho,
cdar_mu_week_eta,
cdar_mu_week_L,
cdar_mu_week_n_steps,
cdar_mu_week_group_idx,
cdar_mu_week_step_idx
);
vector[num_elements(week)] mu = (pop_mu + cdar_mu_week);
}
model {
pop_mu_beta_pop ~ std_normal();
cdar_mu_week_sigma ~ normal(0.0, 0.2);
cdar_mu_week_rho ~ normal(0.8, 0.1);
cdar_mu_week_eta ~ std_normal();
y ~ normal(mu, 0.3);
}
generated quantities {
vector[y_n] y_likelihood = normal_lpdfs(y, mu, 0.3);
vector[y_n] y_gen = normal_vector_rng(y_n, mu, 0.3);
}#= line 0 =# Turing.@model(function brm_model(y, X_mu, terms_mu_1, week, patch, C)
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((:ar, :sd)))}((Distributions.Normal(0.8, 0.1), Distributions.Normal(0.0, 0.2))), NamedTuple{$(QuoteNode((:week, :patch, :C)))}((week, patch, C))))
eta_mu = eta_mu + term_mu_1.effect
mu = eta_mu
begin
for i = Base.eachindex(y)
y[i] ~ Distributions.Normal(mu[i], 0.3)
end
end
(; mu = mu, response = y)
end)HSGP over a model-derived predictor
A one-dimensional hsgp axis may be a sampled linear predictor or formula assignment, not only a raw dataframe column. This makes a latent concentration available to both a linear effect and a residual nonlinear effect in one joint model:
Latent concentration with linear and HSGP effectslatent_hsgp_model = (@brm begin
log(x) ~ 1 + factor(visit) + (1 | xg | subject)
obs_sd ~ Exponential(0.5)
x_obs ~ censored(LogNormal(log(x), obs_sd); lower=lod)
mu ~ 1 + factor(visit) + below_lod + x +
hsgp(x; k=6, domain=(0.02, 4.0), orthogonal_to=:linear) +
(1 + x | grp | subject)
sigma ~ Exponential(1)
y ~ Normal(mu, sigma)
end)((;
visit=repeat([1, 2]; outer=4),
subject=repeat(1:4; inner=2),
below_lod=[1., 1., 1., 1., 0., 0., 0., 0.],
x_obs=[0.25, 0.25, 0.25, 0.25, 0.29, 0.33, 0.41, 0.5],
lod=fill(0.25, 8),
y=[0.9, 1.1, 1.2, 1.4, 1.5, 1.7, 1.9, 2.0],
))BRMI:
visit: data (eltype=Int64, n=8)
subject: data (eltype=Int64, n=8)
log(x) ~ 1 + factor(visit) + (1 | xg | subject)
obs_sd ~ Exponential(0.5)
lod: data (eltype=Float64, n=8)
x_obs ~ censored(LogNormal(log(x), obs_sd); lower=lod)
below_lod: data (eltype=Float64, n=8)
mu ~ 1 + factor(visit) + below_lod + x + hsgp(x; k=6, domain=(0.02, 4.0), orthogonal_to=linear) + ((1 + x) | grp | subject)
sigma ~ Exponential(1)
y ~ Normal(mu, sigma)SBBRMI with data keys = [:below_lod, :col_idx_mu_grp_subject, :lod, :n_subject, :n_terms_grp_subject, :n_terms_xg_subject, :omega2_hsgp_x, :subject_idx, :visit, :visit_idx, :visit_n_levels, :x_obs, :y]
emitted @slic body:
begin
b_xg_subject ~ ranef_correlated_draws(; group_idx = subject_idx, n_groups = n_subject, n_terms = n_terms_xg_subject)
b_grp_subject ~ ranef_correlated_draws(; group_idx = subject_idx, n_groups = n_subject, n_terms = n_terms_grp_subject)
X_log_x = hcat(rep_vector(1.0, num_elements(visit)))
pop_log_x ~ popefs(; X = X_log_x)
cat_log_x_visit ~ _sb_cat(; x = visit_idx, n_levels = visit_n_levels)
r_log_x_xg_subject = b_xg_subject[subject_idx, 1]
log_x = pop_log_x + cat_log_x_visit + r_log_x_xg_subject
x = exp(log_x)
obs_sd ~ exponential(1.0 ./ 0.5)
x_obs ~ censored(lognormal, (log)(x), obs_sd; lower = lod)
X_mu = hcat(rep_vector(1.0, num_elements(below_lod)), below_lod, x)
pop_mu ~ popefs(; X = X_mu)
cat_mu_visit ~ _sb_cat(; x = visit_idx, n_levels = visit_n_levels)
hsgp_x ~ _sb_hsgp_latent_orthogonal(; x = x, omega2 = omega2_hsgp_x, center = 2.01, L = 1.99, rho_lower = 0.9190772477627855)
Z_mu_grp_subject = hcat(rep_vector(1.0, num_elements(subject_idx)), x)
r_mu_grp_subject = rows_dot_product(Z_mu_grp_subject, b_grp_subject[subject_idx, col_idx_mu_grp_subject])
mu = pop_mu + cat_mu_visit + hsgp_x + r_mu_grp_subject
sigma ~ exponential(1.0 ./ 1)
y ~ normal(mu, sigma)
endfunctions {
matrix hcat(vector x) {
int n = dims(x)[1];
return to_matrix(x, n, 1);
}
real lower_clamping_lognormal_lpdf(
vector y,
vector lo,
vector args1,
real args2
) {
return sum(lower_clamping_lognormal_lpdfs(y, lo, args1, args2));
}
vector lower_clamping_lognormal_lpdfs(
vector y,
vector lo,
vector args1,
real args2
) {
int n = dims(y)[1];
return jbroadcasted_lower_clamping_lpdf_lognormal(y, lo, args1, args2);
}
vector jbroadcasted_lower_clamping_lpdf_lognormal(
vector x1,
vector x3,
vector x4,
real x5
) {
int n = dims(x1)[1];
vector[n] rv;
for(i in 1:n) {
rv[i] = lower_clamping_lognormal_lpdf(broadcasted_getindex(x1, i) |
broadcasted_getindex(x3, i),
broadcasted_getindex(x4, i),
x5
);
}
return rv;
}
real lower_clamping_lognormal_lpdf(
real y,
real lo,
real args1,
real args2
) {
array[1] real rv;
rv[1] = negative_infinity();
if((y == lo)) {
rv[1] = lognormal_lcdf(lo | args1, args2);
} else {
if((y > lo)) {
rv[1] = lognormal_lpdf(y | args1, args2);
}
}
return rv[1];
}
real broadcasted_getindex(vector x, int i) {
return x[i];
}
vector lower_clamping_vector_lognormal_rng(
int anontok__1,
vector lo,
vector args1,
real args2
) {
int n = anontok__1;
return jbroadcasted_lower_clamping_cell_rng_lognormal_rng(rep_vector(0.0, n), lo, args1, args2);
}
vector jbroadcasted_lower_clamping_cell_rng_lognormal_rng(
vector x1,
vector x3,
vector x4,
real x5
) {
int n = dims(x1)[1];
vector[n] rv;
for(i in 1:n) {
rv[i] = lower_clamping_cell_lognormal_rng(
broadcasted_getindex(x1, i),
broadcasted_getindex(x3, i),
broadcasted_getindex(x4, i),
x5
);
}
return rv;
}
real lower_clamping_cell_lognormal_rng(
real dummy,
real lo,
real args1,
real args2
) {
return lower_clamping_lognormal_rng(lo, args1, args2);
}
real lower_clamping_lognormal_rng(
real lo,
real args1,
real args2
) {
vector[1] draw;
draw[1] = lognormal_rng(args1, args2);
if((draw[1] < lo)) {
draw[1] = lo;
}
return draw[1];
}
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);
}
matrix brm_hsgp_basis_1d(
vector x,
matrix omega2,
real center,
real L
) {
int n = dims(x)[1];
int m = dims(omega2)[1];
matrix[n, m] PHI;
real inv_sqrt_L = (1.0 / sqrt(L));
for(b in 1:m) {
real omega = sqrt(omega2[b, 1]);
for(i in 1:n) {
PHI[i, b] = (inv_sqrt_L * sin((omega * ((x[i] - center) + L))));
}
}
return PHI;
}
matrix brm_hsgp_orthogonalize_linear(
matrix PHI,
vector x
) {
int n = dims(PHI)[1];
int m = dims(PHI)[2];
if (dims(x)[1] != n) reject("brm_hsgp_orthogonalize_linear: dim mismatch — `x` dim 1 (= ", dims(x)[1], ") does not match `n` (= ", n, "), inferred from `PHI` dim 1. `n` sizes: `PHI` dim 1 (= ", dims(PHI)[1], "), `x` dim 1 (= ", dims(x)[1], ").");
matrix[n, m] out;
vector[dims(x)[1]] x_centered = (x - mean(x));
real x_ss = dot_self(x_centered);
for(b in 1:m) {
vector[dims(PHI)[1]] phi_centered = (col(PHI, b) - mean(col(PHI, b)));
if((x_ss > 1.0e-12)) {
phi_centered = (phi_centered - (x_centered * (dot_product(x_centered, phi_centered) / x_ss)));
}
for(i in 1:n) {
out[i, b] = phi_centered[i];
}
}
return out;
}
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;
}
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);
}
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 n_terms_xg_subject;
int n_subject;
int n_terms_grp_subject;
int visit_n;
array[visit_n] int visit;
int visit_n_levels;
int visit_idx_n;
array[visit_idx_n] int visit_idx;
int subject_idx_n;
array[subject_idx_n] int subject_idx;
int x_obs_n;
vector[x_obs_n] x_obs;
int lod_n;
vector[lod_n] lod;
int below_lod_n;
vector[below_lod_n] below_lod;
int omega2_hsgp_x_m;
int omega2_hsgp_x_n;
matrix[omega2_hsgp_x_m, omega2_hsgp_x_n] omega2_hsgp_x;
int col_idx_mu_grp_subject_n;
array[col_idx_mu_grp_subject_n] int col_idx_mu_grp_subject;
int y_n;
vector[y_n] y;
}
transformed data {
matrix[num_elements(visit), 1] X_log_x = hcat(rep_vector(1.0, num_elements(visit)));
int pop_log_x_n_covariates = 1;
int pop_mu_n_covariates = (2 + 1);
int hsgp_x_n_basis = omega2_hsgp_x_m;
}
parameters {
cholesky_factor_corr[n_terms_xg_subject] b_xg_subject_L;
vector<lower=0.0>[n_terms_xg_subject] b_xg_subject_tau;
vector[(n_terms_xg_subject * n_subject)] b_xg_subject_z_flat;
cholesky_factor_corr[n_terms_grp_subject] b_grp_subject_L;
vector<lower=0.0>[n_terms_grp_subject] b_grp_subject_tau;
vector[(n_terms_grp_subject * n_subject)] b_grp_subject_z_flat;
vector[pop_log_x_n_covariates] pop_log_x_beta_pop;
vector[(visit_n_levels - 1)] cat_log_x_visit_beta;
real<lower=0.0> obs_sd;
vector[pop_mu_n_covariates] pop_mu_beta_pop;
vector[(visit_n_levels - 1)] cat_mu_visit_beta;
real<lower=0.9190772477627855> hsgp_x_rho_iso;
real<lower=0.0> hsgp_x_sigma;
vector[hsgp_x_n_basis] hsgp_x_beta_raw;
real<lower=0.0> sigma;
}
transformed parameters {
matrix[n_terms_xg_subject, n_subject] b_xg_subject_z = to_matrix(b_xg_subject_z_flat, n_terms_xg_subject, n_subject);
matrix[n_subject, n_terms_xg_subject] b_xg_subject = ((diag_pre_multiply(b_xg_subject_tau, b_xg_subject_L) * b_xg_subject_z)');
matrix[n_terms_grp_subject, n_subject] b_grp_subject_z = to_matrix(b_grp_subject_z_flat, n_terms_grp_subject, n_subject);
matrix[n_subject, n_terms_grp_subject] b_grp_subject = ((diag_pre_multiply(b_grp_subject_tau, b_grp_subject_L) * b_grp_subject_z)');
vector[num_elements(visit)] pop_log_x = (X_log_x * pop_log_x_beta_pop);
vector[visit_idx_n] cat_log_x_visit = append_row(0.0, cat_log_x_visit_beta)[visit_idx];
vector[subject_idx_n] r_log_x_xg_subject = b_xg_subject[subject_idx, 1];
vector[num_elements(visit)] log_x = (pop_log_x + cat_log_x_visit + r_log_x_xg_subject);
vector[num_elements(visit)] x = exp(log_x);
matrix[num_elements(visit), (2 + 1)] X_mu = hcat(rep_vector(1.0, num_elements(below_lod)), below_lod, x);
vector[num_elements(visit)] pop_mu = (X_mu * pop_mu_beta_pop);
vector[visit_idx_n] cat_mu_visit = append_row(0.0, cat_mu_visit_beta)[visit_idx];
matrix[num_elements(visit), omega2_hsgp_x_m] hsgp_x_PHI_raw = brm_hsgp_basis_1d(x, omega2_hsgp_x, 2.01, 1.99);
matrix[num_elements(visit), omega2_hsgp_x_m] hsgp_x_PHI = brm_hsgp_orthogonalize_linear(hsgp_x_PHI_raw, x);
vector[1] hsgp_x_rho = rep_vector(hsgp_x_rho_iso, 1);
vector[omega2_hsgp_x_m] hsgp_x_sqrt_spd = brm_hsgp_sqrt_spd(omega2_hsgp_x, hsgp_x_sigma, hsgp_x_rho);
vector[num_elements(visit)] hsgp_x = (hsgp_x_PHI * (hsgp_x_sqrt_spd .* hsgp_x_beta_raw));
matrix[num_elements(visit), 2] Z_mu_grp_subject = hcat(rep_vector(1.0, num_elements(subject_idx)), x);
vector[subject_idx_n] r_mu_grp_subject = rows_dot_product(Z_mu_grp_subject, b_grp_subject[subject_idx, col_idx_mu_grp_subject]);
vector[num_elements(visit)] mu = (pop_mu + cat_mu_visit + hsgp_x + r_mu_grp_subject);
}
model {
b_xg_subject_L ~ lkj_corr_cholesky(1.0);
b_xg_subject_tau ~ std_normal();
b_xg_subject_z_flat ~ std_normal();
b_grp_subject_L ~ lkj_corr_cholesky(1.0);
b_grp_subject_tau ~ std_normal();
b_grp_subject_z_flat ~ std_normal();
pop_log_x_beta_pop ~ std_normal();
cat_log_x_visit_beta ~ std_normal();
obs_sd ~ exponential((1.0 ./ 0.5));
x_obs ~ lower_clamping_lognormal(lod, log(x), obs_sd);
pop_mu_beta_pop ~ std_normal();
cat_mu_visit_beta ~ std_normal();
hsgp_x_rho_iso ~ lognormal(0.0, 1.0);
hsgp_x_sigma ~ lognormal(0.0, 1.0);
hsgp_x_beta_raw ~ std_normal();
sigma ~ exponential((1.0 ./ 1));
y ~ normal(mu, sigma);
}
generated quantities {
vector[x_obs_n] x_obs_likelihood = lower_clamping_lognormal_lpdfs(x_obs, lod, log(x), obs_sd);
vector[x_obs_n] x_obs_gen = lower_clamping_vector_lognormal_rng(x_obs_n, lod, log(x), obs_sd);
vector[y_n] y_likelihood = normal_lpdfs(y, mu, sigma);
vector[y_n] y_gen = normal_vector_rng(y_n, mu, sigma);
}Turing unsupported for this BRM example
BRM backend lowering: random-slope designs support the backend-neutral population column surface (continuous, categorical, transformed, and interaction columns); got `x` in `((1 + x) | grp | subject)`The explicit domain=(lower, upper) is required because sampled x values do not exist when Julia configures the basis. It is the actual compact HSGP interval, so it cannot be combined with the raw-data expansion factor c and remains fixed during reprocess, including with freeze_constants=false. Choose it to cover scientifically plausible posterior support for every latent x: the domain configures the approximation but does not truncate or otherwise constrain x at runtime.
orthogonal_to=:linear centers every basis column and projects it off the current sampled x direction at each draw. Use it when the formula also contains x; the population coefficient then carries the linear association and the HSGP carries only residual nonlinear shape. The option also works for a one-dimensional raw axis. The orthogonality option is ungrouped: multiplying the basis by group-specific weights would no longer preserve the global projection. Model-derived HSGPs are currently one-dimensional, isotropic, and ungrouped; raw-data HSGPs retain their variadic, anisotropic, and group-specific forms when orthogonal_to is omitted.
After fitting, hsgp_population_curve evaluates the combined population exposure contribution on a fixed grid without fabricating assay rows or referring to emitted Stan names:
curve = hsgp_population_curve(
descriptor, constrained_draws, constrained_names,
collect(range(0.02, 4.0; length=100));
predictor=:mu, coefficient=:x, term=:hsgp_x)
curve.linear # beta_x * x, draws × grid
curve.hsgp # orthogonal residual nonlinear contribution, draws × grid
curve.total # the supported total exposure curveconstrained_draws and constrained_names must include transformed parameters (include_tp=true in BridgeStan), because every posterior draw has its own sampled training x. BRM uses those values to replay the exact fitted intercept/linear projection on the new grid; re-orthogonalizing against the grid would define a different curve. Evaluation outside the formula's fixed domain is rejected. This is a population partial effect only: it deliberately excludes the intercept, other covariates, and subject-specific random slopes.
Fixed partial centering of HSGP weights
For a raw, ungrouped squared-exponential HSGP, centeredness chooses the coordinate of each basis weight without changing its physical prior. It may be one scalar shared by all weights or a vector of length prod(k), supplied literally or through a data column. With spectral standard deviation s, unit normal z, and centeredness c in [0,1], BRM samples and reconstructs
u ~ Normal(0, s^c)
w = s^(1-c) u = s zThus c=0 is the default noncentered coordinate and c=1 is centered; intermediate values partially center that frequency. Both StanBlocks and Turing use the same log-spectral-scale calculation and endpoint-safe c=0 path, so a numerically vanished high-frequency scale cannot turn 0 * -Inf into NaN.
select_hsgp_centeredness applies the pilot rule used in the adaptive HSGP case study: rows are pilot draws, columns are basis frequencies, and candidates are fixed before the refit. The selector is not an online warmup controller. Periodic, latent-input, group-specific, and orthogonal_to HSGPs reject nonzero partial centering until those geometries have their own verified coordinate contract.
Interval-censored predictor
Use interval_censored(x; upper=lloq) when x is quantified on some rows and BLOQ on others. Store the measured concentration on quantified rows and the row-specific LLOQ on BLOQ rows; the convention x == lloq identifies BLOQ, so no separate flag is needed. Each BLOQ row allocates one latent covariate value between zero and its LLOQ, and the merged vector enters the ordinary fixed- or random-effect design matrix:
Interval-censored concentration predictorinterval_censored_predictor_model = (@brm begin
y ~ Normal(mu, sigma)
mu ~ 1 + interval_censored(conc; upper=lloq)
effect(mu, conc) ~ Normal(0, 2)
latent(mu, interval_censored(conc)) ~ Normal(0, 5)
sigma ~ Exponential(1)
end)((;
y=[2.1, 2.4, 2.9],
conc=[0.7, 0.25, 1.2],
lloq=[0.25, 0.25, 0.4],
))BRMI:
mu ~ 1 + interval_censored(conc; upper=lloq)
sigma ~ Exponential(1)
y ~ Normal(mu, sigma)
lloq: data (eltype=Float64, n=3)
conc: data (eltype=Float64, n=3)
effect(mu, conc) ~ Normal(0, 2)
effect(term_latent, interval_censored(conc), mu) ~ Normal(0, 5)SBBRMI with data keys = [:Jexact_conc_lloq, :Jinterval_conc_lloq, :conc, :icp_exact_conc_lloq, :icp_lower_conc_lloq, :icp_upper_conc_lloq, :lloq, :y]
emitted @slic body:
begin
sigma ~ exponential(1.0 ./ 1)
interval_censored_conc_lloq ~ _sb_interval_censored_predictor(; x_exact = icp_exact_conc_lloq, x_lower = icp_lower_conc_lloq, x_upper = icp_upper_conc_lloq, Jexact = Jexact_conc_lloq, Jinterval = Jinterval_conc_lloq, x_true_loc = 0, x_true_scale = 5)
X_mu = hcat(rep_vector(1.0, num_elements(conc)), interval_censored_conc_lloq)
pop_mu ~ _popefs_normal_coefs(; X = X_mu, beta_loc = [0.0, 0], beta_scale = [1.0, 2])
y ~ normal_id_glm(X_mu, 0.0, pop_mu, sigma)
mu = X_mu * pop_mu
endfunctions {
real conditioning_normal_lpdf(
vector y,
vector lo,
vector hi,
int args1,
int args2
) {
return sum(conditioning_normal_lpdfs(y, lo, hi, args1, args2));
}
vector conditioning_normal_lpdfs(
vector y,
vector lo,
vector hi,
int args1,
int args2
) {
int n = dims(y)[1];
return jbroadcasted_conditioning_lpdf_normal(y, lo, hi, args1, args2);
}
vector jbroadcasted_conditioning_lpdf_normal(
vector x1,
vector x3,
vector x4,
int x5,
int x6
) {
int n = dims(x1)[1];
vector[n] rv;
for(i in 1:n) {
rv[i] = conditioning_normal_lpdf(broadcasted_getindex(x1, i) |
broadcasted_getindex(x3, i),
broadcasted_getindex(x4, i),
x5,
x6
);
}
return rv;
}
real conditioning_normal_lpdf(
real y,
real lo,
real hi,
int args1,
int args2
) {
array[1] real rv;
rv[1] = negative_infinity();
if((lo >= hi)) {
reject("truncated: lower bound must be less than upper bound");
} else {
if((y >= lo)) {
if((y <= hi)) {
rv[1] = (
normal_lpdf(y | args1, args2) -
log_diff_exp(normal_lcdf_stable(hi, args1, args2), normal_lcdf_stable(lo, args1, args2))
);
}
}
}
return rv[1];
}
real normal_lcdf_stable(
real x,
int loc,
int scale
) {
return (log(erfc(((-(x - loc)) / (scale * sqrt(2.0))))) - log(2.0));
}
real broadcasted_getindex(vector x, int i) {
return x[i];
}
vector mi_merge(
vector y_obs,
vector y_mis,
array[] int Jobs,
array[] int Jmis,
int n
) {
int n_obs = dims(y_obs)[1];
int n_mis = dims(y_mis)[1];
if (dims(Jobs)[1] != n_obs) reject("mi_merge: dim mismatch — `Jobs` dim 1 (= ", dims(Jobs)[1], ") does not match `n_obs` (= ", n_obs, "), inferred from `y_obs` dim 1. `n_obs` sizes: `y_obs` dim 1 (= ", dims(y_obs)[1], "), `Jobs` dim 1 (= ", dims(Jobs)[1], ").");
if (dims(Jmis)[1] != n_mis) reject("mi_merge: dim mismatch — `Jmis` dim 1 (= ", dims(Jmis)[1], ") does not match `n_mis` (= ", n_mis, "), inferred from `y_mis` dim 1. `n_mis` sizes: `y_mis` dim 1 (= ", dims(y_mis)[1], "), `Jmis` dim 1 (= ", dims(Jmis)[1], ").");
vector[n] rv = rep_vector(0.0, n);
if((n_obs > 0)) {
for(i in 1:n_obs) {
rv[Jobs[i]] = y_obs[i];
}
}
if((n_mis > 0)) {
for(i in 1:n_mis) {
rv[Jmis[i]] = y_mis[i];
}
}
return rv;
}
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,
real 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,
real 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,
real 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 Jinterval_conc_lloq_n;
array[Jinterval_conc_lloq_n] int Jinterval_conc_lloq;
int icp_lower_conc_lloq_n;
vector[icp_lower_conc_lloq_n] icp_lower_conc_lloq;
int icp_upper_conc_lloq_n;
vector[icp_upper_conc_lloq_n] icp_upper_conc_lloq;
int Jexact_conc_lloq_n;
array[Jexact_conc_lloq_n] int Jexact_conc_lloq;
int icp_exact_conc_lloq_n;
vector[icp_exact_conc_lloq_n] icp_exact_conc_lloq;
int conc_n;
vector[conc_n] conc;
int y_n;
vector[y_n] y;
}
transformed data {
int interval_censored_conc_lloq_n_interval = num_elements(Jinterval_conc_lloq);
int pop_mu_n_covariates = 2;
}
parameters {
real<lower=0.0> sigma;
vector<lower=icp_lower_conc_lloq, upper=icp_upper_conc_lloq>[interval_censored_conc_lloq_n_interval] interval_censored_conc_lloq_x_interval;
vector[pop_mu_n_covariates] pop_mu_beta_pop;
}
transformed parameters {
vector[(num_elements(Jexact_conc_lloq) + interval_censored_conc_lloq_n_interval)] interval_censored_conc_lloq = mi_merge(
icp_exact_conc_lloq,
interval_censored_conc_lloq_x_interval,
Jexact_conc_lloq,
Jinterval_conc_lloq,
(num_elements(Jexact_conc_lloq) + interval_censored_conc_lloq_n_interval)
);
matrix[(num_elements(Jexact_conc_lloq) + interval_censored_conc_lloq_n_interval), 2] X_mu = hcat(rep_vector(1.0, num_elements(conc)), interval_censored_conc_lloq);
vector[pop_mu_n_covariates] pop_mu = pop_mu_beta_pop;
}
model {
sigma ~ exponential((1.0 ./ 1));
interval_censored_conc_lloq_x_interval ~ conditioning_normal(icp_lower_conc_lloq, icp_upper_conc_lloq, 0, 5);
pop_mu_beta_pop ~ normal([0.0, 0]', [1.0, 2]');
y ~ normal_id_glm(X_mu, 0.0, pop_mu, sigma);
}
generated quantities {
vector[(num_elements(Jexact_conc_lloq) + interval_censored_conc_lloq_n_interval)] y_likelihood = normal_id_glm_lpdfs(y, X_mu, 0.0, pop_mu, sigma);
vector[(num_elements(Jexact_conc_lloq) + interval_censored_conc_lloq_n_interval)] y_gen = normal_id_glm_vector_rng(y_n, X_mu, 0.0, pop_mu, sigma);
vector[(num_elements(Jexact_conc_lloq) + interval_censored_conc_lloq_n_interval)] mu = (X_mu * pop_mu);
}Turing unsupported for this BRM example
BRM preparation: `conc` is not a population coefficient of `mu`. Available labels: Intercept.The term's default latent prior is Normal(0, 1); use latent(<lp|:>, interval_censored(x)) ~ Normal(location, scale) to set it. The population slope keeps the ordinary address effect(<lp>, x). Reusing the same term as a fixed effect and a random slope shares one latent covariate vector. reprocess rebuilds the quantified/BLOQ split from new x and lloq columns.
This is a continuous latent-covariate model, not an LLOQ/2 substitution. For a BLOQ assay row the lower bound defaults to zero and can be changed with a numeric lower=. Values and LLOQs must be finite; x must equal its LLOQ on BLOQ rows and exceed it on quantified rows, and at least one row must be BLOQ. The predictor form is SBBRMI-only and is distinct from the response-likelihood form documented under Likelihoods.
gp and hsgp are distinct terms with no compatibility alias. Both are direct predictor summands carrying their own latent draws and hyperparameters, so neither contributes a beta_pop coefficient. jitter belongs only to gp; k, c and by only to hsgp. Both support cov=:exp_quad (the default) and cov=:periodic; see Periodic covariance for the latter's contract.
Periodic covariance
cov=:periodic selects Stan's periodic kernel
k(x, x') = σ² exp(−2 sin²(π |x − x'| / period) / ρ²)over exactly one axis. period is a required numeric formula constant on the axis's own scale (hours, days, radians — whatever the column carries), and rho / sigma keep their meaning and their length_scale(…) / sd(…) addresses. It is an explicit assumption that the effect repeats with that period; clock-time or dosing-phase evidence is what makes a 24-hour period on elapsed time readable as shared diurnal variation.
Periodic (24 h) time effect over hours since first doseperiodic_time_model = (@brm begin
mu ~ 1 + hsgp(hours_since_dose; k=8, cov=:periodic, period=24.0)
length_scale(:, hsgp(hours_since_dose)) ~ LogNormal(0, 0.5)
sd(:, hsgp(hours_since_dose)) ~ Normal(0, 0.3)
sigma ~ Exponential(1)
y ~ Normal(mu, sigma)
end)((;
hours_since_dose=[0.5, 1.0, 2.0, 4.0, 8.0, 12.0, 24.5, 48.0, 72.0, 168.0],
y=[0.2, 0.9, 1.1, 0.8, 0.4, -0.6, 0.3, 0.1, -0.2, 0.0],
))BRMI:
hours_since_dose: data (eltype=Float64, n=10)
mu ~ 1 + hsgp(hours_since_dose; k=8, cov=periodic, period=24.0)
effect(term_length_scale, hsgp(hours_since_dose), :) ~ LogNormal(0, 0.5)
effect(term_sd, hsgp(hours_since_dose), :) ~ Normal(0, 0.3)
sigma ~ Exponential(1)
y ~ Normal(mu, sigma)SBBRMI with data keys = [:PHI_hsgp_hours_since_dose, :harmonics_hsgp_hours_since_dose, :hours_since_dose, :rho_lower_hsgp_hours_since_dose, :y]
configured submodels:
_sb_hsgp_periodic_configured_1 = Base.merge(BayesianRegressionModels._sb_hsgp_periodic, quote
rho_iso ~ lognormal(0.0, 0.5; lower = 0.0)
sigma ~ normal(0.0, 0.3; lower = 0.0)
end)
emitted @slic body:
begin
X_mu = hcat(rep_vector(1.0, num_elements(hours_since_dose)))
pop_mu ~ popefs(; X = X_mu)
hsgp_hours_since_dose ~ _sb_hsgp_periodic_configured_1(; PHI = PHI_hsgp_hours_since_dose, harmonics = harmonics_hsgp_hours_since_dose, rho_lower = rho_lower_hsgp_hours_since_dose)
mu = pop_mu + hsgp_hours_since_dose
sigma ~ exponential(1.0 ./ 1)
y ~ normal(mu, sigma)
endfunctions {
matrix hcat(vector x) {
int n = dims(x)[1];
return to_matrix(x, n, 1);
}
vector brm_hsgp_periodic_sqrt_spd(
vector harmonics,
real sigma,
real rho
) {
int m = dims(harmonics)[1];
vector[m] rv;
real a = (1.0 / (rho * rho));
real base = (log(sigma) + (0.5 * (log(2.0) - a)));
for(b in 1:m) {
rv[b] = exp((base + (0.5 * log_modified_bessel_first_kind(harmonics[b], a))));
}
return rv;
}
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 hours_since_dose_n;
vector[hours_since_dose_n] hours_since_dose;
int harmonics_hsgp_hours_since_dose_n;
vector[harmonics_hsgp_hours_since_dose_n] harmonics_hsgp_hours_since_dose;
int PHI_hsgp_hours_since_dose_m;
int PHI_hsgp_hours_since_dose_n;
matrix[PHI_hsgp_hours_since_dose_m, PHI_hsgp_hours_since_dose_n] PHI_hsgp_hours_since_dose;
int y_n;
vector[y_n] y;
}
transformed data {
matrix[num_elements(hours_since_dose), 1] X_mu = hcat(rep_vector(1.0, num_elements(hours_since_dose)));
int pop_mu_n_covariates = 1;
int hsgp_hours_since_dose_n_basis = harmonics_hsgp_hours_since_dose_n;
}
parameters {
vector[pop_mu_n_covariates] pop_mu_beta_pop;
real<lower=0.0> hsgp_hours_since_dose_rho_iso;
real<lower=0.0> hsgp_hours_since_dose_sigma;
vector[hsgp_hours_since_dose_n_basis] hsgp_hours_since_dose_beta_raw;
real<lower=0.0> sigma;
}
transformed parameters {
vector[num_elements(hours_since_dose)] pop_mu = (X_mu * pop_mu_beta_pop);
vector[harmonics_hsgp_hours_since_dose_n] hsgp_hours_since_dose_sqrt_spd = brm_hsgp_periodic_sqrt_spd(
harmonics_hsgp_hours_since_dose,
hsgp_hours_since_dose_sigma,
hsgp_hours_since_dose_rho_iso
);
vector[PHI_hsgp_hours_since_dose_m] hsgp_hours_since_dose = (PHI_hsgp_hours_since_dose * (hsgp_hours_since_dose_sqrt_spd .* hsgp_hours_since_dose_beta_raw));
vector[num_elements(hours_since_dose)] mu = (pop_mu + hsgp_hours_since_dose);
}
model {
pop_mu_beta_pop ~ std_normal();
hsgp_hours_since_dose_rho_iso ~ lognormal(0.0, 0.5);
hsgp_hours_since_dose_sigma ~ normal(0.0, 0.3);
hsgp_hours_since_dose_beta_raw ~ std_normal();
sigma ~ exponential((1.0 ./ 1));
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);
}#= line 0 =# Turing.@model(function brm_model(y, X_mu, terms_mu_1, hours_since_dose)
beta_pop ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 1))
eta_mu = X_mu * beta_pop
term_mu_1 ~ DynamicPPL.to_submodel(BayesianRegressionModelsTuringExt._brm_turing_hsgp_ncp_model(terms_mu_1, Base.length(y), NamedTuple{$(QuoteNode((:rho, :sigma)))}((Distributions.LogNormal(0, 0.5), Distributions.Normal(0, 0.3))), NamedTuple{$(QuoteNode((:hours_since_dose,)))}((hours_since_dose,))))
eta_mu = eta_mu + term_mu_1.effect
mu = eta_mu
sigma ~ Distributions.Exponential(1)
begin
for i = Base.eachindex(y)
y[i] ~ Distributions.Normal(mu[i], sigma)
end
end
(; mu = mu, sigma = sigma, response = y)
end)gp(x; cov=:periodic, period=…)is the exact kernel, lowered to Stan's nativegp_periodic_covplusjitter.hsgp(x; k, cov=:periodic, period=…)is the Hilbert-space approximation of Riutort-Mayol et al. (2023): withw0 = 2π / periodthe basis iscos(j·w0·x)andsin(j·w0·x)forj = 1:k—2kfunctions, emitted as dataPHI_hsgp_<x>— and each harmonic's spectral weight isσ·√(2·e^{−a}·I_j(a))witha = 1/ρ², computed inside Stan throughlog_modified_bessel_first_kindso a small length scale cannot overflow.The term contains no constant. The
j = 0harmonic is excluded, so every basis function has zero mean over a period and the term behaves like a0 +term: a nonzero mean of the periodic function belongs to the formula intercept (or whichever declared quantity owns it), never to this term, and nothing is added to the consumer's intercept silently. Count2kfunctions against a basis budget (16 fork = 8).Truncation is the only approximation, and it is measured. The
k-harmonic expansion is exactly periodic but approximates the exp-sine-squared kernel with a tailσ²·2·e^{−a}·Σ_{j>k} I_j(a)that shrinks asrhogrows. Dropped fraction of the (non-constant) variance:
k (harmonics) | basis functions | default floor on rho | Riutort-Mayol B.6 (3.72/k) | dropped variance fraction at rho = 0.3 / 0.47 / 0.6 / 0.8 / 1.0 / 1.5 |
|---|---|---|---|---|
| 4 | 8 | 1.931 | 0.930 | 2.0e-1 / 4.4e-2 / 1.3e-2 / 2.0e-3 / 4.1e-4 / 1.9e-5 |
| 8 | 16 | 0.614 | 0.465 | 1.3e-2 / 2.4e-4 / 1.2e-5 / 2.2e-7 / 8.0e-9 / 1.5e-11 |
| 16 | 32 | 0.278 | 0.233 | 3.1e-6 / 1.2e-10 / 1.5e-13 / 3.0e-17 / 3.0e-20 / 8.7e-26 |
| 32 | 64 | 0.135 | 0.116 | 4.3e-17 / 1.8e-27 / 1.0e-33 / 2.2e-41 / 1.9e-47 / 1.3e-58 |
The default floor (next section) is the rho at which the k-th harmonic's amplitude is 1/100 of the first's; it is stricter than the paper's k ≥ 3.72/rho rule. At k = 8 the maximum covariance error over a full period is 9e-6·σ² at the floor and 1e-2·σ² at rho = 0.3, where an explicit length_scale override would take you knowingly.
rhois dimensionless: it is the kernel's length scale on the unit circle, in chord units (2·sin(π·τ/period)is the chord spanned by a lagτ), never in the axis's own units —rho ≈ 1means features of roughly one radian, i.e.period/(2π)on the axis. Statelength_scale(…)priors in those units. The raw axis values go straight into the cosines and sines, so a large elapsed time needs no consumer-sidemod(x, period).The periodic basis needs no boundary factor and no domain:
c,domain,orthogonal_to,by, andiso=falseare refused by name, and a model-derived axis is not supported. Cosines and sines already have zero mean over a period and are not collinear with a linearx.reprocessrebuilds the cosine/sine columns from the new axis. Nothing here is fitted, sofreeze_constants=trueandfalseare identical, a constant prediction axis is valid, and rows far outside the training range are the point rather than an extrapolation hazard.brm_term_coordinatesresolves:length_scale,:sd, and:basis_weightsexactly as for the exp-quad basis;:basis_weightshas2kcoordinates, andterm_drawszeroes all of them.The validity floor below applies with the periodic weights substituted:
rho_lower_hsgp_<x>is the length scale at whichI_k(a)/I_1(a) = 100⁻², a function ofkalone.StanBlocks floor. The weight helper calls Stan's
log_modified_bessel_first_kind, registered in StanBlocks frombec23bc3c52303ebde60a026af48c435e4c81330(devibe, 2026-09-07). On an older StanBlocks the periodichsgpform fails at transpile withCould not find log_modified_bessel_first_kind …; the exactgp(…; cov=:periodic)form needs nothing new.
Column transforms
Applied to a raw column before it enters the design matrix. The constants are fitted on the training frame and frozen, so reprocess replays them rather than re-deriving them — see the replay contract below.
| Term | What it does |
|---|---|
zscale(x) | subtract the mean, divide by the SD |
center(x) | subtract the mean |
standardize(x) | mean/SD standardisation |
Grouping factors
Used on the right of | in a random-effect block.
| Term | What it does |
|---|---|
gr(g) | ordinary grouping factor / strata |
mm(g1, g2, …; weights=(w1, w2, …), normalize=true) | multi-membership — one coefficient block shared across two or more levels per row |
mm needs at least two group columns; omitting weights gives exact equal weights 1/M. Supplied weights must be present, real, finite, nonnegative and sum to something positive on every row. It does not combine with |ID|, cv_groups or centered_groups, which error explicitly.
Exact total coefficients
The StanBlocks backend uses total coefficients by default for eligible independent random effects. For example, in
Exact total intercepts with a fixed slopebuilder = @brm begin
mu ~ 1 + x + (1 | subject_effect | subject)
effect(mu, Intercept) ~ Normal(0, 3)
sd(:, subject_effect) ~ Normal(0, 1)
y ~ Normal(mu, 1)
end
data = (; subject=[1, 1, 2, 2, 3, 3], x=[0., 1., 0., 1., 0., 1.],
y=[0.2, 0.8, -0.1, 0.9, 0.4, 1.2])
total_intercept_model = builder(data)BRMI:
x: data (eltype=Float64, n=6)
subject: data (eltype=Int64, n=6)
mu ~ 1 + x + (1 | subject_effect | subject)
effect(mu, Intercept) ~ Normal(0, 3)
effect(sd, subject_effect) ~ Normal(0, 1)
y ~ Normal(mu, 1)SBBRMI with data keys = [:subject, :total_A_mu, :total_group_mu, :total_location_mu, :total_ng_mu, :total_nk_mu, :total_np_mu, :total_precision_mu, :x, :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_mu ~ _brm_total_scales_configured_1(; n = total_nk_mu)
total_mu::matrix[total_ng_mu, total_nk_mu] ~ brm_total(total_scale_mu, total_A_mu, total_location_mu, total_precision_mu)
population_mu = brm_total_recover_rng(total_mu, total_scale_mu, total_A_mu, total_location_mu, total_precision_mu)
deviation_mu = brm_total_deviations(total_mu, total_A_mu * population_mu)
total_Z_mu = hcat(rep_vector(1.0, num_elements(total_group_mu)))
X_mu = hcat(x)
pop_mu ~ popefs(; X = X_mu)
mu = rows_dot_product(total_mu[total_group_mu, :], total_Z_mu) + pop_mu
y ~ normal(mu, 1)
endfunctions {
// 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,
int scale
) {
return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
vector x1,
vector x2,
int 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,
int 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,
int 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_mu;
int total_nk_mu;
int total_A_mu_m;
int total_A_mu_n;
matrix[total_A_mu_m, total_A_mu_n] total_A_mu;
int total_location_mu_n;
vector[total_location_mu_n] total_location_mu;
int total_precision_mu_n;
vector[total_precision_mu_n] total_precision_mu;
int total_group_mu_n;
array[total_group_mu_n] int total_group_mu;
int x_n;
vector[x_n] x;
int y_n;
vector[y_n] y;
}
transformed data {
matrix[num_elements(total_group_mu), 1] total_Z_mu = hcat(rep_vector(1.0, num_elements(total_group_mu)));
matrix[x_n, 1] X_mu = hcat(x);
int pop_mu_n_covariates = 1;
}
parameters {
vector<lower=0.0>[1] total_scale_mu_tau;
matrix[total_ng_mu, total_nk_mu] total_mu;
vector[pop_mu_n_covariates] pop_mu_beta_pop;
}
transformed parameters {
vector<lower=0.0>[1] total_scale_mu = total_scale_mu_tau;
vector[x_n] pop_mu = (X_mu * pop_mu_beta_pop);
vector[num_elements(total_group_mu)] mu = (rows_dot_product(total_mu[total_group_mu, :], total_Z_mu) + pop_mu);
}
model {
total_scale_mu_tau ~ brm_vector_prior_ca4b8a1c1bc116d6(0.0, 1.0);
total_mu ~ brm_total(total_scale_mu, total_A_mu, total_location_mu, total_precision_mu);
pop_mu_beta_pop ~ std_normal();
y ~ normal(mu, 1);
}
generated quantities {
vector[total_precision_mu_n] population_mu = brm_total_recover_rng(total_mu, total_scale_mu, total_A_mu, total_location_mu, total_precision_mu);
matrix[total_ng_mu, total_A_mu_m] deviation_mu = brm_total_deviations(total_mu, (total_A_mu * population_mu));
vector[y_n] y_likelihood = normal_lpdfs(y, mu, 1);
vector[y_n] y_gen = normal_vector_rng(y_n, mu, 1);
}#= line 0 =# Turing.@model(function brm_model(y, X_mu, group_effects_mu_1)
beta_pop ~ Distributions.product_distribution([Distributions.Normal(0, 3), Distributions.Normal(0, 1)])
eta_mu = X_mu * beta_pop
group_effect_1 = Base.zeros(Base.length(y))
group_1_1 ~ DynamicPPL.to_submodel(BRM.turing_group_effect(group_effects_mu_1, (Distributions.Normal(0, 1),), nothing))
group_effect_1 = group_effect_1 + group_1_1.effect
eta_mu = eta_mu + group_effect_1
mu = eta_mu
begin
for i = Base.eachindex(y)
y[i] ~ Distributions.Normal(mu[i], 1)
end
end
(; mu = mu, response = y)
end)the sampled group intercept is total[j] = intercept + deviation[j]. The fixed slope of x remains explicit. BRM integrates out the shared intercept using its original prior; this preserves the likelihood and posterior of the totals and hyperparameters. All groups have the same status, and the sampler uses one total per group without an extra mean direction.
With intercept ~ Normal(m, s) and independent deviation[j] ~ Normal(0, tau), the exact induced prior is
The density uses group sums and a small population-coefficient precision matrix, so its work is linear in the number of groups for a fixed number of coefficients. Student-t population priors (TDist or LocationScale(..., TDist(...))) use an exact Gaussian scale mixture; Flat() denotes an improper uniform population prior. Marginalization preserves the supported population prior specified in the formula.
total_effect_blocks(sb) reports the selected blocks. total_groups=() opts out; total_groups=[:subject] requires that group to be eligible. Explicit centered_groups=[:subject] selects conventional centered deviations.
Automatic selection currently requires one grouping structure per predictor, independent margins, supported population priors, and a known relation between the population and random-effect columns. Raw numeric columns, pure data expressions, and centered/scaled population columns with matching raw random columns are supported. Correlated, crossed, stratified, multi-membership and R2D2 blocks use the conventional representation. A request for totals on an unsupported group raises an error.
Automatic selection retains the conventional representation if some ordinary group blocks require it. A single adaptive wrapper currently combines totals with S2Z contrasts, or ordinary random effects with HSGP weights, but not both pairs.
WarmupHMC and recovery
BRM discovers each total's scale and supplies independent centering controls:
using Random, Enzyme, WarmupHMC, BridgeStan, StanBlocks
using DifferentiationInterface: AutoEnzyme
sb = SBBRMI(total_intercept_model)
problem = StanBlocks.stan_instantiate(sb.model)
backend = AutoEnzyme(;
mode=Enzyme.set_runtime_activity(Enzyme.Reverse),
function_annotation=Enzyme.Const)
adaptive = adaptive_centering_problem(sb, problem, backend;
centeredness=0.0)
fit = WarmupHMC.adaptive_warmup_mcmc(Xoshiro(1), adaptive;
n_draws=2000, nonlinear_adapt=true)centeredness=1 means model-scale totals; 0 means scaled totals. Intermediate values may differ by group and coefficient. The marginal totals remain correlated under their exact prior, including at the NCP endpoint. Online selection uses WarmupHMC's default position-gradient loss. For a fixed fit, use nonlinear_adapt=false.
For a post-hoc refit, call select_total_centeredness(sb, draws, names) on a pilot in the compiled model frame, then pass its centeredness vector to adaptive_centering_problem. The default criterion=:position needs positions; criterion=:gradient additionally needs matching compiled-frame gradients. Both inputs are draws × coordinates. Selection makes no additional gradient calls; include the pilot's cost when comparing full workflows.
names = BridgeStan.param_unc_names(problem.model)
draws = permutedims(fit.posterior_position)
recovered = recover_population_draws(sb, draws, names; rng=Xoshiro(2))
recovered[:mu].population # original shared intercept, conditionally recovered
recovered[:mu].totals # sampled subject-specific intercepts
recovered[:mu].deviations # totals minus the recovered shared interceptThe compiled model also exposes the recovered coefficients and deviations in generated quantities. Recovery adds conditional randomness, so compare sampler efficiency on consistent scientific quantities and distinguish recovered means from sampled totals. The pupil study explains this distinction.
population_draws(...; groups=:subject, rng=...) substitutes a recovered population draw for each group's total. For new groups, use the reusable generative_plan(builder, data) form, rebuild with generative_plan(plan, new_data), and call transport_draws. Existing groups retain their totals; new groups share one recovered population draw per posterior draw. resample=:subject redraws existing groups too. Total blocks require frozen preprocessing for replay and use this transport route for resampling; reprocess(...; resample_groups=...) is not supported for them. Improper population priors have no prior-predictive distribution.
Posterior-preserving sum-to-zero (S2Z) effects
SBBRMI(brmi; s2z_groups=[:g], s2z_rho=...) opts a grouping factor into the S2Z construction of brms PR #1919. Each coefficient's J group effects become J - 1 free orthonormal (Helmert) contrasts plus a block mean. The mean is integrated into the population coefficients exactly, as for totals, and is recovered afterwards. The current scope is one independent Gaussian grouping structure per predictor with J >= 2, a population column for every group column, and Normal or Flat() population priors (Student-t priors are not yet supported).
s2z_rho sets the compiled frame per group and coefficient: a scalar, one weight per coefficient, or a J × K matrix in [0, 1]. Intermediate values use Sean's projected partial map. select_s2z_rho chooses these weights from a pilot with brms's Fisher rule.
Sean's rule at every warm-up window
adaptive_centering_problem(sb, problem, backend; s2z_rule=:fisher) keeps the contrast coordinates and Sean's projected partial map, with one weight per group and coefficient, and re-selects those weights at every restarting WarmupHMC window using Sean's rule from brms PR #1919. At each window boundary, every evidence draw is mapped to the compiled frame. The per-row expected information is evaluated there at the draw's fitted values through BridgeStan and accumulated per group through the design. It is turned into Sean's per-draw weights, rescaled at that draw's tau, and the new weight of each cell is the median across draws.
The evidence is the retained pool by default. With nonlinear_evidence=:nuts_weighted or :all_good_leaves, WarmupHMC hands over an equally sized (recording_target) weight-proportional sample of the window's NUTS leaves. No Pathfinder precursor or importance sampling is involved. The compiled s2z_rho is the target frame and the starting point, and centeredness (a scalar, or one J × K matrix per block) overrides the start. The rule currently supports Gaussian identity, Bernoulli or binomial logit, and Poisson log likelihoods. It cannot yet share a wrapper with totals, ordinary, HSGP or cdar blocks.
sb = SBBRMI(brmi; s2z_groups=[:g], s2z_rho=0.0)
problem = StanBlocks.stan_instantiate(sb.model)
adaptive = adaptive_centering_problem(sb, problem, AutoEnzyme(); s2z_rule=:fisher)
fit = WarmupHMC.adaptive_warmup_mcmc(Xoshiro(1), adaptive; n_draws=2000)Per-group coordinates and WarmupHMC
Deprecated
s2z_coordinates=:groups is unrequested and unvalidated — it keeps working, but only pending further exploration/research. (An Opus model came up with it in a fever dream.) Prefer the default :contrasts.
s2z_coordinates=:groups samples one coordinate per group instead of the J - 1 contrasts: independent cells s_j ~ N(0, tau^(2c_j)), with w = s ./ tau.^c and deviations tau * (w - mean(w)). This is Sean's per-group projected map with one extra dimension: mean(w) is an independent N(0, 1/J) auxiliary that never reaches the likelihood, so the posterior is unchanged. Here s2z_rho is each group's power-interpolation centeredness c (0 = noncentered, 1 = centered; default 0).
Because every group is an independent scalar cell, adaptive_centering_problem gives each group and coefficient its own online centering control, starting from the compiled c:
sb = SBBRMI(brmi; s2z_groups=[:g], s2z_coordinates=:groups)
problem = StanBlocks.stan_instantiate(sb.model)
adaptive = adaptive_centering_problem(sb, problem, AutoEnzyme())
fit = WarmupHMC.adaptive_warmup_mcmc(Xoshiro(1), adaptive;
n_draws=2000, nonlinear_adapt=true)
names = BridgeStan.param_unc_names(problem.model)
recovered = recover_s2z_draws(sb, permutedims(fit.posterior_position), names)For a post-hoc refit, select_s2z_centeredness(sb, draws, names; criterion) scores the same cells from compiled-frame pilot draws, using the same losses as select_total_centeredness. Pass its centeredness to a fresh adaptive_centering_problem and fit with nonlinear_adapt=false, or recompile with s2z_rho=reshape(centeredness, J, K).
Starting from the noncentered frame is the robust default. A fully centered start on a weakly identified scale can spend warmup in the funnel before adaptation moves away from it.
Contrast coordinates and WarmupHMC
With the default contrast coordinates, compile an endpoint frame: s2z_rho=0 (standard-normal contrasts) or s2z_rho=1 (centered contrasts); a vector such as [0, 1] sets the endpoint per coefficient. The same wrapper and selector then treat every free contrast as a scalar cell with zero location and scale tau_k. Interior weights have no per-contrast equivalent and are refused.
These controls belong to contrasts, not groups. Contrast r puts weight r/(r+1) on group r+1 and the remainder on groups 1:r. With unbalanced groups and a weakly identified scale, a contrast that mixes data-rich and data-poor groups can settle between their preferred centerings and leave divergent transitions that per-group coordinates avoid.
S2Z and total cells can share one wrapper, with totals first. Neither can yet be combined with ordinary, HSGP or cdar cells.
Response-level wrappers
| Term | What it does |
|---|---|
mi(y) | brms-style observed/imputed split — missing rows become parameters drawn from the same family |
weighted(y, w) | typed observation weights on the likelihood |
truncated, censored, interval_censored | truncation and censoring — see Likelihoods |
mm's membership weights are part of the random-effect contribution and are a different thing from weighted(...), which scales the likelihood.
Group-local kernels
| Term | What it does |
|---|---|
kernel(args…) do … | a per-subject model cell — the PMX/PKPD surface |
ragged(x, group) | group a flat secondary row axis into the ragged per-subject view a kernel cell slices |
ragged's x may be a linear predictor declared in the same @brm block or a raw flat data column; group names, for every row of that frame, which subject the row belongs to. See the multi-axis population PK kernel for a runnable example whose subject and observation columns have different lengths.
Downstream grouped terms: the reusable-term pattern
A parametric shape with group-varying parameters — a single-peak transient, a saturating sigmoid — belongs in the downstream package whose science needs it, not in BRM core. BRM core promotes such a term only once two or more downstream packages have shipped the same curves; until then the package ships its own term through the supported group-block seams, first-class rather than as an escape hatch. The worked example below is a transient / saturating pair.
The recipe, for a term used in nested predictor position (mu ~ 1 + transient(logt; series)), where it runs on both backends:
Declare the marker in the downstream module (
function transient end).Declare one structured-latent field per parameter group with the fields form of
_sb_term_group_block. Thegroupspec names the grouping column ((; kwarg=:series)reads it from the call'sseries=keyword);prior=:correlated_normaldraws the per-group parameters through BRM's non-centered LKJ block on both backends.Route the nested summand to your emitter:
_sb_is_direct_term(::typeof(transient)) = true.Prepare once, replay frozen:
_brm_prepares_termis true;_brm_prepare_termdelegates to_brm_prepare_structured_termand bakes the axis values the native effect needs into state;_brm_replay_termrefreshes indices and axis values against the fitted levels.Assemble per-row values natively for Turing in
_brm_native_structured_effect— pure Julia math over the block row.Emit the Stan contribution in
_sb_predictor_term!, reusing StanBlocks builtins (e.g.biomarker_time_response) and threading the preallocated block via_sb_find_group_block, so the downstream module ships no custom Stan code. (Extend, never shadow:import BayesianRegressionModels: _sb_term_group_block, _sb_is_direct_term, _brm_prepares_term, _brm_prepare_term, _brm_replay_term, _brm_native_structured_effect, _sb_predictor_term!.)
module DownstreamTerms
using StanBlocks
import BayesianRegressionModels
const BRM = BayesianRegressionModels
import BayesianRegressionModels: _sb_term_group_block, _sb_is_direct_term,
_brm_prepares_term, _brm_prepare_term, _brm_replay_term,
_brm_native_structured_effect, _sb_predictor_term!
function transient end
_sb_is_direct_term(::typeof(transient)) = true
_sb_term_group_block(::typeof(transient)) = (; fields=[
(; name=:transient, n_per_group=3, group=(; kwarg=:series),
prior=:correlated_normal),
])
# Independent Julia reference math, shared by the native effect below.
bump_math(x, loc, ls, mag) = begin
xi = (x - loc) * exp(ls)
s = 1 / (1 + exp(-xi))
sm = 1 / (1 + exp(xi))
exp(log(s) + log(sm)) * mag
end
_brm_prepares_term(::BRM.ExprColumn{typeof(transient)}) = true
function _brm_prepare_term(term::BRM.ExprColumn{typeof(transient)}, target,
context)
base = BRM._brm_prepare_structured_term(term, target, context,
_sb_term_group_block(transient, term))
xkey, xraw = BRM._brm_term_data(
:transient, only(BRM.getargs(term)), context)
BRM._BRMPreparedTerm(base.callable, base.source,
merge(base.state, (; xname=xkey, x=collect(Float64, xraw))),
base.dependencies)
end
function _brm_replay_term(::typeof(transient), training, fresh, context)
fields = map(training.state.fields) do field
raw = context.data[field.source]
merge(field, (; idx=BRM._brm_apply_levels(field.levels, raw)))
end
xraw = context.data[training.state.xname]
BRM._BRMPreparedTerm(training.callable, training.source,
merge(training.state,
(; fields=Tuple(fields), x=collect(Float64, xraw))),
training.dependencies)
end
function _brm_native_structured_effect(
term::BRM._BRMPreparedTerm{typeof(transient)}, block, field)
st = term.state
[bump_math(st.x[i], block[field.idx[i], 1], block[field.idx[i], 2],
block[field.idx[i], 3]) for i in eachindex(field.idx)]
end
function _sb_predictor_term!(stmts, data, ::typeof(transient), t;
target::Symbol, group_block_lookup=Dict(),
kwargs...)
xname, xraw = BRM._sb_inner_data(:transient, only(BRM.getargs(t)))
data[xname] = collect(Float64, xraw)
info = BRM._sb_find_group_block(transient, t, group_block_lookup)
isnothing(info) && error("sbimpl: `transient` found no allocated block")
(; block_name, idx_name) = info
loc = Symbol(:transient_, target, :_loc)
ls = Symbol(:transient_, target, :_ls)
mag = Symbol(:transient_, target, :_mag)
col = Symbol(:transient_, target, :_, xname)
push!(stmts, :($loc = $(block_name)[$(idx_name), 1]))
push!(stmts, :($ls = $(block_name)[$(idx_name), 2]))
push!(stmts, :($mag = $(block_name)[$(idx_name), 3]))
push!(stmts, :($col = biomarker_time_response($xname, $loc, $ls, $mag)))
col
end
end- Fit through
SBBRMI(brmi; mod=DownstreamTerms)orTuringBRMI(brmi). The Stan path lowers to stanc-clean code; the Turing path draws the same non-centered LKJ block natively and assembles bit-exact effects —test/downstream_group_block_term.jlverifies both against independent Julia math.
A second shape follows the same seven steps. A saturating 0-to-1 dose multiplier is two per-group parameters (loc, log_slope) whose Stan emit exps the log-sigmoid builtin and whose native effect is one line of Julia:
function saturating end
_sb_is_direct_term(::typeof(saturating)) = true
_sb_term_group_block(::typeof(saturating)) = (; fields=[
(; name=:saturating, n_per_group=2, group=(; kwarg=:series),
prior=:correlated_normal),
])
sigmoid_math(x, loc, ls) = 1 / (1 + exp(-((x - loc) * exp(ls))))
# ... same prep/replay shape as above, with `saturating` for `transient` ...
function _brm_native_structured_effect(
term::BRM._BRMPreparedTerm{typeof(saturating)}, block, field)
st = term.state
[sigmoid_math(st.x[i], block[field.idx[i], 1], block[field.idx[i], 2])
for i in eachindex(field.idx)]
end
function _sb_predictor_term!(stmts, data, ::typeof(saturating), t;
target::Symbol, group_block_lookup=Dict(),
kwargs...)
xname, xraw = BRM._sb_inner_data(:saturating, only(BRM.getargs(t)))
data[xname] = collect(Float64, xraw)
info = BRM._sb_find_group_block(saturating, t, group_block_lookup)
isnothing(info) && error("sbimpl: `saturating` found no allocated block")
(; block_name, idx_name) = info
loc = Symbol(:saturating_, target, :_loc)
ls = Symbol(:saturating_, target, :_ls)
col = Symbol(:saturating_, target, :_, xname)
push!(stmts, :($loc = $(block_name)[$(idx_name), 1]))
push!(stmts, :($ls = $(block_name)[$(idx_name), 2]))
push!(stmts, :($col = exp(biomarker_dose_response($xname, $loc, $ls))))
col
endShapes compose two ways. Additively, as ordinary predictor summands:
mu ~ 1 + transient(logt; series) + saturating(logd; series)Multiplicatively — baseline plus bump times response — through intercept-free named predictors combined by assignment (write *: the formula layer is element-wise by intent and sbimpl dots it — a literal .* is evaluated at parse time and fails):
brmi = (@brm df begin
sigma ~ Exponential(1)
base ~ Normal(0, 1)
bump ~ 0 + transient(logt; series)
resp ~ 0 + saturating(logd; series)
mu = base + bump * resp
y ~ Normal(mu, sigma)
end)
sb = SBBRMI(brmi; mod=DownstreamTerms)
tb = TuringBRMI(brmi)Each shape owns its own per-group hierarchy (two LKJ blocks here): parameters correlate within a shape, not across shapes. A model that needs the bump and the sigmoid parameters jointly correlated (all six in one covariance) wants one joint term with n_per_group=6 instead, or the shared-bucket + kernel(...) composition; which one fits is decided per model. Term-internal prior statements (sd(mu, transient(logt))) are not addressed yet: the hierarchical scales keep their shared defaults, and naming a new term in a prior address needs core registration alongside _TERM_HEADS.
Verified contract (test/downstream_group_block_term.jl pins all of it on both backends): stanc-clean Stan lowering, Turing effects bit-exact versus independent Julia math with finite joint densities, brm_descriptor, and reprocess / restan_data on fitted group levels. Replaying a group level the fit never saw stays refused by the shared ranef floor — the fitted model has no coordinate for it, on any backend, for core grouped terms alike.
Fixed-one contribution: offset(x)
offset adds x directly to a population-level linear predictor with coefficient one. It never allocates a beta_pop column. The argument may be a raw-data expression, as in an exposure offset, or an already-declared sampled scalar:
Sampled-value offsetoffset_model = (@brm begin
log_ka_pop ~ Normal(-2.08, 1.0)
sigma ~ Exponential(1.0)
eta_ka ~ offset(log_ka_pop) + (1 | p | subject)
concentration ~ Normal(exp(eta_ka), sigma)
end)((;
subject=[1, 1, 2, 2],
concentration=[0.2, 0.4, 0.7, 0.5],
))BRMI:
log_ka_pop ~ Normal(-2.08, 1.0)
sigma ~ Exponential(1.0)
subject: data (eltype=Int64, n=4)
eta_ka ~ offset(log_ka_pop) + (1 | p | subject)
concentration ~ Normal(exp(eta_ka), sigma)SBBRMI with data keys = [:concentration, :n_subject, :n_terms_p_subject, :subject_idx]
emitted @slic body:
begin
b_p_subject ~ ranef_correlated_draws(; group_idx = subject_idx, n_groups = n_subject, n_terms = n_terms_p_subject)
log_ka_pop ~ normal(-2.08, 1.0)
sigma ~ exponential(1.0 ./ 1.0)
r_eta_ka_p_subject = b_p_subject[subject_idx, 1]
eta_ka = log_ka_pop + r_eta_ka_p_subject
concentration ~ normal((exp)(eta_ka), sigma)
endfunctions {
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 n_terms_p_subject;
int n_subject;
int subject_idx_n;
array[subject_idx_n] int subject_idx;
int concentration_n;
vector[concentration_n] concentration;
}
transformed data {
}
parameters {
cholesky_factor_corr[n_terms_p_subject] b_p_subject_L;
vector<lower=0.0>[n_terms_p_subject] b_p_subject_tau;
vector[(n_terms_p_subject * n_subject)] b_p_subject_z_flat;
real log_ka_pop;
real<lower=0.0> sigma;
}
transformed parameters {
matrix[n_terms_p_subject, n_subject] b_p_subject_z = to_matrix(b_p_subject_z_flat, n_terms_p_subject, n_subject);
matrix[n_subject, n_terms_p_subject] b_p_subject = ((diag_pre_multiply(b_p_subject_tau, b_p_subject_L) * b_p_subject_z)');
vector[subject_idx_n] r_eta_ka_p_subject = b_p_subject[subject_idx, 1];
vector[subject_idx_n] eta_ka = (log_ka_pop + r_eta_ka_p_subject);
}
model {
b_p_subject_L ~ lkj_corr_cholesky(1.0);
b_p_subject_tau ~ std_normal();
b_p_subject_z_flat ~ std_normal();
log_ka_pop ~ normal(-2.08, 1.0);
sigma ~ exponential((1.0 ./ 1.0));
concentration ~ normal(exp(eta_ka), sigma);
}
generated quantities {
vector[concentration_n] concentration_likelihood = normal_lpdfs(concentration, exp(eta_ka), sigma);
vector[concentration_n] concentration_gen = normal_vector_rng(concentration_n, exp(eta_ka), sigma);
}Turing unsupported for this BRM example
BRM backend lowering: direct-BRMI `offset(...)` currently requires a pure expression over raw numeric data columnsUse protect for literal data transformations whose resulting column should receive an estimated population coefficient. protect(log_ka_pop) is therefore not an alternative spelling: protect materializes raw-data expressions, while offset preserves a model-value reference and fixes its coefficient to one. Offsets are population-level terms and are rejected inside (... | group) random-effects terms.
Penalized smooth: s(x)
s adds a penalized one-dimensional thin-plate regression spline to a linear predictor. It is available only in the SBBRMI StanBlocks backend.
Penalized smoothx_smooth = collect(range(-2, 2; length=50))
smooth_model = (@brm begin
y ~ Normal(mu, sigma)
mu ~ s(x)
sigma ~ Exponential(1)
end)((; x=x_smooth, y=sin.(x_smooth)))BRMI:
mu ~ s(x)
sigma ~ Exponential(1)
y ~ Normal(mu, sigma)
x: data (eltype=Float64, n=50)SBBRMI with data keys = [:Xnull_x, :Zpen_x, :x, :y]
emitted @slic body:
begin
sigma ~ exponential(1.0 ./ 1)
s_x ~ BayesianRegressionModels._sb_s_generic(; Xnull = Xnull_x, Zpen = Zpen_x)
mu = s_x
y ~ normal(mu, sigma)
endfunctions {
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 Zpen_x_n;
int Xnull_x_m;
int Xnull_x_n;
matrix[Xnull_x_m, Xnull_x_n] Xnull_x;
int Zpen_x_m;
matrix[Zpen_x_m, Zpen_x_n] Zpen_x;
int y_n;
vector[y_n] y;
}
transformed data {
int s_x_n_pen = Zpen_x_n;
}
parameters {
real<lower=0.0> sigma;
vector[2] s_x_b_fixed;
vector<lower=0.0>[1] s_x_sd_pen;
vector[s_x_n_pen] s_x_b_pen_raw;
}
transformed parameters {
vector[s_x_n_pen] s_x_b_pen = (s_x_sd_pen[1] * s_x_b_pen_raw);
vector[Xnull_x_m] s_x = ((Xnull_x * s_x_b_fixed) + (Zpen_x * s_x_b_pen));
vector[Xnull_x_m] mu = s_x;
}
model {
sigma ~ exponential((1.0 ./ 1));
s_x_sd_pen ~ std_normal();
s_x_b_pen_raw ~ std_normal();
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);
}#= line 0 =# Turing.@model(function brm_model(y, terms_mu_1, x)
sigma ~ Distributions.Exponential(1)
term_mu_1 ~ DynamicPPL.to_submodel(BRM.turing_term_model(terms_mu_1, Base.length(y), NamedTuple{$(QuoteNode((:sd,)))}((Distributions.Normal(0.0, 1.0),)), NamedTuple{$(QuoteNode((:x,)))}((x,))))
eta_mu = term_mu_1.effect
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 supported public call has exactly one numeric predictor and no keyword arguments. Training values must be finite and contain at least 10 unique values. Internally, s(x) uses a fixed rank-10 basis: two unpenalized null-space columns for {1, x} and eight penalty-whitened range columns. The range coefficients share a smoothing standard deviation whose default prior is a standard half-normal; sd(<lp|:>, s(x)) ~ Exponential(scale) replaces it (see Term-internal priors). The term contributes this complete smooth directly to the linear predictor, so it does not receive an additional population coefficient.
For prediction or posterior replay on new data, the default reprocess(sb, new_df) and restan_data(sb, new_df) calls evaluate x against the frozen training centers and basis. Passing freeze_constants=false re-estimates the basis from the new data and therefore has fresh-fit rather than prediction semantics.
Difference from bs(...)
Bambi/Formulae formulas such as bs(x, knots=knots) create a deterministic B-spline design matrix whose dimension and knots are controlled by the formula. BayesianRegressionModels' s(x) instead represents a penalized thin-plate smooth with the fixed rank and smoothing prior described above. They are not a one-for-one syntax translation: s(x; k=...), s(x; knots=...), bs(...), and t2(...) are not alternative spellings of this term.
Tensor-product smooth: t2(x, z)
t2 adds a two-margin tensor-product smooth to an SBBRMI linear predictor. Its public defaults follow the brms/mgcv t2 catalogue surface, with Julia tuples for per-margin options:
Tensor-product smootht2_model = (@brm begin
loc ~ 1 + t2(area, yearc;
k=(5, 5), basis=(:cr, :cr), full=false)
rent ~ Normal(loc, sigma)
sigma ~ Exponential(1)
end)((;
area=repeat(collect(1.0:5.0); inner=5),
yearc=repeat(collect(-2.0:2.0); outer=5),
rent=collect(range(0.1, 2.5; length=25)),
))BRMI:
area: data (eltype=Float64, n=25)
yearc: data (eltype=Float64, n=25)
loc ~ 1 + t2(area, yearc; k=(5, 5), basis=(:cr, :cr), full=false)
sigma ~ Exponential(1)
rent ~ Normal(loc, sigma)SBBRMI with data keys = [:Xfixed_t2_loc_area_yearc, :Znr_t2_loc_area_yearc, :Zrn_t2_loc_area_yearc, :Zrr_t2_loc_area_yearc, :area, :rent, :yearc]
emitted @slic body:
begin
X_loc = hcat(rep_vector(1.0, num_elements(area)))
pop_loc ~ popefs(; X = X_loc)
t2_loc_area_yearc ~ BayesianRegressionModels._sb_t2_generic(; Xfixed = Xfixed_t2_loc_area_yearc, Zrr = Zrr_t2_loc_area_yearc, Zrn = Zrn_t2_loc_area_yearc, Znr = Znr_t2_loc_area_yearc)
loc = pop_loc + t2_loc_area_yearc
sigma ~ exponential(1.0 ./ 1)
rent ~ normal(loc, sigma)
endfunctions {
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 area_n;
vector[area_n] area;
int Zrr_t2_loc_area_yearc_n;
int Zrn_t2_loc_area_yearc_n;
int Znr_t2_loc_area_yearc_n;
int Xfixed_t2_loc_area_yearc_m;
int Xfixed_t2_loc_area_yearc_n;
matrix[Xfixed_t2_loc_area_yearc_m, Xfixed_t2_loc_area_yearc_n] Xfixed_t2_loc_area_yearc;
int Zrr_t2_loc_area_yearc_m;
matrix[Zrr_t2_loc_area_yearc_m, Zrr_t2_loc_area_yearc_n] Zrr_t2_loc_area_yearc;
int Zrn_t2_loc_area_yearc_m;
matrix[Zrn_t2_loc_area_yearc_m, Zrn_t2_loc_area_yearc_n] Zrn_t2_loc_area_yearc;
int Znr_t2_loc_area_yearc_m;
matrix[Znr_t2_loc_area_yearc_m, Znr_t2_loc_area_yearc_n] Znr_t2_loc_area_yearc;
int rent_n;
vector[rent_n] rent;
}
transformed data {
matrix[num_elements(area), 1] X_loc = hcat(rep_vector(1.0, num_elements(area)));
int pop_loc_n_covariates = 1;
int t2_loc_area_yearc_n_rr = Zrr_t2_loc_area_yearc_n;
int t2_loc_area_yearc_n_rn = Zrn_t2_loc_area_yearc_n;
int t2_loc_area_yearc_n_nr = Znr_t2_loc_area_yearc_n;
}
parameters {
vector[pop_loc_n_covariates] pop_loc_beta_pop;
vector[3] t2_loc_area_yearc_b_fixed;
vector<lower=0.0>[3] t2_loc_area_yearc_sd_pen;
vector[t2_loc_area_yearc_n_rr] t2_loc_area_yearc_b_rr_raw;
vector[t2_loc_area_yearc_n_rn] t2_loc_area_yearc_b_rn_raw;
vector[t2_loc_area_yearc_n_nr] t2_loc_area_yearc_b_nr_raw;
real<lower=0.0> sigma;
}
transformed parameters {
vector[num_elements(area)] pop_loc = (X_loc * pop_loc_beta_pop);
vector[t2_loc_area_yearc_n_rr] t2_loc_area_yearc_b_rr = (t2_loc_area_yearc_sd_pen[1] * t2_loc_area_yearc_b_rr_raw);
vector[t2_loc_area_yearc_n_rn] t2_loc_area_yearc_b_rn = (t2_loc_area_yearc_sd_pen[2] * t2_loc_area_yearc_b_rn_raw);
vector[t2_loc_area_yearc_n_nr] t2_loc_area_yearc_b_nr = (t2_loc_area_yearc_sd_pen[3] * t2_loc_area_yearc_b_nr_raw);
vector[Xfixed_t2_loc_area_yearc_m] t2_loc_area_yearc = (
(Xfixed_t2_loc_area_yearc * t2_loc_area_yearc_b_fixed) +
(Zrr_t2_loc_area_yearc * t2_loc_area_yearc_b_rr) +
(Zrn_t2_loc_area_yearc * t2_loc_area_yearc_b_rn) +
(Znr_t2_loc_area_yearc * t2_loc_area_yearc_b_nr)
);
vector[num_elements(area)] loc = (pop_loc + t2_loc_area_yearc);
}
model {
pop_loc_beta_pop ~ std_normal();
t2_loc_area_yearc_sd_pen ~ std_normal();
t2_loc_area_yearc_b_rr_raw ~ std_normal();
t2_loc_area_yearc_b_rn_raw ~ std_normal();
t2_loc_area_yearc_b_nr_raw ~ std_normal();
sigma ~ exponential((1.0 ./ 1));
rent ~ normal(loc, sigma);
}
generated quantities {
vector[rent_n] rent_likelihood = normal_lpdfs(rent, loc, sigma);
vector[rent_n] rent_gen = normal_vector_rng(rent_n, loc, sigma);
}#= line 0 =# Turing.@model(function brm_model(y, X_loc, terms_loc_1, area, yearc)
beta_pop ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 1))
eta_loc = X_loc * beta_pop
term_loc_1 ~ DynamicPPL.to_submodel(BRM.turing_term_model(terms_loc_1, Base.length(y), NamedTuple{$(QuoteNode((:sd,)))}(((Distributions.Normal(0.0, 1.0), Distributions.Normal(0.0, 1.0), Distributions.Normal(0.0, 1.0)),)), NamedTuple{$(QuoteNode((:area, :yearc)))}((area, yearc))))
eta_loc = eta_loc + term_loc_1.effect
loc = eta_loc
sigma ~ Distributions.Exponential(1)
begin
for i = Base.eachindex(y)
y[i] ~ Distributions.Normal(loc[i], sigma)
end
end
(; loc = loc, sigma = sigma, response = y)
end)Both predictors must be finite numeric columns with at least the corresponding number of unique values in k. Each k entry must be an integer greater than
- The current implementation accepts only cubic-regression-spline margins
(basis=(:cr, :cr)) and full=false; unsupported values and unknown keywords are rejected while the BRMI is built.
Each marginal basis is split into a two-dimensional null space and a penalty-whitened range space. Their tensor product has three unpenalized null×null columns after the intercept constraint, plus independently scaled range×range, range×null, and null×range blocks. At the default k=(5, 5), those penalized blocks have 9, 6, and 6 coefficients. Each block has its own smoothing scale, defaulting to a standard half-normal and addressable individually as sd(<lp|:>, t2(x, z), <block>) with <block> one of rr, rn, nr (see Term-internal priors). The complete smooth is added directly without an extra population coefficient.
reprocess(sb, new_df) evaluates all four blocks against the frozen training knots, penalty decomposition, and intercept constraint. Passing freeze_constants=false re-estimates them from new_df. Like s, t2 is implemented by both the StanBlocks and Turing backends; it remains unavailable to VBRMI.
Term-internal priors
Some terms own parameters that no coefficient address can reach. s(x)'s smoothing scale, mo(c)'s Dirichlet increments, me(x, sd)'s latent true covariate, interval_censored(x; upper=lloq)'s bounded latent values, and a Gaussian process's length scale and amplitude all live inside the term's own submodel. A dar trajectory likewise owns its persistence and innovation scale. None is a beta_pop column or a grouping-factor margin. They are addressed by naming the term itself in the target slot, in the same head-position grammar the rest of the prior surface uses:
<quantity>(<linear predictor | :>, <term>[, <component>]) ~ <distribution>Term-internal priorsterm_prior_model = (@brm begin
y ~ Normal(mu, sigma)
mu ~ 1 + s(age) + t2(x, z) + mo(dose) + me(w_obs, 0.3) + hsgp(conc; k=5)
sigma ~ Exponential(1)
sd(:, s(age)) ~ Exponential(2) # smoothing scale
sd(mu, t2(x, z), rr) ~ Exponential(3) # one tensor penalty
simplex(:, mo(dose)) ~ Dirichlet(2) # monotonic increments
latent(mu, 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:44.0),
x=repeat(collect(1.0:5.0); inner=5),
z=repeat(collect(-2.0:2.0); outer=5),
dose=repeat(1:5; inner=5),
w_obs=collect(range(1, 4; length=25)),
conc=collect(range(-2, 2; length=25)),
y=collect(range(-1, 1; length=25)),
))BRMI:
mu ~ 1 + s(age) + t2(x, z) + mo(dose) + me(w_obs, 0.3) + hsgp(conc; k=5)
sigma ~ Exponential(1)
y ~ Normal(mu, sigma)
age: data (eltype=Float64, n=25)
x: data (eltype=Float64, n=25)
z: data (eltype=Float64, n=25)
dose: data (eltype=Int64, n=25)
w_obs: data (eltype=Float64, n=25)
conc: data (eltype=Float64, n=25)
effect(term_sd, s(age), :) ~ Exponential(2)
effect(term_sd, t2(x,z), mu, rr) ~ Exponential(3)
effect(term_simplex, mo(dose), :) ~ Dirichlet(2)
effect(term_latent, me(w_obs), mu) ~ Normal(0, 5)
effect(term_length_scale, hsgp(conc), :) ~ Uniform(0.84, 2)
effect(term_sd, hsgp(conc), :) ~ Normal(0, 0.5)SBBRMI with data keys = [:PHI_hsgp_conc, :Xfixed_t2_mu_x_z, :Xnull_age, :Znr_t2_mu_x_z, :Zpen_age, :Zrn_t2_mu_x_z, :Zrr_t2_mu_x_z, :age, :conc, :dose, :dose_idx, :omega2_hsgp_conc, :rho_lower_hsgp_conc, :sd_w_obs, :w_obs, :x, :y, :z]
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_t2_generic_configured_1 = Base.merge(BayesianRegressionModels._sb_t2_generic, quote
sd_pen ~ (ValueFamily(brm_vector_prior_7e7ad5ab57b1e50c))(0.3333333333333333, 0.0, 1.0, 0.0, 1.0; n = 3)
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
sigma ~ exponential(1.0 ./ 1)
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)
t2_mu_x_z ~ _sb_t2_generic_configured_1(; Xfixed = Xfixed_t2_mu_x_z, Zrr = Zrr_t2_mu_x_z, Zrn = Zrn_t2_mu_x_z, Znr = Znr_t2_mu_x_z)
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 + t2_mu_x_z + hsgp_conc
y ~ normal(mu, sigma)
endfunctions {
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);
}
// value UDF brm_vector_prior_7e7ad5ab57b1e50c_lpdf
real brm_vector_prior_7e7ad5ab57b1e50c_lpdf(
vector x,
real arg_1,
real arg_2,
real arg_3,
real arg_4,
real arg_5
) {
if((x[1] < 0.0)) {
return negative_infinity();
}
if((x[2] < 0.0)) {
return negative_infinity();
}
if((x[3] < 0.0)) {
return negative_infinity();
}
return (
(exponential_lpdf(x[1] | arg_1) + normal_lpdf(x[2] | arg_2, arg_3)) +
normal_lpdf(x[3] | arg_4, arg_5)
);
}
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 Zrr_t2_mu_x_z_n;
int Zrn_t2_mu_x_z_n;
int Znr_t2_mu_x_z_n;
int Xfixed_t2_mu_x_z_m;
int Xfixed_t2_mu_x_z_n;
matrix[Xfixed_t2_mu_x_z_m, Xfixed_t2_mu_x_z_n] Xfixed_t2_mu_x_z;
int Zrr_t2_mu_x_z_m;
matrix[Zrr_t2_mu_x_z_m, Zrr_t2_mu_x_z_n] Zrr_t2_mu_x_z;
int Zrn_t2_mu_x_z_m;
matrix[Zrn_t2_mu_x_z_m, Zrn_t2_mu_x_z_n] Zrn_t2_mu_x_z;
int Znr_t2_mu_x_z_m;
matrix[Znr_t2_mu_x_z_m, Znr_t2_mu_x_z_n] Znr_t2_mu_x_z;
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 t2_mu_x_z_n_rr = Zrr_t2_mu_x_z_n;
int t2_mu_x_z_n_rn = Zrn_t2_mu_x_z_n;
int t2_mu_x_z_n_nr = Znr_t2_mu_x_z_n;
int hsgp_conc_n_basis = omega2_hsgp_conc_m;
int hsgp_conc_n_axes = omega2_hsgp_conc_n;
}
parameters {
real<lower=0.0> sigma;
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;
vector[3] t2_mu_x_z_b_fixed;
vector<lower=0.0>[3] t2_mu_x_z_sd_pen;
vector[t2_mu_x_z_n_rr] t2_mu_x_z_b_rr_raw;
vector[t2_mu_x_z_n_rn] t2_mu_x_z_b_rn_raw;
vector[t2_mu_x_z_n_nr] t2_mu_x_z_b_nr_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[t2_mu_x_z_n_rr] t2_mu_x_z_b_rr = (t2_mu_x_z_sd_pen[1] * t2_mu_x_z_b_rr_raw);
vector[t2_mu_x_z_n_rn] t2_mu_x_z_b_rn = (t2_mu_x_z_sd_pen[2] * t2_mu_x_z_b_rn_raw);
vector[t2_mu_x_z_n_nr] t2_mu_x_z_b_nr = (t2_mu_x_z_sd_pen[3] * t2_mu_x_z_b_nr_raw);
vector[Xfixed_t2_mu_x_z_m] t2_mu_x_z = (
(Xfixed_t2_mu_x_z * t2_mu_x_z_b_fixed) +
(Zrr_t2_mu_x_z * t2_mu_x_z_b_rr) +
(Zrn_t2_mu_x_z * t2_mu_x_z_b_rn) +
(Znr_t2_mu_x_z * t2_mu_x_z_b_nr)
);
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 + t2_mu_x_z + hsgp_conc);
}
model {
sigma ~ exponential((1.0 ./ 1));
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();
t2_mu_x_z_sd_pen ~ brm_vector_prior_7e7ad5ab57b1e50c(0.3333333333333333, 0.0, 1.0, 0.0, 1.0);
t2_mu_x_z_b_rr_raw ~ std_normal();
t2_mu_x_z_b_rn_raw ~ std_normal();
t2_mu_x_z_b_nr_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, sigma);
}
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, sigma);
vector[y_n] y_gen = normal_vector_rng(y_n, mu, sigma);
}#= line 0 =# Turing.@model(function brm_model(y, X_mu, terms_mu_1, age, terms_mu_2, x, z, terms_mu_3, dose, terms_mu_4, w_obs, terms_mu_5, conc)
sigma ~ Distributions.Exponential(1)
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((:sd,)))}(((Distributions.Exponential(3), Distributions.Normal(0.0, 1.0), Distributions.Normal(0.0, 1.0)),)), NamedTuple{$(QuoteNode((:x, :z)))}((x, z))))
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((:simplex,)))}((Distributions.Dirichlet(Base.vect(2.0, 2.0, 2.0, 2.0)),)), NamedTuple{$(QuoteNode((:dose,)))}((dose,))))
eta_mu = eta_mu + term_mu_3.effect
term_mu_4 ~ DynamicPPL.to_submodel(BRM.turing_term_model(terms_mu_4, Base.length(y), NamedTuple{$(QuoteNode((:latent,)))}((Distributions.Normal(0, 5),)), NamedTuple{$(QuoteNode((:w_obs,)))}((w_obs,))))
eta_mu = eta_mu + term_mu_4.effect
term_mu_5 ~ DynamicPPL.to_submodel(BayesianRegressionModelsTuringExt._brm_turing_hsgp_ncp_model(terms_mu_5, 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_5.effect
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)| head | term | parameter | default |
|---|---|---|---|
sd | s(x) | smoothing SD | half-standard-normal |
sd | t2(x, z) | one of the rr / rn / nr penalty SDs | half-standard-normal |
simplex | mo(c), mo1(c) | Dirichlet concentration | Dirichlet(1) |
latent | me(x, sd) | latent true covariate | Normal(0, 1) |
latent | interval_censored(x; upper=lloq) | latent covariate on BLOQ rows | truncated Normal(0, 1) |
length_scale | gp(x…), hsgp(x…) | GP length scale rho | LogNormal(0, 1) |
sd | gp(x…), hsgp(x…) | GP marginal amplitude sigma | LogNormal(0, 1) |
ar | dar(time) | bounded persistence beta | truncated Normal(0.5, 0.2) on [0, 1] |
sd | dar(time) | innovation scale sigma | half-Normal(0, 0.2) |
Spell the term the way the formula does, minus numeric and keyword arguments. me(w_obs, 0.3) is addressed as me(w_obs) and t2(x, z; k=(5,5)) as t2(x, z) — the address names a term, not a call. A : in the predictor slot is THE DEFAULT: it is the base layer that a statement naming a concrete predictor overrides, exactly as on the coefficient and grouping-factor surfaces. Specificity counts concrete slots, so an exact tie is an error rather than a silent winner.
An address that cannot be honoured is refused by name while the model is lowered — an unknown term, a term with no such parameter, a t2 component slot that is missing or names no penalty block, or two terms in one predictor spelled identically. term_priors(brmi) returns the captured statements (class, term, predictor, component, family, arguments, keywords, expression) for formula-level provenance.
Only model-scale quantities are exposed
The standardized raw innovations these submodels sample — b_pen_raw, z, beta_raw — stay iid standard normal and are deliberately NOT configurable. For a smooth, Cov(f_pen | sd) = sd² Zpen Zpen', so a scale on the raw coefficients would simply duplicate the smoothing SD and confound the two. The configurable parameters are the ones that mean something on the model's own scale.
Current limits
sdons(x)/t2(x, z)acceptsExponential(scale)only. In contrast, an addressed shared grouping-factor SD also accepts the zero-centered half-Normal spellingNormal(0, scale).Distributions.Exponentialis scale-parameterized while Stan'sexponential_lpdftakes a rate; BRM performs the conversion.simplexacceptsDirichlet(a)(one concentration, broadcast over every increment) orDirichlet(a₁, …, a_{K-1})(one per increment of aK-level factor). The dimension comes from the data, so theDirichlet(dimension, concentration)spelling used for a standalone simplex parameter is not accepted here.latentacceptsNormal(location, scale). The observation likelihoodx_obs ~ Normal(x_true, sd)is never configurable.length_scaleandsdon agp/hsgpterm acceptLogNormal,InverseGamma,Gamma,Exponential,NormalandUniform, with numeric hyperparameters. Both parameters are strictly positive, so every family exceptUniformis emitted truncated at zero (Normal(0, s)is therefore a half-normal);Uniform(a, b)additionally declares<lower=a, upper=b>so the declaration and the density's support agree. An anisotropic term (iso=false) has one length scale per axis and the statement sets all of them; the isotropic form has a single shared one. Withby=gthe term keeps one shared length scale and amplitude unless hyper-predictor statements predict them per group (Hyper-predictors).ar(:, dar(time))acceptsNormal,Beta, orUniform; every declaration stays within[0, 1].sd(:, dar(time))accepts the same positive-scale families as a GP amplitude. The olderar(time)term's transformed autocorrelation still has no prior address.
The backend compatibility floor for configured gp / hsgp term priors is StanBlocks 10529af04d42a330df383864059c2b61a11d9480. These statements splice a configured SlicModel value into BRM's generated model body. Earlier StanBlocks revisions, including c30d3a158ae6c996dee2423023ab6b35d2756fc9, trace that value before its PHI / omega2 keyword data are bound and fail with Could not find omega2 .... Co-pin BRM revisions containing this surface with 10529af or later. Because StanBlocks is unregistered and both revisions identify as version 0.1.5, the commit SHA—not the package version—is the effective compatibility check.
Omitting an explicit HSGP prior is not a model-preserving workaround. In particular, removing length_scale(:, hsgp(x)) ~ LogNormal(0, 1) restores the default approximation-validity floor described below, changing the parameter support and the emitted Stan program.
Hyper-predictors
A bare grouped term shares one length scale and one amplitude across all groups — only the basis weights vary per group. When each group needs its own smoothness or amplitude (per-biomarker hypers), predict the hypers with hyper-predictor statements: log(length_scale(hsgp(x))) ~ 1 + (1 | g) and log(sd(hsgp(x))) ~ 1 + (1 | g). Both hypers use a log link, and each group evaluates its own linear predictor to eta_rho[g] / eta_sigma[g], consumed as rho_g = max(exp(eta_rho[g]), rho_lower) and sigma_g = exp(eta_sigma[g]).
The hyper linear predictor accepts
1and(1 | g)only. Smooths are refused, population slopes need level-grid covariate semantics that are not decided yet, and the random effect must be exactly(1 | g)over the addressed term's ownby=grouping. An ungrouped term takes an intercept-only scalar (~ 1); a random effect without grouping is refused.Defaults reproduce today's predictive prior: the intercept carries
Normal(0, 1)on the log scale (that isLogNormal(0, 1)on the natural scale), the random-effect SD carries half-Normal(0, 1), and the group deviations are non-centered. An explicitlength_scale(:, hsgp(x)) ~ .../sd(:, hsgp(x)) ~ ...statement retargets from the shared scalar to the hyper-LP intercept and sets its own support, exactly as term priors do.The approximation-validity floor (next section) applies per group: each
rho_gis maximised with the same domain-derived floor the bare term declares, so no group can silently leave the kernel the basis approximates.Scope is deliberately narrow:
hsgponly (gpis refused by name), isotropic (iso=true) only, and neither periodic nor model-derived axes. Lowering is SBBRMI-only.Posterior addressing: the sampled hyper coefficients resolve through
brm_term_coordinatesunder the:length_scale_intercept/:length_scale_ranef_sd/:length_scale_ranef_zroles (and the:sd_*mirrors) — one coordinate for the intercept and SD, one per group for the deviations. The:length_scale/:sdroles themselves name no sampled carrier on a predicted hyper and redirect to these roles; the per-group values are deterministic transforms with no role.After fitting,
hsgp_boundary_checkreports per-group posterior-mean-rhoover the domain margin (flagged at 1.0) and the per-group probability the hyper value sits below the validity floor — the executable form of the sizing guidance below.
Hyper recovery needs a wider domain than latent recovery
A good latent fit does not imply the hypers recover: fitting the smooth only needs the basis to span the data, while hyper posteriors need the HSGP prior covariance to match the kernel on the data grid — a strictly stronger requirement. The Dirichlet boundary pins the prior toward zero within roughly one length scale of ±L, so when the plausible rho approaches the domain margin (c-1)·max|x-mean(x)| the prior variance collapses at the data edges and the rho posterior drags low with tight, confident-wrong credible intervals (snag hyper-predictor-d453ecf9: truth rho=2.1 recovered as 0.76±0.12 at c=1.5, k=10 on a 6-point grid, while the latent RMSE was 0.07).
Size (c, k) for hyper recovery from the plausible hyper range before fitting:
Margin clears the largest plausible
rho:(c-1)·max|x-mean(x)| ≳ 2·rho_max. The defaultc=1.5leaves a margin of half the data half-range; oncerho_maxapproaches that margin, widenc.Floor clears the smallest plausible
rho: raisingcraisesLand the floor(4L/π)·√(log(100)/(k²-1))with it, so raisekalongside until the floor sits belowrho_minwith headroom (targetrho_min/2).Widen until the posterior stabilises: intermediate domains can overshoot (same case:
c=2.5recovers3.35,c=3.0recovers3.86,c=4.0recovers2.75against the exact-GP2.71). Do not stop at the firstcthat moves the posterior.
Two further posterior signatures to read correctly:
Mass with
etabelowlog(rho_lower)is likelihood-flat — the floor binds bymax()inside the basis function, sorho_vec = exp(eta)reports prior mass there as posterior. CheckP(rho_g < rho_lower); if substantial, widenkto lower the floor and re-fit.With few groups the half-
Normal(0, 1)hyper SD shrinks group hypers together (expected hierarchical pooling, not a bug). A chain stuck at hugerhowith a vanished smooth is therho→∞degenerate tail all HSGP hypers share — checkRhat: it is non-convergence, not information.
hsgp bounds its length scale by default
An HSGP with k basis functions over a domain of half-width L = c·max|x − mean(x)| stops approximating the kernel it was asked for once the length scale falls below
(4L/π)·√(log(100)/(k² − 1))and the failure is silent: the model transpiles, passes stanc, samples, and returns finite draws that simply are not that Gaussian process. On the default LogNormal(0, 1) at L = 1.5 that region holds 42.9 % of the prior mass at k = 5, 18.8 % at k = 10 and 5.7 % at k = 20.
Every hsgp term therefore declares rho with that floor as its lower bound. The density is unchanged — only the support moves. Exact gp has no basis truncation and is untouched.
The rule behind the number is that the k-th basis function's spectral amplitude has fallen to 1/100 of the first's. A cov=:periodic term applies the same rule to its own weights, q_j ∝ √(e^{−a} I_j(a)) with a = 1/ρ²: its floor is the length scale solving I_k(a)/I_1(a) = 100⁻², found by bisection on the scaled Bessel functions. It depends on k only — there is no data-derived L — so every reprocess reproduces it exactly.
The floor is emitted as data (rho_lower_hsgp_<axes>), not as a literal, because L comes from the covariate: reprocess(sb, df2; freeze_constants=false) re-derives it alongside PHI / omega2 while the Stan source stays byte-identical. Isotropic spellings share one rho across axes and take the strictest per-axis floor; iso=false bounds each axis separately. k = 1 puts the formula at infinity, so that degenerate basis stays unbounded.
An explicit length_scale statement replaces the whole declaration, floor included:
mu ~ hsgp(x; k=5, c=1.5) # real<lower=rho_lower_hsgp_x>
length_scale(:, hsgp(x)) ~ LogNormal(0, 1) # real<lower=0.0> — floor goneThat is deliberate. It is how a pre-floor posterior is reproduced, and it is what keeps Uniform(a, b) self-consistent — an unconditional floor would overwrite lower=a while the density stayed uniform(a, b), leaving the declaration and the support disagreeing. The consequence is that the guarantee is default-on, not absolute: an override with mass below the floor restores the silent approximation error, with no warning. Either keep the floor in the overriding bounds (Uniform(0.84, 2)) or accept the error knowingly.
This changes the posterior of existing unedited hsgp models
Models that never named length_scale sampled an unbounded rho before this default landed. If you version models by their emitted Stan plus data, treat the change as a new model version rather than a refresh of the old one.
R²-induced variance decomposition: effect(lp, :) ~ r2d2(...)
r2d2 replaces the independent per-coefficient priors on a linear predictor's population block with a joint prior induced by a prior on that predictor's coefficient of determination. It is a separate statement in the formula block — addressed with effect, never with a per-column effect(lp, coef) ~ Normal(...) — and is implemented by both the SBBRMI and TuringBRMI backends.
R2D2 population and group decompositionr2d2_model = (@brm begin
log_CL ~ 1 + wt + age + (1 | p | subject)
log_V ~ 1 + wt + (1 | p | subject)
effect(log_CL, :) ~ r2d2(R2=Beta(1, 1), tau_bsv=0.5)
effect(log_V, :) ~ r2d2(R2=Beta(2, 3), tau_bsv=0.25)
conc ~ Normal(exp(log_CL - log_V) * time, 1)
end)((;
wt=[55.0, 65.0, 75.0, 85.0], age=[21.0, 38.0, 55.0, 29.0],
subject=[1, 1, 2, 2], time=[0.5, 1.0, 1.5, 2.0],
conc=[0.2, 0.4, 0.3, 0.1],
))BRMI:
wt: data (eltype=Float64, n=4)
age: data (eltype=Float64, n=4)
subject: data (eltype=Int64, n=4)
log_CL ~ 1 + wt + age + (1 | p | subject)
log_V ~ 1 + wt + (1 | p | subject)
effect(log_CL, :) ~ r2d2(; R2=Beta(1, 1), tau_bsv=0.5)
effect(log_V, :) ~ r2d2(; R2=Beta(2, 3), tau_bsv=0.25)
time: data (eltype=Float64, n=4)
conc ~ Normal((exp(-(log_CL, log_V)) * time), 1)SBBRMI with data keys = [:age, :conc, :n_subject, :n_terms_p_subject, :r2d2_log_CL_alpha, :r2d2_log_CL_beta_loc, :r2d2_log_CL_fallback, :r2d2_log_CL_share_idx, :r2d2_log_CL_tau_bsv, :r2d2_log_V_alpha, :r2d2_log_V_beta_loc, :r2d2_log_V_fallback, :r2d2_log_V_share_idx, :r2d2_log_V_tau_bsv, :subject_idx, :time, :wt]
emitted @slic body:
begin
r2d2_log_CL_R2 ~ beta(1.0, 1.0)
r2d2_log_CL_phi ~ dirichlet(r2d2_log_CL_alpha)
r2d2_log_V_R2 ~ beta(2.0, 3.0)
r2d2_log_V_phi ~ dirichlet(r2d2_log_V_alpha)
b_p_subject_r2d2_tau = [sqrt((1.0 - r2d2_log_CL_R2) * r2d2_log_CL_tau_bsv ^ 2), sqrt((1.0 - r2d2_log_V_R2) * r2d2_log_V_tau_bsv ^ 2)]
b_p_subject ~ ranef_correlated_draws_r2d2(; group_idx = subject_idx, n_groups = n_subject, n_terms = n_terms_p_subject, tau = b_p_subject_r2d2_tau, lkj_eta = 1.0)
X_log_CL = hcat(rep_vector(1.0, num_elements(wt)), wt, age)
r2d2_log_CL_varx = brm_col_variances(X_log_CL, (dims(X_log_CL))[1], (dims(X_log_CL))[2])
r2d2_log_CL_beta_scale = brm_r2d2_scale(r2d2_log_CL_share_idx, r2d2_log_CL_fallback, r2d2_log_CL_varx, r2d2_log_CL_phi, r2d2_log_CL_R2, r2d2_log_CL_tau_bsv, (dims(X_log_CL))[2], 2)
pop_log_CL ~ _popefs_normal(; X = X_log_CL, beta_loc = r2d2_log_CL_beta_loc, beta_scale = r2d2_log_CL_beta_scale)
r_log_CL_p_subject = b_p_subject[subject_idx, 1]
log_CL = pop_log_CL + r_log_CL_p_subject
X_log_V = hcat(rep_vector(1.0, num_elements(wt)), wt)
r2d2_log_V_varx = brm_col_variances(X_log_V, (dims(X_log_V))[1], (dims(X_log_V))[2])
r2d2_log_V_beta_scale = brm_r2d2_scale(r2d2_log_V_share_idx, r2d2_log_V_fallback, r2d2_log_V_varx, r2d2_log_V_phi, r2d2_log_V_R2, r2d2_log_V_tau_bsv, (dims(X_log_V))[2], 1)
pop_log_V ~ _popefs_normal(; X = X_log_V, beta_loc = r2d2_log_V_beta_loc, beta_scale = r2d2_log_V_beta_scale)
r_log_V_p_subject = b_p_subject[subject_idx, 2]
log_V = pop_log_V + r_log_V_p_subject
conc ~ normal((exp)((-)(log_CL, log_V)) .* time, 1)
endfunctions {
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);
}
vector brm_col_variances(
matrix X,
int m,
int n
) {
vector[n] rv = rep_vector(1.0, n);
for(k in 1:n) {
rv[k] = variance(col(X, k));
}
return rv;
}
vector brm_r2d2_scale(
array[] int share_idx,
vector fallback,
vector varx,
vector phi,
real R2,
real tau_bsv,
int n,
int k
) {
vector[n] rv = rep_vector(0.0, n);
for(j in 1:n) {
if((share_idx[j] == 0)) {
rv[j] = fallback[j];
} else {
rv[j] = sqrt(((phi[share_idx[j]] * R2 * (tau_bsv ^ 2)) / varx[j]));
}
}
return rv;
}
vector normal_lpdfs(
vector obs,
vector loc,
int scale
) {
return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
vector x1,
vector x2,
int 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,
int 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,
int b
) {
int n = anontok__1;
if((n == 0)) {
vector[n] rv;
return rv;
} else {
return to_vector(normal_rng(a, b));
}
}
}
data {
int r2d2_log_CL_alpha_n;
vector[r2d2_log_CL_alpha_n] r2d2_log_CL_alpha;
int r2d2_log_V_alpha_n;
vector[r2d2_log_V_alpha_n] r2d2_log_V_alpha;
real r2d2_log_CL_tau_bsv;
real r2d2_log_V_tau_bsv;
int n_terms_p_subject;
int n_subject;
int age_n;
int wt_n;
vector[wt_n] wt;
vector[age_n] age;
int r2d2_log_CL_share_idx_n;
array[r2d2_log_CL_share_idx_n] int r2d2_log_CL_share_idx;
int r2d2_log_CL_fallback_n;
vector[r2d2_log_CL_fallback_n] r2d2_log_CL_fallback;
int r2d2_log_CL_beta_loc_n;
vector[r2d2_log_CL_beta_loc_n] r2d2_log_CL_beta_loc;
int subject_idx_n;
array[subject_idx_n] int subject_idx;
int r2d2_log_V_share_idx_n;
array[r2d2_log_V_share_idx_n] int r2d2_log_V_share_idx;
int r2d2_log_V_fallback_n;
vector[r2d2_log_V_fallback_n] r2d2_log_V_fallback;
int r2d2_log_V_beta_loc_n;
vector[r2d2_log_V_beta_loc_n] r2d2_log_V_beta_loc;
int conc_n;
vector[conc_n] conc;
int time_n;
vector[time_n] time;
}
transformed data {
matrix[age_n, (2 + 1)] X_log_CL = hcat(rep_vector(1.0, num_elements(wt)), wt, age);
vector[(2 + 1)] r2d2_log_CL_varx = brm_col_variances(X_log_CL, age_n, (2 + 1));
int pop_log_CL_n_covariates = (2 + 1);
matrix[wt_n, 2] X_log_V = hcat(rep_vector(1.0, num_elements(wt)), wt);
vector[2] r2d2_log_V_varx = brm_col_variances(X_log_V, wt_n, 2);
int pop_log_V_n_covariates = 2;
}
parameters {
real<lower=0, upper=1> r2d2_log_CL_R2;
simplex[r2d2_log_CL_alpha_n] r2d2_log_CL_phi;
real<lower=0, upper=1> r2d2_log_V_R2;
simplex[r2d2_log_V_alpha_n] r2d2_log_V_phi;
cholesky_factor_corr[n_terms_p_subject] b_p_subject_L;
vector[(n_terms_p_subject * n_subject)] b_p_subject_z_flat;
vector[pop_log_CL_n_covariates] pop_log_CL_beta_pop;
vector[pop_log_V_n_covariates] pop_log_V_beta_pop;
}
transformed parameters {
vector[2] b_p_subject_r2d2_tau = [
sqrt(((1.0 - r2d2_log_CL_R2) * (r2d2_log_CL_tau_bsv ^ 2))),
sqrt(((1.0 - r2d2_log_V_R2) * (r2d2_log_V_tau_bsv ^ 2)))
]';
matrix[n_terms_p_subject, n_subject] b_p_subject_z = to_matrix(b_p_subject_z_flat, n_terms_p_subject, n_subject);
matrix[n_subject, n_terms_p_subject] b_p_subject = ((diag_pre_multiply(b_p_subject_r2d2_tau, b_p_subject_L) * b_p_subject_z)');
vector[(2 + 1)] r2d2_log_CL_beta_scale = brm_r2d2_scale(
r2d2_log_CL_share_idx,
r2d2_log_CL_fallback,
r2d2_log_CL_varx,
r2d2_log_CL_phi,
r2d2_log_CL_R2,
r2d2_log_CL_tau_bsv,
(2 + 1),
2
);
vector[age_n] pop_log_CL = (X_log_CL * pop_log_CL_beta_pop);
vector[subject_idx_n] r_log_CL_p_subject = b_p_subject[subject_idx, 1];
vector[age_n] log_CL = (pop_log_CL + r_log_CL_p_subject);
vector[2] r2d2_log_V_beta_scale = brm_r2d2_scale(
r2d2_log_V_share_idx,
r2d2_log_V_fallback,
r2d2_log_V_varx,
r2d2_log_V_phi,
r2d2_log_V_R2,
r2d2_log_V_tau_bsv,
2,
1
);
vector[wt_n] pop_log_V = (X_log_V * pop_log_V_beta_pop);
vector[subject_idx_n] r_log_V_p_subject = b_p_subject[subject_idx, 2];
vector[wt_n] log_V = (pop_log_V + r_log_V_p_subject);
}
model {
r2d2_log_CL_R2 ~ beta(1.0, 1.0);
r2d2_log_CL_phi ~ dirichlet(r2d2_log_CL_alpha);
r2d2_log_V_R2 ~ beta(2.0, 3.0);
r2d2_log_V_phi ~ dirichlet(r2d2_log_V_alpha);
b_p_subject_L ~ lkj_corr_cholesky(1.0);
b_p_subject_z_flat ~ std_normal();
pop_log_CL_beta_pop ~ normal(r2d2_log_CL_beta_loc, r2d2_log_CL_beta_scale);
pop_log_V_beta_pop ~ normal(r2d2_log_V_beta_loc, r2d2_log_V_beta_scale);
conc ~ normal((exp((log_CL - log_V)) .* time), 1);
}
generated quantities {
vector[conc_n] conc_likelihood = normal_lpdfs(conc, (exp((log_CL - log_V)) .* time), 1);
vector[conc_n] conc_gen = normal_vector_rng(conc_n, (exp((log_CL - log_V)) .* time), 1);
}#= line 0 =# Turing.@model(function brm_model(y, X_log_CL, X_log_V, group_effects_log_V_1, group_effects_log_CL_1, time)
r2d2_log_CL ~ DynamicPPL.to_submodel(BayesianRegressionModelsTuringExt._brm_r2d2_population(X_log_CL, $(QuoteNode([0, 1, 2])), Distributions.Beta(1.0, 1.0), 1.0, 0.5, [Distributions.Normal(0, 1), Distributions.Normal(0, 1), Distributions.Normal(0, 1)]))
beta_pop = r2d2_log_CL.beta
eta_log_CL = X_log_CL * beta_pop
group_effect_1 = Base.zeros(Base.length(y))
eta_log_CL = eta_log_CL + group_effect_1
r2d2_log_V ~ DynamicPPL.to_submodel(BayesianRegressionModelsTuringExt._brm_r2d2_population(X_log_V, $(QuoteNode([0, 1])), Distributions.Beta(2.0, 3.0), 1.0, 0.25, [Distributions.Normal(0, 1), Distributions.Normal(0, 1)]))
beta_pop_log_V = r2d2_log_V.beta
eta_log_V = X_log_V * beta_pop_log_V
group_effect_2 = Base.zeros(Base.length(y))
eta_log_V = eta_log_V + group_effect_2
shared_group_1 ~ DynamicPPL.to_submodel(BayesianRegressionModelsTuringExt._brm_shared_group_effect_model((group_effects_log_V_1, group_effects_log_CL_1), (nothing, nothing), (r2d2_log_V.residual_scale, r2d2_log_CL.residual_scale)))
eta_log_V = eta_log_V + shared_group_1.effects[1]
eta_log_CL = eta_log_CL + shared_group_1.effects[2]
log_CL = eta_log_CL
log_V = eta_log_V
begin
for i = Base.eachindex(y)
y[i] ~ Distributions.Normal(Base.exp(log_CL[i] - log_V[i]) * time[i], 1)
end
end
(; log_V = log_V, log_CL = log_CL, response = y)
end)effect(lp, :) addresses every population coefficient of lp at once, and effect(:, :) addresses them across every predictor — which for an r2d2 decomposition resolves only when the model has exactly one population predictor, since the decomposition is per-predictor. An r2d2 statement is deliberately invisible to effect_priors — it names no single labelled column — and is read back with r2d2_priors instead.
What it emits
Writing tau_bsv for the predictor's total scale, the decomposition is
R2 ~ Beta(a, b) # parameter
phi ~ Dirichlet(alpha) # simplex
beta_scale[k] = sqrt(phi[k] * R2 * tau_bsv^2 / Var(x_k)) # transformed
tau_resid = sqrt((1 - R2) * tau_bsv^2) # transformedbeta_scale is injected into the shipped _popefs_normal seam, so the population block still samples as beta_pop ~ normal(beta_loc, beta_scale); tau_resid becomes the random effect's standard deviation. The column variances Var(x_k) depend on data alone and hoist to Stan's transformed-data block.
The intercept is excluded from the simplex. Its design column is constant, so Var(x_k) = 0 there, and an intercept is a location rather than explained variance. A predictor whose only non-intercept column count is one emits a simplex[1]; it is deterministically [1] and adds no sampler dimensions.
The random effect plays the residual role
The inciting shape is population PK: log_CL is a latent per-subject parameter with no residual term of its own, so its subject random effect is the unexplained half of tau_bsv^2. That is why r2d2 derives the random-effect SD rather than sampling it, and why an sd(...) statement on the same block is rejected — the decomposition already determines those scales.
Keywords
| keyword | default | meaning |
|---|---|---|
R2 | required | Beta(a, b) prior on the coefficient of determination |
tau_bsv | sampled half-standard-normal | the predictor's total scale; pass a number to fix it |
alpha | 1 | Dirichlet concentration over the non-intercept columns |
tau_bsv has no data-derived default. A latent per-subject parameter has no observed response to derive a scale from, so an omitted tau_bsv becomes a sampled real<lower=0> with a half-standard-normal prior. Fix it whenever you have a defensible scale — it is the quantity the whole decomposition is relative to.
Current limits
All of these fail loudly rather than silently sampling something else:
The random effect playing the residual role must be a single intercept —
(1 | g)or one margin per predictor inside a(1 | ID | g)bucket. Splitting(1 - R2) * tau_bsv^2over several margins needs a second simplex that this decomposition does not build.A shared
(… | ID | g)bucket is all-or-nothing: either every margin's predictor carries anr2d2statement or none does.cor(:, ID)composes with a shared-ID decomposition and controls its LKJeta; the marginal scales remain derived.Non-centred shared-ID
resample_groupsreplay is supported. Adaptive centering (centered_groups), plain-group R2D2 resampling, stratifiedgr(g, by=b)groups, andmm(...)multi-membership terms remain unsupported.A column that also carries its own
effect(lp, coef) ~ Normal(loc, scale)statement is dropped from the simplex and keeps that explicit prior. The default-layer spellingeffect(:, coef) ~ Normal(...)excludes that column in every predictor it reaches, exactly like the predictor-specific one; the remaining columns and the random-effect residual are still decomposed.Excluding every non-intercept column of an
r2d2-scoped predictor this way is refused. There is then nothing left to allocate, and the only consistent emission would dropR2/phiand fix the random-effect scale at the baretau_bsvwith no prior — a silently different model. Either keep at least one column unaddressed, or move the decomposition onto the random-effect scale withsd(lp, ID) ~ r2d2(reference_scale=...)(the random-effect R2D2M2/ICC form below composes with per-column Normal priors; a shared bucket is all-or-nothing, so switch the whole bucket). A predictor with no non-intercept population column at all (log_ka ~ 1 + (1 | p | g), forced into anr2d2statement by the all-or-nothing rule) is the one legitimate zero-share shape: nothing to explain, so the wholetau_bsvis its random-effect scale.
Random-effect R2D2M2 and per-margin ICC: sd(...) ~ r2d2(...)
Use the same marker on an sd address when the decomposition is over the marginal variances of a shared random-effect block rather than over one linear predictor's population coefficients:
Joint random-effect R2D2M2joint_ranef_r2d2 = (@brm begin
sigma_pk ~ Exponential(1)
sigma_qt ~ Exponential(1)
log_Vc ~ 1 + (1 | p | subject)
log_k10 ~ 1 + (1 | p | subject)
qt_base ~ 1 + (1 | p | subject)
qt_slope ~ 1 + (1 | p | subject)
sd(:, p) ~ r2d2(mean_R2=0.5, prec_R2=2,
concentration=1, reference_scale=sigma_pk)
sd(qt_base, p) ~ r2d2(reference_scale=sigma_qt)
sd(qt_slope, p) ~ r2d2(reference_scale=sigma_qt)
cor(:, p) ~ LKJCholesky(4, 2)
end)((; subject=[1, 1, 2, 2]))BRMI:
sigma_pk ~ Exponential(1)
sigma_qt ~ Exponential(1)
subject: data (eltype=Int64, n=4)
log_Vc ~ 1 + (1 | p | subject)
log_k10 ~ 1 + (1 | p | subject)
qt_base ~ 1 + (1 | p | subject)
qt_slope ~ 1 + (1 | p | subject)
effect(sd, p) ~ r2d2(; mean_R2=0.5, prec_R2=2, concentration=1, reference_scale=sigma_pk)
effect(sd, p, qt_base) ~ r2d2(; reference_scale=sigma_qt)
effect(sd, p, qt_slope) ~ r2d2(; reference_scale=sigma_qt)
effect(cor, p) ~ LKJCholesky(4, 2)StanBlocks unsupported for this BRM example
sbimpl: this `@brm` declares no observation — no `response ~ distribution(...)` statement binds data, and none is present without data either. Every `@brm` needs an observation statement; for prior draws keep the statement and omit the response column from the data — dropping the statement is not supported.Stan emission unavailable because StanBlocks lowering failed.
sbimpl: this `@brm` declares no observation — no `response ~ distribution(...)` statement binds data, and none is present without data either. Every `@brm` needs an observation statement; for prior draws keep the statement and omit the response column from the data — dropping the statement is not supported.Turing unsupported for this BRM example
Turing backend: direct execution requires at least one observed likelihoodThe block-wide statement samples one global R² and one Dirichlet simplex over the four marginal variances. Margin j is reconstructed as
tau[j] = reference_scale[j] * sqrt(phi[j] * R2 / (1 - R2))so the simplex allocates explained-variance odds in observation-scale units. The two more-specific statements change only the QT margins' reference scale; they do not create extra R² parameters. This is the explicit multi-response form: BRM does not guess which observation scale belongs to a latent predictor. The LKJ factor and the non-centred b = diag_pre_multiply(tau, L) * z construction are unchanged.
R2=Beta(a,b) and mean_R2=/prec_R2= are equivalent spellings. alpha= and concentration= are aliases for the symmetric Dirichlet concentration. reference_scale= is required and may be a positive numeric formula constant or an earlier sampled scalar parameter.
For independent per-margin ICC priors, omit the block-wide statement and address margins separately:
Independent per-margin ICC priorsicc_ranef_r2d2 = (@brm begin
sigma_pk ~ Exponential(1)
sigma_qt ~ Exponential(1)
log_Vc ~ 1 + (1 | p | subject)
qt_base ~ 1 + (1 | p | subject)
sd(log_Vc, p) ~ r2d2(R2=Beta(1, 1), reference_scale=sigma_pk)
sd(qt_base, p) ~ r2d2(mean_R2=0.5, prec_R2=2,
reference_scale=sigma_qt)
cor(:, p) ~ LKJCholesky(2, 2)
end)((; subject=[1, 1, 2, 2]))BRMI:
sigma_pk ~ Exponential(1)
sigma_qt ~ Exponential(1)
subject: data (eltype=Int64, n=4)
log_Vc ~ 1 + (1 | p | subject)
qt_base ~ 1 + (1 | p | subject)
effect(sd, p, log_Vc) ~ r2d2(; R2=Beta(1, 1), reference_scale=sigma_pk)
effect(sd, p, qt_base) ~ r2d2(; mean_R2=0.5, prec_R2=2, reference_scale=sigma_qt)
effect(cor, p) ~ LKJCholesky(2, 2)StanBlocks unsupported for this BRM example
sbimpl: this `@brm` declares no observation — no `response ~ distribution(...)` statement binds data, and none is present without data either. Every `@brm` needs an observation statement; for prior draws keep the statement and omit the response column from the data — dropping the statement is not supported.Stan emission unavailable because StanBlocks lowering failed.
sbimpl: this `@brm` declares no observation — no `response ~ distribution(...)` statement binds data, and none is present without data either. Every `@brm` needs an observation statement; for prior draws keep the statement and omit the response column from the data — dropping the statement is not supported.Turing unsupported for this BRM example
Turing backend: direct execution requires at least one observed likelihoodEach address then owns its own R². A one-margin group emits a deterministic simplex[1], so R2 = tau² / (tau² + reference_scale²) is exactly the ICC. Margins not addressed by either statement keep their ordinary half-standard-Normal scale prior. A block may not mix R2D2 addresses with explicit direct-scale sd(...) priors; use separate blocks when both constructions are required.
Random-effect R2D2 is SBBRMI-only. It supports ordinary non-centred shared-ID blocks, cor(:, ID), and reprocess(...; resample_groups=[group]). Centered, stratified, and multi-membership blocks fail loudly.
Joint R2D2M2 budget over coefficients, contrasts and random effects: include=
When several linear predictors share one correlated block AND one covariate right-hand side, brms' R2D2M2 puts ONE global R² and ONE Dirichlet over the union of the population coefficients (continuous and categorical) and the random-effect variances. Spell that by adding include= to the block-wide statement:
Joint R2D2M2 budget with include=joint_budget_r2d2 = (@brm begin
sigma_pk ~ Exponential(1)
sigma_qt ~ Exponential(1)
log_Vc ~ 1 + wt + indication + (1 | p | subject)
log_k10 ~ 1 + wt + indication + (1 | p | subject)
qt_base ~ 1 + wt + indication + (1 | p | subject)
sd(:, p) ~ r2d2(mean_R2=0.5, prec_R2=2, concentration=1,
reference_scale=sigma_pk,
include=(:population, :contrasts))
sd(qt_base, p) ~ r2d2(reference_scale=sigma_qt)
cor(:, p) ~ LKJCholesky(3, 2)
end)((;
subject=[1, 1, 2, 2, 3, 3],
wt=[-1.0, 0.5, 0.2, -0.3, 1.1, -0.6],
indication=[1, 2, 1, 2, 2, 1],
))BRMI:
sigma_pk ~ Exponential(1)
sigma_qt ~ Exponential(1)
wt: data (eltype=Float64, n=6)
indication: data (eltype=Int64, n=6)
subject: data (eltype=Int64, n=6)
log_Vc ~ 1 + wt + indication + (1 | p | subject)
log_k10 ~ 1 + wt + indication + (1 | p | subject)
qt_base ~ 1 + wt + indication + (1 | p | subject)
effect(sd, p) ~ r2d2(; mean_R2=0.5, prec_R2=2, concentration=1, reference_scale=sigma_pk, include=(:population, :contrasts))
effect(sd, p, qt_base) ~ r2d2(; reference_scale=sigma_qt)
effect(cor, p) ~ LKJCholesky(3, 2)StanBlocks unsupported for this BRM example
sbimpl: this `@brm` declares no observation — no `response ~ distribution(...)` statement binds data, and none is present without data either. Every `@brm` needs an observation statement; for prior draws keep the statement and omit the response column from the data — dropping the statement is not supported.Stan emission unavailable because StanBlocks lowering failed.
sbimpl: this `@brm` declares no observation — no `response ~ distribution(...)` statement binds data, and none is present without data either. Every `@brm` needs an observation statement; for prior draws keep the statement and omit the response column from the data — dropping the statement is not supported.Turing unsupported for this BRM example
Turing backend: direct execution requires at least one observed likelihoodinclude= names which population components of every predictor slicing |p| join the block's single R²/Dirichlet: :population is the non-intercept continuous beta_pop columns, :contrasts the categorical treatment-contrast coefficients (cat_* blocks), and :ranef the margins, which are always allocated and may be listed for readability. The simplex is ordered margins first, then each scoped predictor's population columns and contrast blocks in formula order. A component of predictor m is measured in m's own margin reference — its Intercept margin, or its single margin — so a margin keeps
tau[j] = reference_scale[m] * sqrt(phi[j] * R2 / (1 - R2))and a coefficient or contrast takes
beta_scale[k] = reference_scale[m] * sqrt(phi[k] * R2 / ((1 - R2) * Var(x_k)))which keeps the whole-predictor form's design-column variance adjustment (a contrast's dummy column has variance p * (1 - p) for level frequency p). Intercepts stay outside and keep their ordinary or explicitly overridden prior. Per-margin reference_scale overrides work exactly as above; the example gives qt_base's margin, coefficients, and contrast the QT residual scale.
With include=, reference_scale= becomes optional. An omitted margin reference is a sampled half-standard-normal parameter, so the statement then allocates latent between-subject variation: nothing observed anchors the unit, R2 is not an outcome R², and the sampled reference and R2 are identified by the data only through their product reference_scale² * R2 / (1 - R2). Prefer explicit references whenever a margin's scale is known.
Under a joint budget a per-column effect(lp, coef) ~ Normal(...) or effect(lp, categorical) ~ Normal(...) statement inside the scope is refused: it would silently pull that coefficient out of the simplex. A scoped predictor may not also carry effect(lp, :) ~ r2d2(...), and include= is accepted on the block-wide statement only — not on a per-margin override, not on the ICC form. ranef_effect_priors reports the joint statement with its include keyword. The emitted carriers keep their names (pop_<lp>_beta_pop, cat_<lp>_<col>_beta, the block's derived tau), so brm_population_effect_coordinates, brm_ranef_sd_coordinates, and reprocess(...; resample_groups=[group]) are unchanged.
Bounded scalar parameter priors
A non-data scalar prior may add finite numeric lower and/or upper declaration bounds. This is the direct spelling for a fitted positive scale with the same fixed-hyperparameter Normal kernel as a Stan real<lower=0> parameter:
Bounded scalar parameter priorbounded_scalar = (@brm begin
sigma ~ Normal(0, 2; lower=0.0)
mu ~ 1 + x
y ~ Normal(mu, sigma)
end)((;
x=[-1.0, 0.5, 2.0],
y=[0.2, 1.1, -0.4],
))BRMI:
sigma ~ Normal(0, 2; lower=0.0)
x: data (eltype=Float64, n=3)
mu ~ 1 + x
y ~ Normal(mu, sigma)SBBRMI with data keys = [:x, :y]
emitted @slic body:
begin
sigma ~ normal(0, 2; lower = 0.0)
X_mu = hcat(rep_vector(1.0, num_elements(x)), x)
pop_mu ~ _popefs_coefs(; X = X_mu)
y ~ normal_id_glm(X_mu, 0.0, pop_mu, sigma)
mu = X_mu * pop_mu
endfunctions {
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,
real 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,
real 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,
real 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 x_n;
vector[x_n] x;
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;
}
transformed parameters {
vector[pop_mu_n_covariates] pop_mu = pop_mu_beta_pop;
}
model {
sigma ~ normal(0, 2);
pop_mu_beta_pop ~ std_normal();
y ~ normal_id_glm(X_mu, 0.0, pop_mu, sigma);
}
generated quantities {
vector[x_n] y_likelihood = normal_id_glm_lpdfs(y, X_mu, 0.0, pop_mu, sigma);
vector[x_n] y_gen = normal_id_glm_vector_rng(y_n, X_mu, 0.0, pop_mu, sigma);
vector[x_n] mu = (X_mu * pop_mu);
}#= line 0 =# Turing.@model(function brm_model(y, X_mu)
sigma ~ BayesianRegressionModelsTuringExt._brm_constrained_kernel(Distributions.Normal(0, 2); lower = 0.0)
beta_pop ~ Distributions.product_distribution(Base.fill(Distributions.Normal(), 2))
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)Only lower and upper are accepted, both bounds and all distribution arguments must be numeric formula constants, and two bounds must satisfy lower < upper. BRM emits the bound on the Stan declaration and the ordinary family kernel in the model block, exactly as in hand-written Stan. This surface is SBBRMI-only; a bound with sampled hyperparameters needs an explicit normalized parameterization and is rejected.
Scalar horseshoe prior: coef ~ Horseshoe(...)
Horseshoe attaches a scalar Carvalho–Polson–Scott shrinkage hierarchy to an explicitly named coefficient:
Scalar horseshoe priorhorseshoe_model = (@brm begin
beta_sparse ~ Horseshoe(local_scale=0.5, global_scale=0.1)
sigma ~ Exponential(1)
mu = beta_sparse * x
y ~ Normal(mu, sigma)
end)((; x=[-1.0, 0.5, 2.0], y=[0.2, 1.1, -0.4]))BRMI:
beta_sparse ~ Horseshoe(; local_scale=0.5, global_scale=0.1)
sigma ~ Exponential(1)
x: data (eltype=Float64, n=3)
:mu = beta_sparse * x
y ~ Normal(mu, sigma)SBBRMI with data keys = [:x, :y]
emitted @slic body:
begin
beta_sparse ~ _sb_horseshoe_scaled(; local_scale = 0.5, global_scale = 0.1)
sigma ~ exponential(1.0 ./ 1)
mu = beta_sparse .* x
y ~ normal(mu, sigma)
endfunctions {
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 y_n;
vector[y_n] y;
}
transformed data {
}
parameters {
real beta_sparse_raw;
real<lower=0.0> beta_sparse_lambda;
real<lower=0.0> beta_sparse_tau;
real<lower=0.0> sigma;
}
transformed parameters {
real beta_sparse = (beta_sparse_raw * beta_sparse_lambda * beta_sparse_tau);
vector[x_n] mu = (beta_sparse .* x);
}
model {
beta_sparse_raw ~ std_normal();
beta_sparse_lambda ~ cauchy(0.0, 0.5);
beta_sparse_tau ~ cauchy(0.0, 0.1);
sigma ~ exponential((1.0 ./ 1));
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);
}#= line 0 =# Turing.@model(function brm_model(y, x)
beta_sparse ~ DynamicPPL.to_submodel(BayesianRegressionModels._brm_turing_horseshoe_prior(; local_scale = 0.5, global_scale = 0.1))
sigma ~ Distributions.Exponential(1)
mu = [beta_sparse * x[i] for i = Base.eachindex(y)]
begin
for i = Base.eachindex(y)
y[i] ~ Distributions.Normal(mu[i], sigma)
end
end
(; sigma = sigma, beta_sparse = beta_sparse, mu = mu, response = y)
end)The emitted non-centered hierarchy is
raw ~ Normal(0, 1)
lambda ~ half-Cauchy(0, local_scale)
tau ~ half-Cauchy(0, global_scale)
beta_sparse = raw * lambda * tauA sampled coefficient enters the linear predictor through an assignment (mu = beta_sparse * x), not a ~ formula term. A ~ right-hand side is formula-style only — an intercept plus population/group summands whose coefficients SBBRMI allocates — so mu ~ 0 + beta_sparse * x is rejected: beta_sparse is your own sampled scalar, not a data column, and the .* product is what you want to write directly. Use = for any linear predictor that multiplies a sampled parameter by data.
local_scale and global_scale are optional positive finite formula constants. Both default to 1.0; the no-keyword spelling Horseshoe() retains the historical emission byte for byte. Literal arithmetic such as local_scale=1/2 is evaluated with Julia semantics. Positional arguments, unknown keywords, nonnumeric values, zero, negative, and non-finite scales are rejected during SBBRMI construction.
The current surface is scalar and each coef ~ Horseshoe(...) call owns its own (raw, lambda, tau) triple. Consequently tau is “global” only within that scalar hierarchy; it is not shared across several coefficients. A genuinely shared global scale requires a vector/group horseshoe declaration, which this marker does not imply. The standardized raw draw is intentionally not configurable because its scale would duplicate lambda/tau.
Like other sbimpl prior surfaces, this is implemented by SBBRMI; do not assume VBRMI has the same marker-specific lowering.
Simplex-valued parameter: s ~ Dirichlet(...)
A non-data left-hand side with a Dirichlet right-hand side declares a simplex-valued parameter — Stan's simplex[K] — rather than a linear predictor. Like every other scalar-parameter prior declaration it is addressable by name anywhere later in the formula block, kernel(...) cells included. It is implemented by both SBBRMI and TuringBRMI.
The multi-axis population PK kernel shows the executable base shape: ordinary formula parameters feed a kernel(...) cell, while ragged(...) attaches a secondary observation frame to the subject axis.
The inciting shape is a per-event multiplier: diet scales dose bioavailability inside the cell, so what the model needs is a K-simplex parameter indexed by an ordinal level — not a per-observation design column.
What it emits
For diet_share ~ Dirichlet(3, 1.0), the generated Stan declares the concentration vector as data, the parameter as simplex[3], and the model statement as diet_share ~ dirichlet(diet_share_alpha). Complete emissions are shown only through build-generated comparisons rather than copied Stan fences.
The concentration is registered as data under <name>_alpha, which is what sizes the simplex; the name is reserved, so a collision is rejected rather than overwritten. A simplex[K] costs K - 1 unconstrained coordinates and reports K constrained ones.
Because the concentration is data, Stan drops the Dirichlet log-normalizer — a function of alpha alone — from target, exactly as it does for any other ~ statement with data-only hyperparameters. The sampled distribution is unchanged.
Accepted spellings
Only the two genuine Distributions.Dirichlet constructors, with Julia's parameterization preserved:
| spelling | meaning |
|---|---|
Dirichlet(alpha) | concentration vector literal, e.g. Dirichlet([2.0, 1.5, 3.0]) |
Dirichlet(K, a) | symmetric: K components, each concentration a |
@brm is a macro over the formula block, so a bare Julia symbol on the right-hand side is parsed as a formula local, not interpolated — Dirichlet(3, alpha) with a captured alpha reaches the backend as a column carrier and is rejected by name. Spell the concentration as a literal; there is no $ escape ($ outside a quote is a Julia syntax error, so the macro never sees it). This matches the rest of the formula surface — r2d2's alpha= and effect(:, x) ~ Normal(0, 0.25) are literals for the same reason.
Concentrations must be finite and strictly positive, and K >= 2: a one-element simplex is deterministically [1.0], so there is no parameter to sample.
Difference from mo(c) / mo1(c)
mo and mo1 are population linear-predictor terms. Their Dirichlet increments live inside a submodel and are not addressable from the formula; what a bare share ~ mo1(c) hands back is a per-row monotonic contrast vector of length nrow(df). Indexing that by an ordinal level is a silent double indirection — it transpiles and stanc-checks, but computes something else. Use Dirichlet when you want the simplex itself.
Current limits
The left-hand side must be a non-data name. A data-backed response with a
Dirichletright-hand side is a simplex-valued likelihood, which has no density/pointwise/predictive support and is rejected by the family path.The concentration is a constant, not a hyperprior:
Dirichlet(alpha)withalphaa sampled parameter is not supported.Dirichlet(K)— no concentration — is not aDistributions.jlconstructor and is rejected; writeDirichlet(K, 1.0)for the flat case.Not available to
VBRMI.s,t2,r2d2, and this simplex surface are supported by bothSBBRMIandTuringBRMI.
Vector-valued parameter: x ~ MvNormal(...)
A non-data left-hand side with an MvNormal right-hand side declares a vector-valued parameter — Stan's vector[n] — rather than a linear predictor. It is addressable by name anywhere later in the formula block: as an argument of a @deffun in a top-level assignment, inside a custom family call, or in a kernel(...) cell. SBBRMI emits Stan's vector-parameter form (decision 187g4va), while TuringBRMI lowers the retained generic callable directly.
The inciting shape is a latent path whose innovations the formula wants to state directly — the log-reproduction-number random walk of a renewal model written with explicit innovations, eps ~ MvNormal(zeros(T - 1), 1.0), or a vector of per-group seeds, log_I0 ~ MvNormal(seed_mean, 0.25 * I). (The Epidemic renewal models page states the same walk with the rw(time) term and the seeds as cell means.) Before this seam the only way to get such a vector was a one-cell kernel(...) with a dummy grouping random effect.
Random walk on a top-level vector parameterusing StanBlocks
StanBlocks.@deffun begin
rw_path(init::real, sig::real, eps::vector[K])::vector[K + 1] =
append_row(init, init + sig * cumulative_sum(eps))
end
random_walk = (@brm begin
sig ~ Normal(0.0, 0.05; lower=0.0)
init ~ Normal(0.0, 1.0)
eps ~ MvNormal(zeros(length(time) - 1), 1.0) # vector[T-1] of iid N(0, 1)
walk = rw_path(init, sig, eps)
y ~ Normal(walk, 0.3)
end)((;
time=collect(1.0:8),
y=[0.4, 1.1, 0.9, 1.6, 1.2, 2.0, 1.8, 2.5],
))BRMI:
sig ~ Normal(0.0, 0.05; lower=0.0)
init ~ Normal(0.0, 1.0)
time: data (eltype=Float64, n=8)
eps ~ MvNormal(zeros(-(length(time), 1)), 1.0)
:walk = rw_path(init, sig, eps)
y ~ Normal(walk, 0.3)SBBRMI with data keys = [:eps_mu, :eps_n, :time, :y]
emitted @slic body:
begin
sig ~ normal(0.0, 0.05; lower = 0.0)
init ~ normal(0.0, 1.0)
eps::vector[eps_n] ~ normal(eps_mu, 1.0)
walk = (Main.__atexample__1152.rw_path)(init, sig, eps)
y ~ normal(walk, 0.3)
endfunctions {
vector rw_path(
real init,
real sig,
vector eps
) {
int K = dims(eps)[1];
return append_row(init, (init + (sig * cumulative_sum(eps))));
}
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 eps_n;
int eps_mu_n;
vector[eps_mu_n] eps_mu;
int y_n;
vector[y_n] y;
}
transformed data {
}
parameters {
real<lower=0.0> sig;
real init;
vector[eps_n] eps;
}
transformed parameters {
vector[(eps_n + 1)] walk = rw_path(init, sig, eps);
}
model {
sig ~ normal(0.0, 0.05);
init ~ normal(0.0, 1.0);
eps ~ normal(eps_mu, 1.0);
y ~ normal(walk, 0.3);
}
generated quantities {
vector[y_n] y_likelihood = normal_lpdfs(y, walk, 0.3);
vector[y_n] y_gen = normal_vector_rng(y_n, walk, 0.3);
}#= line 0 =# Turing.@model(function brm_model(y, time)
sig ~ BayesianRegressionModelsTuringExt._brm_constrained_kernel(Distributions.Normal(0.0, 0.05); lower = 0.0)
init ~ Distributions.Normal(0.0, 1.0)
eps ~ Distributions.MvNormal(Base.zeros(Base.length(time) - 1), 1.0)
walk = Main.__atexample__1152.rw_path(init, sig, eps)
begin
for i = Base.eachindex(y)
y[i] ~ Distributions.Normal(walk, 0.3)
end
end
(; sig = sig, init = init, eps = eps, walk = walk, response = y)
end)What it emits
For eps ~ MvNormal(zeros(length(time) - 1), 1.0) the generated Stan declares the dimension as data (<name>_n), the mean as data (<name>_mu), the parameter as vector[eps_n] eps, and the model statement as eps ~ normal(eps_mu, 1.0) — Stan's vectorised univariate normal, no matrix. A full covariance emits multi_normal(<name>_mu, <name>_scale). The three data names are reserved, so a collision is rejected rather than overwritten; a data column used as the mean is referenced by its own name instead of being copied.
A vector nothing later reads is a generated quantity, not a sampler parameter: StanBlocks' activity analysis lowers an unused declaration, exactly as it does for an unused random effect. A top-level assignment's right-hand side is Julia-evaluated by the macro, so Stan builtins (append_row, cumulative_sum) reach it through a @deffun, as in the example.
Accepted spellings
Only the genuine Distributions.MvNormal constructors, with Julia's parameterization preserved (the σ forms are deprecated in Distributions.jl but keep their Distributions meaning — a standard deviation):
| spelling | meaning | emitted |
|---|---|---|
MvNormal(mu, Σ::Matrix) | mean and covariance | multi_normal |
MvNormal(mu, Diagonal(v)) | variances on the diagonal | vectorised normal with sqrt.(v) |
MvNormal(mu, λ * I) | covariance λ·I | vectorised normal with sqrt(λ) |
MvNormal(Σ) | zero mean | multi_normal |
MvNormal(mu, σ::Real) | isotropic, standard deviation σ | vectorised normal |
MvNormal(mu, σ::Vector) | standard deviations | vectorised normal |
MvNormal(n, σ::Real) | zero mean, n components | vectorised normal |
The mean is a numeric vector literal, a data column, or a data-only expression over data columns (zeros(length(time) - 1), fill(c, k)) — evaluated in Julia at build time and shipped as data. The scale may be a literal, a data column, or a sampled scalar parameter (MvNormal(zeros(k), sig) emits normal(eps_mu, sig)). The dimension comes from the mean, from a vector/matrix scale, or from the integer form; when the mean itself involves a parameter and the scale is a scalar, there is no data-determinable size and the statement is rejected with a message saying so.
Current limits
The left-hand side must be a non-data name. A vector-valued response keeps its own spelling (
MvNormalCholesky,[y1, y2] ~ ...);MvNormalis not added to the family table.No keywords: bounds on a multivariate normal parameter are not supported.
A parameter-bearing covariance (an
LKJCovarianceFactorproduct, a sampled matrix) is not admitted here; use the scalar-scale form or compose the vector inside a@deffun.Not available to
VBRMI.TuringBRMIretains the generic callable, so the same statement lowers there through its own path;Dirichlet,s, andr2d2are likewise supported by both executable backends.
Function-valued arguments: lambdas and do blocks
A top-level assignment may hand a function to a higher-order @deffun: written inline as the first argument, or as a trailing do block, which is the same thing — @brm folds the block in as the first positional argument, as kernel(...) do … end already does. Higher-order functions therefore take their function argument first.
The body is Stan-side code, resolved by StanBlocks as a closure. It may read
its own parameters (
lbelow),sampled model parameters (
rate) — they are passed into the generated Stan function as arguments, so they are estimated, not frozen in, andshared data vectors by name: a free name that is a data field and not a formula name is registered as Stan data, the way a
kernel(...)cell body captures one.
The example is a distributed lag whose weights decay at a sampled rate — a filter that a data vector of weights could not express, because the weights depend on a parameter.
A do block that reads a sampled parameterusing StanBlocks
StanBlocks.@deffun begin
# Y[t] = sum over lags l = 0 .. L-1 of f(l) * x[t - l]
lagged(f, x::vector[T], L::int)::vector[T] = begin
Y::vector[T]
for t in 1:T
acc = 0.0
for l in 0:(L - 1)
if t - l >= 1
acc += f(l) * x[t - l]
end
end
Y[t] = acc
end
Y
end
end
decaying_lag = (@brm begin
rate ~ Normal(1.0, 0.5; lower=0.0)
log_mu ~ 1 + rw(time)
Y = lagged(exp(log_mu), 3) do l
exp(-rate * l) # reads the sampled `rate`
end
y ~ Poisson(Y)
end)((;
time=collect(1.0:12),
y=[3, 4, 6, 5, 8, 9, 12, 11, 15, 18, 17, 21],
))BRMI:
rate ~ Normal(1.0, 0.5; lower=0.0)
time: data (eltype=Float64, n=12)
log_mu ~ 1 + rw(time)
:Y = lagged((l,)->begin
#= brm-docs-example.jl:22 =#
exp(-rate * l)
end, exp(log_mu), 3)
y ~ Poisson(Y)SBBRMI with data keys = [:rw_log_mu_time_n_steps, :rw_log_mu_time_time_idx, :time, :y]
emitted @slic body:
begin
rate ~ normal(1.0, 0.5; lower = 0.0)
X_log_mu = hcat(rep_vector(1.0, num_elements(time)))
pop_log_mu ~ popefs(; X = X_log_mu)
rw_log_mu_time ~ _sb_rw1(; n_steps = rw_log_mu_time_n_steps, time_idx = rw_log_mu_time_time_idx)
log_mu = pop_log_mu + rw_log_mu_time
Y = (Main.__atexample__1184.lagged)(((l,)->begin
#= brm-docs-example.jl:22 =#
exp(-rate * l)
end), (exp)(log_mu), 3)
y ~ poisson(Y)
endfunctions {
matrix hcat(vector x) {
int n = dims(x)[1];
return to_matrix(x, n, 1);
}
vector random_walk_rows(
real sigma,
vector z,
array[] int idx
) {
int N = dims(idx)[1];
vector[(dims(z)[1] + 1)] x = random_walk_path(sigma, z);
vector[N] out;
for(i in 1:N) {
out[i] = x[idx[i]];
}
return out;
}
vector random_walk_path(
real sigma,
vector z
) {
int n = dims(z)[1];
return append_row(0.0, (sigma * cumulative_sum(z)));
}
vector lagged_closure_1(
real rate,
vector x,
int L
) {
int T = dims(x)[1];
vector[T] Y;
for(t in 1:T) {
real acc = 0.0;
for(l in 0:(L - 1)) {
if(((t - l) >= 1)) {
acc += (exp(((-rate) * l)) * x[(t - l)]);
}
}
Y[t] = acc;
}
return Y;
}
vector poisson_lpmfs(
array[] int obs,
vector lambda
) {
return jbroadcasted_poisson_lpmfs(obs, lambda);
}
vector jbroadcasted_poisson_lpmfs(
array[] int x1,
vector x2
) {
int n = dims(x1)[1];
vector[n] rv;
for(i in 1:n) {
rv[i] = poisson_lpmfs(broadcasted_getindex(x1, i), broadcasted_getindex(x2, i));
}
return rv;
}
real poisson_lpmfs(int args1, real args2) {
return poisson_lpmf(args1 | args2);
}
int broadcasted_getindex(array[] int x, int i) {
return x[i];
}
real broadcasted_getindex(vector x, int i) {
return x[i];
}
array[] int poisson_int_rng(
int anontok__1,
vector p
) {
int n = anontok__1;
return poisson_rng(p);
}
}
data {
int time_n;
vector[time_n] time;
int rw_log_mu_time_n_steps;
int rw_log_mu_time_time_idx_n;
array[rw_log_mu_time_time_idx_n] int rw_log_mu_time_time_idx;
int y_n;
array[y_n] int y;
}
transformed data {
matrix[num_elements(time), 1] X_log_mu = hcat(rep_vector(1.0, num_elements(time)));
int pop_log_mu_n_covariates = 1;
}
parameters {
real<lower=0.0> rate;
vector[pop_log_mu_n_covariates] pop_log_mu_beta_pop;
real<lower=0.0> rw_log_mu_time_sigma;
vector[(rw_log_mu_time_n_steps - 1)] rw_log_mu_time_z;
}
transformed parameters {
vector[num_elements(time)] pop_log_mu = (X_log_mu * pop_log_mu_beta_pop);
vector[rw_log_mu_time_time_idx_n] rw_log_mu_time = random_walk_rows(rw_log_mu_time_sigma, rw_log_mu_time_z, rw_log_mu_time_time_idx);
vector[num_elements(time)] log_mu = (pop_log_mu + rw_log_mu_time);
vector[num_elements(time)] Y = lagged_closure_1(rate, exp(log_mu), 3);
}
model {
rate ~ normal(1.0, 0.5);
pop_log_mu_beta_pop ~ std_normal();
rw_log_mu_time_sigma ~ normal(0.0, 0.2);
rw_log_mu_time_z ~ std_normal();
y ~ poisson(Y);
}
generated quantities {
vector[y_n] y_likelihood = poisson_lpmfs(y, Y);
array[y_n] int y_gen = poisson_int_rng(y_n, Y);
}Turing unsupported for this BRM example
Turing backend: assignment `Y` calls `lagged`, which defines no Julia methods (a Stan `@deffun` function has no Julia implementation). Row-dependent assignments lower to row-wise Julia comprehensions, so this call cannot execute. Express the computation with row-wise Julia code or fit this model with the Stan backend.In the generated Stan the closure is a function of its own (lagged_closure_1) and the call site passes what the body captured: Y = lagged_closure_1(rate, exp(log_mu), 3). The inline spelling lagged(l -> exp(-rate * l), exp(log_mu), 3) emits the same program.
The callee may live in another module — a package that ships @deffun operators — as long as the calling module imports it and is passed as mod=. This is a StanBlocks feature: TuringBRMI refuses the assignment by name, as it does for every Stan-only callee.