Air pollution: regional effects and automatic totals

Result

This study applies BRM's built-in total-coefficient construction to two supported versions of the AIR model from the Stan S2Z discussion. It retains the population effects and their priors, and compares ordinary brms, brms sum-to-zero (S2Z), and BRM totals on the same scientific quantities within each model.

For regional intercepts, S2Z CP + WHMC gave 279× the native-NCP baseline's total-gradient efficiency; fixed BRM totals gave 166× and both online total variants 128×. With independent intercepts and slopes, fixed total CP gave 179×, online gradient 140× and S2Z CP + WHMC 95.5×. Thus fixed total CP was about 1.87× as efficient as S2Z CP + WHMC in the latter pilot, while S2Z CP was about 1.68× as efficient as total CP in the intercept-only pilot.

These are one-chain exploratory comparisons. They show that centering is important after marginalization and that pilot cost can dominate post-hoc adaptation. They do not establish one uniformly best parameterization or adaptation loss.

Data and model variants

There are 6,003 observations of log ground-level particulate concentration (log_pm25) and a log satellite predictor (log_sat). We preserve row order and use cluster_region, with six group sizes 27, 3,534, 1,696, 418, 315 and 13. The uneven information per region makes a common centering choice potentially restrictive.

Both variants retain the population intercept and slope:

# Regional intercepts
bf(log_pm25 ~ log_sat + (1 | region))

# Independent regional intercepts and slopes
bf(log_pm25 ~ log_sat + (1 + log_sat || region))

The second formula deliberately removes the correlation between random intercepts and slopes. The correlated source model is outside the current automatic-total planner's supported scope. Removing a covariance parameter changes the posterior; every comparison row here uses the same explicitly simplified model. No population effects are deleted.

Let x denote log_sat. The observation model is

$$y_n\sim N({\mu }_n,{\sigma }^2),{\mu }_n={\beta }_0+{\beta }_1(x_n-\overline{x})+a_{j[n]}+b_{j[n]}{x}_n.$$

The intercept-only variant omits b[j]. The population intercept is defined at the observation-weighted mean predictor; random slopes use raw x. Priors match the generated brms models:

• Population intercept: Student-t(3, location 2.8, scale 2.5).

• Population slope: flat.

• Independent regional deviations: zero-mean Gaussian, with a separate SD for each included term.

• Group SDs and residual SD: half-Student-t(3, 0, 2.5).

The source captures also retain cluster_log_region and super_region specifications. The measurements on this page use cluster_region only.

BRM models

The authoring panes below read the exact fitted declarations. The backend views are regenerated during the docs build. Sampling uses StanBlocks/BridgeStan and WHMC; the generated Turing pane is for inspection only.

Regional intercepts

AIR with regional intercepts

using BayesianRegressionModels, Distributions, JSON, DelimitedFiles

function air_intercept_brm_model()
    reference = JSON.parsefile(joinpath(pkgdir(BayesianRegressionModels),
        "research", "air_total_effects", "reference", "cluster_region",
        "intercept_only", "ordinary_ncp.json"))
    data = (;log_pm25=Float64.(reference["Y"]),
             log_sat=Float64.(getindex.(reference["X"], 2)),
             region=Int.(reference["J_1"]))
    
    builder = @brm begin
        mu ~ 1 + center(log_sat) + (1 | regional_intercept | region)
        effect(mu, Intercept) ~ LocationScale(2.8,2.5,TDist(3))
        effect(mu, center_log_sat) ~ Flat()
        sd(:, regional_intercept) ~ LocationScale(0.,2.5,TDist(3))
        sigma ~ LocationScale(0.,2.5,TDist(3);lower=0.)
        log_pm25 ~ Normal(mu,sigma)
    end
    
    builder(data)
end

BRMI:
  log_sat: data (eltype=Float64, n=6003)
  region: data (eltype=Int64, n=6003)
  mu ~ 1 + center(log_sat) + (1 | regional_intercept | region)
  effect(mu, Intercept) ~ AffineDistribution(2.8, 2.5, TDist(3))
  effect(mu, center_log_sat) ~ Flat()
  effect(sd, regional_intercept) ~ AffineDistribution(0.0, 2.5, TDist(3))
  sigma ~ AffineDistribution(0.0, 2.5, TDist(3); lower=0.0)
  log_pm25 ~ Normal(mu, sigma)

SBBRMI with data keys = [:center_log_sat, :log_pm25, :log_sat, :region, :total_A_mu, :total_group_mu, :total_location_mu, :total_mixture_shape_mu, :total_ng_mu, :total_nk_mu, :total_nm_mu, :total_np_mu, :total_precision_mu]
configured submodels:
_brm_total_scales_configured_1 = Base.merge(BayesianRegressionModels._brm_total_scales, quote
            tau ~ (ValueFamily(brm_vector_prior_b28d90387aac6d90))(3.0, 0.0, 2.5; n = 1)
        end)
_popefs_generic_configured_1 = Base.merge(BayesianRegressionModels._popefs_generic, quote
            beta_pop::vector[n_covariates] ~ (ValueFamily(brm_vector_prior_9127091398de311e))(; n = 1)
        end)
emitted @slic body:
begin
    total_scale_mu ~ _brm_total_scales_configured_1(; n = total_nk_mu)
    total_mixture_mu::vector[total_nm_mu] ~ gamma(total_mixture_shape_mu, total_mixture_shape_mu; lower = 0.0)
    total_conditional_precision_mu = [total_precision_mu[1] * total_mixture_mu[1]]
    total_mu::matrix[total_ng_mu, total_nk_mu] ~ brm_total(total_scale_mu, total_A_mu, total_location_mu, total_conditional_precision_mu)
    population_mu = brm_total_recover_rng(total_mu, total_scale_mu, total_A_mu, total_location_mu, total_conditional_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(center_log_sat)
    pop_mu ~ _popefs_generic_configured_1(; X = X_mu)
    mu = rows_dot_product(total_mu[total_group_mu, :], total_Z_mu) + pop_mu
    sigma ~ student_t(3, 0.0, 2.5; lower = 0.0)
    log_pm25 ~ normal(mu, sigma)
end

functions {
// value UDF brm_vector_prior_b28d90387aac6d90_lpdf
real brm_vector_prior_b28d90387aac6d90_lpdf(
    vector x,
    real arg_1,
    real arg_2,
    real arg_3
) {
    if((x[1] < 0.0)) {
        return negative_infinity();
    }
    return student_t_lpdf(x[1] | arg_1, arg_2, arg_3);
}
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);
}
// value UDF brm_vector_prior_9127091398de311e_lpdf
real brm_vector_prior_9127091398de311e_lpdf(vector x) {
    return 0.0;
}
vector normal_lpdfs(
    vector obs,
    vector loc,
    real scale
) {
    return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
    vector x1,
    vector x2,
    real x3
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = normal_lpdfs(broadcasted_getindex(x1, i), broadcasted_getindex(x2, i), x3);
    }
    return rv;
}
real normal_lpdfs(
    real args1,
    real args2,
    real args3
) {
    return normal_lpdf(args1 | args2, args3);
}
real broadcasted_getindex(vector x, int i) {
    return x[i];
}
vector normal_vector_rng(
    int anontok__1,
    vector a,
    real b
) {
    int n = anontok__1;
    if((n == 0)) {
        vector[n] rv;
        return rv;
    } else {
        return to_vector(normal_rng(a, b));
    }
}
}
data {
    int total_nm_mu;
    int total_mixture_shape_mu_n;
    vector[total_mixture_shape_mu_n] total_mixture_shape_mu;
    int total_precision_mu_n;
    vector[total_precision_mu_n] total_precision_mu;
    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_group_mu_n;
    array[total_group_mu_n] int total_group_mu;
    int center_log_sat_n;
    vector[center_log_sat_n] center_log_sat;
    int log_pm25_n;
    vector[log_pm25_n] log_pm25;
}
transformed data {
    matrix[num_elements(total_group_mu), 1] total_Z_mu = hcat(rep_vector(1.0, num_elements(total_group_mu)));
    matrix[center_log_sat_n, 1] X_mu = hcat(center_log_sat);
    int pop_mu_n_covariates = 1;
}
parameters {
    vector<lower=0.0>[1] total_scale_mu_tau;
    vector<lower=0.0>[total_nm_mu] total_mixture_mu;
    matrix[total_ng_mu, total_nk_mu] total_mu;
    vector[pop_mu_n_covariates] pop_mu_beta_pop;
    real<lower=0.0> sigma;
}
transformed parameters {
    vector<lower=0.0>[1] total_scale_mu = total_scale_mu_tau;
    vector[1] total_conditional_precision_mu = [(total_precision_mu[1] * total_mixture_mu[1])]';
    vector[center_log_sat_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_b28d90387aac6d90(3.0, 0.0, 2.5);
    total_mixture_mu ~ gamma(total_mixture_shape_mu, total_mixture_shape_mu);
    total_mu ~ brm_total(total_scale_mu, total_A_mu, total_location_mu, total_conditional_precision_mu);
    pop_mu_beta_pop ~ brm_vector_prior_9127091398de311e();
    sigma ~ student_t(3, 0.0, 2.5);
    log_pm25 ~ normal(mu, sigma);
}
generated quantities {
    vector[1] population_mu = brm_total_recover_rng(
        total_mu,
        total_scale_mu,
        total_A_mu,
        total_location_mu,
        total_conditional_precision_mu
    );
    matrix[total_ng_mu, total_A_mu_m] deviation_mu = brm_total_deviations(total_mu, (total_A_mu * population_mu));
    vector[log_pm25_n] log_pm25_likelihood = normal_lpdfs(log_pm25, mu, sigma);
    vector[log_pm25_n] log_pm25_gen = normal_vector_rng(log_pm25_n, mu, sigma);
}

#= line 0 =# Turing.@model(function brm_model(y, callable_2, X_mu, group_effects_mu_1, callable_3, callable_1)
        beta_pop ~ Distributions.product_distribution([callable_2(2.8, 2.5, Distributions.TDist(3)), BayesianRegressionModels.Flat()])
        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, (callable_3(0.0, 2.5, Distributions.TDist(3)),), nothing))
        group_effect_1 = group_effect_1 + group_1_1.effect
        eta_mu = eta_mu + group_effect_1
        mu = eta_mu
        sigma ~ BayesianRegressionModelsTuringExt._brm_constrained_kernel(callable_1(0.0, 2.5, Distributions.TDist(3)); lower = 0.0)
        begin
            for i = Base.eachindex(y)
                y[i] ~ Distributions.Normal(mu[i], sigma)
            end
        end
        (; mu = mu, sigma = sigma, response = y)
    end)

BRM integrates the population intercept and samples alpha[j] = beta0 + a[j]. The likelihood is alpha[j] + beta1*(x-xbar); the population slope remains explicit. For reporting, the regional intercept at raw x=0 is A[j] = alpha[j] - xbar*beta1.

Independent regional intercepts and slopes

AIR with independent regional intercepts and slopes

using BayesianRegressionModels, Distributions, JSON, DelimitedFiles

function air_independent_brm_model()
    reference = JSON.parsefile(joinpath(pkgdir(BayesianRegressionModels),
        "research", "air_total_effects", "reference", "cluster_region",
        "independent", "ordinary_ncp.json"))
    data = (;log_pm25=Float64.(reference["Y"]),
             log_sat=Float64.(getindex.(reference["X"], 2)),
             region=Int.(reference["J_1"]))
    
    builder = @brm begin
        mu ~ 1 + center(log_sat) + (1 | regional_intercept | region) + (0 + log_sat | regional_slope | region)
        effect(mu, Intercept) ~ LocationScale(2.8,2.5,TDist(3))
        effect(mu, center_log_sat) ~ Flat()
        sd(:, regional_intercept) ~ LocationScale(0.,2.5,TDist(3))
        sd(:, regional_slope) ~ LocationScale(0.,2.5,TDist(3))
        sigma ~ LocationScale(0.,2.5,TDist(3);lower=0.)
        log_pm25 ~ Normal(mu,sigma)
    end
    
    builder(data)
end

BRMI:
  log_sat: data (eltype=Float64, n=6003)
  region: data (eltype=Int64, n=6003)
  mu ~ 1 + center(log_sat) + (1 | regional_intercept | region) + ((0 + log_sat) | regional_slope | region)
  effect(mu, Intercept) ~ AffineDistribution(2.8, 2.5, TDist(3))
  effect(mu, center_log_sat) ~ Flat()
  effect(sd, regional_intercept) ~ AffineDistribution(0.0, 2.5, TDist(3))
  effect(sd, regional_slope) ~ AffineDistribution(0.0, 2.5, TDist(3))
  sigma ~ AffineDistribution(0.0, 2.5, TDist(3); lower=0.0)
  log_pm25 ~ Normal(mu, sigma)

SBBRMI with data keys = [:log_pm25, :log_sat, :region, :total_A_mu, :total_group_mu, :total_location_mu, :total_mixture_shape_mu, :total_ng_mu, :total_nk_mu, :total_nm_mu, :total_np_mu, :total_precision_mu]
configured submodels:
_brm_total_scales_configured_1 = Base.merge(BayesianRegressionModels._brm_total_scales, quote
            tau ~ (ValueFamily(brm_vector_prior_eecc99814bef47ac))(3.0, 0.0, 2.5, 3.0, 0.0, 2.5; n = 2)
        end)
emitted @slic body:
begin
    total_scale_mu ~ _brm_total_scales_configured_1(; n = total_nk_mu)
    total_mixture_mu::vector[total_nm_mu] ~ gamma(total_mixture_shape_mu, total_mixture_shape_mu; lower = 0.0)
    total_conditional_precision_mu = [total_precision_mu[1] * total_mixture_mu[1], total_precision_mu[2]]
    total_mu::matrix[total_ng_mu, total_nk_mu] ~ brm_total(total_scale_mu, total_A_mu, total_location_mu, total_conditional_precision_mu)
    population_mu = brm_total_recover_rng(total_mu, total_scale_mu, total_A_mu, total_location_mu, total_conditional_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)), log_sat)
    mu = rows_dot_product(total_mu[total_group_mu, :], total_Z_mu)
    sigma ~ student_t(3, 0.0, 2.5; lower = 0.0)
    log_pm25 ~ normal(mu, sigma)
end

functions {
// value UDF brm_vector_prior_eecc99814bef47ac_lpdf
real brm_vector_prior_eecc99814bef47ac_lpdf(
    vector x,
    real arg_1,
    real arg_2,
    real arg_3,
    real arg_4,
    real arg_5,
    real arg_6
) {
    if((x[1] < 0.0)) {
        return negative_infinity();
    }
    if((x[2] < 0.0)) {
        return negative_infinity();
    }
    return (student_t_lpdf(x[1] | arg_1, arg_2, arg_3) + student_t_lpdf(x[2] | arg_4, arg_5, arg_6));
}
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,
    vector y
) {
    int n = dims(x)[1];
    if (dims(y)[1] != n) reject("hcat: dim mismatch — `y` dim 1 (= ", dims(y)[1], ") does not match `n` (= ", n, "), inferred from `x` dim 1. `n` sizes: `x` dim 1 (= ", dims(x)[1], "), `y` dim 1 (= ", dims(y)[1], ").");
    return append_col(x, y);
}
vector normal_lpdfs(
    vector obs,
    vector loc,
    real scale
) {
    return jbroadcasted_normal_lpdfs(obs, loc, scale);
}
vector jbroadcasted_normal_lpdfs(
    vector x1,
    vector x2,
    real x3
) {
    int n = dims(x1)[1];
    vector[n] rv;
    for(i in 1:n) {
        rv[i] = normal_lpdfs(broadcasted_getindex(x1, i), broadcasted_getindex(x2, i), x3);
    }
    return rv;
}
real normal_lpdfs(
    real args1,
    real args2,
    real args3
) {
    return normal_lpdf(args1 | args2, args3);
}
real broadcasted_getindex(vector x, int i) {
    return x[i];
}
vector normal_vector_rng(
    int anontok__1,
    vector a,
    real b
) {
    int n = anontok__1;
    if((n == 0)) {
        vector[n] rv;
        return rv;
    } else {
        return to_vector(normal_rng(a, b));
    }
}
}
data {
    int total_nm_mu;
    int total_mixture_shape_mu_n;
    vector[total_mixture_shape_mu_n] total_mixture_shape_mu;
    int total_precision_mu_n;
    vector[total_precision_mu_n] total_precision_mu;
    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 log_sat_n;
    int total_group_mu_n;
    array[total_group_mu_n] int total_group_mu;
    vector[log_sat_n] log_sat;
    int log_pm25_n;
    vector[log_pm25_n] log_pm25;
}
transformed data {
    matrix[log_sat_n, 2] total_Z_mu = hcat(rep_vector(1.0, num_elements(total_group_mu)), log_sat);
}
parameters {
    vector<lower=0.0>[2] total_scale_mu_tau;
    vector<lower=0.0>[total_nm_mu] total_mixture_mu;
    matrix[total_ng_mu, total_nk_mu] total_mu;
    real<lower=0.0> sigma;
}
transformed parameters {
    vector<lower=0.0>[2] total_scale_mu = total_scale_mu_tau;
    vector[2] total_conditional_precision_mu = [(total_precision_mu[1] * total_mixture_mu[1]), total_precision_mu[2]]';
    vector[log_sat_n] mu = rows_dot_product(total_mu[total_group_mu, :], total_Z_mu);
}
model {
    total_scale_mu_tau ~ brm_vector_prior_eecc99814bef47ac(3.0, 0.0, 2.5, 3.0, 0.0, 2.5);
    total_mixture_mu ~ gamma(total_mixture_shape_mu, total_mixture_shape_mu);
    total_mu ~ brm_total(total_scale_mu, total_A_mu, total_location_mu, total_conditional_precision_mu);
    sigma ~ student_t(3, 0.0, 2.5);
    log_pm25 ~ normal(mu, sigma);
}
generated quantities {
    vector[2] population_mu = brm_total_recover_rng(
        total_mu,
        total_scale_mu,
        total_A_mu,
        total_location_mu,
        total_conditional_precision_mu
    );
    matrix[total_ng_mu, total_A_mu_m] deviation_mu = brm_total_deviations(total_mu, (total_A_mu * population_mu));
    vector[log_pm25_n] log_pm25_likelihood = normal_lpdfs(log_pm25, mu, sigma);
    vector[log_pm25_n] log_pm25_gen = normal_vector_rng(log_pm25_n, mu, sigma);
}

#= line 0 =# Turing.@model(function brm_model(y, callable_2, X_mu, group_effects_mu_1, callable_3, group_effects_mu_2, callable_4, callable_1)
        beta_pop ~ Distributions.product_distribution([callable_2(2.8, 2.5, Distributions.TDist(3)), BayesianRegressionModels.Flat()])
        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, (callable_3(0.0, 2.5, Distributions.TDist(3)),), nothing))
        group_effect_1 = group_effect_1 + group_1_1.effect
        group_1_2 ~ DynamicPPL.to_submodel(BRM.turing_group_effect(group_effects_mu_2, (callable_4(0.0, 2.5, Distributions.TDist(3)),), nothing))
        group_effect_1 = group_effect_1 + group_1_2.effect
        eta_mu = eta_mu + group_effect_1
        mu = eta_mu
        sigma ~ BayesianRegressionModelsTuringExt._brm_constrained_kernel(callable_1(0.0, 2.5, Distributions.TDist(3)); lower = 0.0)
        begin
            for i = Base.eachindex(y)
                y[i] ~ Distributions.Normal(mu[i], sigma)
            end
        end
        (; mu = mu, sigma = sigma, response = y)
    end)

BRM integrates both population coefficients and samples

$$A_j={\beta }_0-\overline{x}{\beta }_1+a_j,B_j={\beta }_1+b_j,{\mu }_n=A_{j[n]}+B_{j[n]}{x}_n.$$

In both variants, a Gamma(3/2, rate=3/2) precision multiplier represents the Student-t intercept exactly as a conditional Gaussian. The resulting Gaussian integral is evaluated analytically in O(J) for fixed coefficient dimension. Original population parameters remain recoverable from their exact conditional distribution.

The intercept-only target has 10 sampled dimensions in ordinary brms, totals and S2Z: the one integrated population coefficient is replaced by a mixture variable. The intercept-and-slope target has 17 ordinary dimensions and 16 total/S2Z dimensions.

Sampling and adaptation

using StanBlocks, BridgeStan, WarmupHMC, Enzyme, Random
using DifferentiationInterface: AutoEnzyme

sb = SBBRMI(air_independent_brm_model(); mod=@__MODULE__)
problem = StanBlocks.stan_instantiate(sb.model)
names = BridgeStan.param_unc_names(problem.model)
adaptive = adaptive_centering_problem(sb, problem, AutoEnzyme(); centeredness=0.0)
fit = adaptive_warmup_mcmc(Xoshiro(1), adaptive;
    n_draws=2000, monitor_ess=true, nonlinear_adapt=true)

draws = permutedims(fit.posterior_position)
recovered = recover_population_draws(sb, draws, names; rng=Xoshiro(404))

The centering family uses each regional total directly at CP and scales it about its prespecified reference location at NCP. The integrated prior remains coupled, so this NCP is not a whitening transformation. centeredness=1.0 selects CP; either endpoint stays fixed with nonlinear_adapt=false.

Each post-hoc arm uses select_total_centeredness on the same 2,000-draw total-NCP pilot, followed by a fresh fit. Both position/Jacobian and position–gradient losses use the built-in grid 0:0.1:1. Saved gradients are transported into the compiled model frame before gradient-based selection. Both losses are also tested online during WHMC warmup. The harness selects the existing position-loss weights locally; the default online loss uses position and gradient.

brms S2Z auto uses its branch's Pathfinder/Fisher precursor to choose fixed weights. Native Stan and WHMC receive the same resolved target and weights; the precursor is charged to both workflows. This is distinct from online adaptation.

Every arm retains 2,000 draws from one chain, seed 1. Native Stan uses 1,000 warmup iterations, target acceptance 0.8 and maximum tree depth 10. WHMC uses its adaptive warmup and Pathfinder initialization. All representations start at the same pooled physical coefficients, with SDs estimated from regional regressions and mixture precision one.

Gradient counts measure actual evaluations. The native harness adds a zero-contribution C++ counter to the model, leaving the NUTS sampler unchanged, checks retained increments against leapfrogs plus one, and records final counts for all sampling and precursor processes. WHMC counts density-and-gradient requests at its target wrapper. Gradient count is a work proxy; different targets need not take identical time per gradient.

One scientific scope for each model

Each model has its own ordinary brms NCP + native Stan baseline. Both efficiency columns are relative to that baseline: minimum bulk ESS divided by sampling gradients, and minimum bulk ESS divided by total gradients. Total costs include initialization, warmup and required pilots. The different model posteriors are not pooled into one minimum.

Regional intercepts: 10 quantities

The common quantities are the population intercept and slope, group SD, residual SD, and six regional total intercepts at raw predictor zero.

Method	Total gradients	Sampling efficiency	Total efficiency
brms NCP · Native Stan	1,060,905	1×	1×
brms NCP · WHMC	110,573	1.52×	1.94×
brms S2Z CP · Native Stan	117,637	21.4×	17.4×
brms S2Z CP · WHMC	16,068	234×	279×
brms S2Z NCP · Native Stan	603,455	1.18×	1.02×
brms S2Z NCP · WHMC	33,426	4.53×	5.15×
brms S2Z auto · Native Stan	109,652	31.8×	23.6×
brms S2Z auto · WHMC	17,512	151×	190×
BRM total NCP · WHMC	62,909	3.51×	4.03×
BRM total CP · WHMC	16,013	139×	166×
BRM total post-hoc position · WHMC	78,922	139×	33.7×
BRM total post-hoc gradient · WHMC	78,922	139×	33.7×
BRM total online position · WHMC	16,286	107×	128×
BRM total online gradient · WHMC	16,286	107×	128×

Both post-hoc losses select full centering for all six totals. They therefore produce the same fit as fixed CP, with the additional pilot cost. Both online losses also select CP and produce identical results to each other. S2Z CP + WHMC has the largest observed efficiency in this pilot; totals do not dominate this comparison.

Regional intercepts and slopes: 17 quantities

The common quantities are the two population coefficients, two group SDs, residual SD, six total intercepts and six total slopes.

Method	Total gradients	Sampling efficiency	Total efficiency
brms NCP · Native Stan	2,459,692	1×	1×
brms NCP · WHMC	1,351,111	1.1×	1.44×
brms S2Z CP · Native Stan	851,478	1.79×	1.61×
brms S2Z CP · WHMC	30,584	82.3×	95.5×
brms S2Z NCP · Native Stan	1,681,113	0.976×	0.912×
brms S2Z NCP · WHMC	485,549	1.55×	2×
brms S2Z auto · Native Stan	883,701	1.85×	1.64×
brms S2Z auto · WHMC	157,936	5.46×	7.07×
BRM total NCP · WHMC	1,178,130	0.577×	0.598×
BRM total CP · WHMC	27,183	142×	179×
BRM total post-hoc position · WHMC	1,215,894	22.9×	0.857×
BRM total post-hoc gradient · WHMC	1,202,655	122×	2.83×
BRM total online position · WHMC	53,202	14.3×	17.5×
BRM total online gradient · WHMC	20,137	130×	140×

The post-hoc workflows pay for a 1,178,130-gradient NCP pilot, so good retained-sampling efficiency can coexist with low total efficiency. Fixed CP and online adaptation avoid that separate pilot. Online position and gradient select different configurations here.

Saved-draw geometry

Each figure reuses the same 2,000 draws across its columns, with CP as the visualization baseline. Rows select distinct coordinates with minimum inferred centeredness, closest to 0.5, and maximum. The total figures use post-hoc position fits; their ACP column is checked against saved sampler coordinates. The S2Z figures use auto-WHMC fits and include the weighted shift needed to reconstruct the region-labelled Q*z coordinates. Those J contrast values represent J−1 independent directions.

Regional intercepts

Regional intercepts and slopes

The first set shows total coefficients; the S2Z set shows zero-sum contrasts. These are coordinate views of saved posterior draws, not predictive intervals or separate posterior fits.

Diagnostics, recovery and limits

The intercept-only ordinary NCP fits had 6 native and 68 WHMC divergences; S2Z NCP had 35 and 22, and total NCP had six. Other intercept-only arms had zero. Ordinary and S2Z NCP under WHMC also show the largest mean discrepancies from total CP: 4.52 and 5.89 combined estimated MCSEs for group SD. Their within-chain split R-hat maxima are 1.037 and 1.043. These weak arms need longer or replicated convergence checks.

For independent intercepts and slopes, ordinary NCP had 20 native and six WHMC divergences, S2Z NCP 37 and 60, total NCP eleven, and post-hoc gradient two. All other arms had zero. Native ordinary NCP hit maximum tree depth on 1,375 of 2,000 retained transitions. Mean differences from total CP are below three combined MCSEs except two S2Z-auto WHMC quantities, with a maximum of 3.56. These are descriptive checks across many quantities, not proof of convergence.

Removing the recovered population coefficients leaves nine invariant quantities for the intercept-only model and 15 for the independent model. The latter still gives 169× total efficiency for total CP, 132× for online gradient and 90.4× for S2Z CP + WHMC against that baseline scope. All twelve total-model minimum ESS values are unchanged across ten recovery seeds.

Population coefficients are recovered conditionally for both marginalized representations. Recovery adds genuine posterior variation and can increase ESS while also adding noise to a posterior-mean estimate. We retain per-quantity MCSE, compare the invariant quantities excluding stochastic recovery, and repeat total-model recovery with ten seeds. The same full scientific scope is used for every row in each main table.

Both automatic targets passed 25 independent density, gradient, coordinate and recovery checks. Ordinary brms and every S2Z setting passed their own density/gradient checks, including exact preservation of totals under generated-quantity recovery. BRM omits constant half-Student-t normalizers for bounded SD declarations; adding (K+1)*log(2) aligns the normalized densities without changing gradients or posteriors. Saved source-coordinate plots are audited against the retained fits.

The two rejected extreme Pathfinder trial points in the independent-total post-hoc fits were checked at higher precision: the conditional precision is mathematically positive definite, but its condition numbers exceed reliable double precision. The target adapter rejects those numerical proposals just as native Stan does. They do not change the target or the retained-draw scope.

This study does not establish a replicated ranking across chains/seeds, behavior of the original correlated-effects model, or controlled wall-time speedups.

Inspect and reproduce

• BRM models, priors and independent reference density.

• Density, gradient and recovery audit.

• Six total fits and ordinary WHMC baseline.

• Native Stan driver.

• S2Z WHMC and analysis.

• S2Z coordinate and generated-recovery audit.

• Saved-draw validation.

• Numerical proposal audit.

• Full raw-fit archive manifest.

• Reproduction instructions.

• cluster-intercept: comparison table.

• cluster-intercept: per-quantity ESS and MCSE.

• cluster-intercept: recovery sensitivity.

• cluster-intercept, ordinary_ncp: Stan.

• cluster-intercept, ordinary_ncp: data.

• cluster-intercept, ordinary_ncp: gradient counts.

• cluster-intercept, s2z_cp: Stan.

• cluster-intercept, s2z_cp: data.

• cluster-intercept, s2z_cp: gradient counts.

• cluster-intercept, s2z_ncp: Stan.

• cluster-intercept, s2z_ncp: data.

• cluster-intercept, s2z_ncp: gradient counts.

• cluster-intercept, s2z_auto: Stan.

• cluster-intercept, s2z_auto: data.

• cluster-intercept, s2z_auto: gradient counts.

• cluster-independent: comparison table.

• cluster-independent: per-quantity ESS and MCSE.

• cluster-independent: recovery sensitivity.

• cluster-independent, ordinary_ncp: Stan.

• cluster-independent, ordinary_ncp: data.

• cluster-independent, ordinary_ncp: gradient counts.

• cluster-independent, s2z_cp: Stan.

• cluster-independent, s2z_cp: data.

• cluster-independent, s2z_cp: gradient counts.

• cluster-independent, s2z_ncp: Stan.

• cluster-independent, s2z_ncp: data.

• cluster-independent, s2z_ncp: gradient counts.

• cluster-independent, s2z_auto: Stan.

• cluster-independent, s2z_auto: data.

• cluster-independent, s2z_auto: gradient counts.

• Automatic intercept_only: Stan.

• Automatic intercept_only: data.

• Automatic independent: Stan.

• Automatic independent: data.

Data revision: 2afb605f81cf0bdadc6476df865e0337aaacf183; brms PR #1919 revision: 73cf607889879cb2a55f50b88d8141d76ff43279; repaired WHMC: 7aed40b18bd4cdabb75330f285d5c9b355575ab9; Julia 1.10.11, CmdStan 2.39.0 and CmdStanR 0.9.0. BRM uses the built-in total planner with the retained-flat-prior correction. The source, raw fits, native CSVs, resolved weights and gradient counts are preserved in the linked manifests.