Wastewater-based Rt inference (EpiSewer / ww-inference-model)

This page ports the semi-mechanistic renewal model behind EpiSewer and the CDC ww-inference-model onto the @brm formula surface, and renders it in the standard four-pane view: the exact BRM authoring source, the StanBlocks model it emits, the generated Stan, and the selected Turing model. The build reads the declaration directly from the checked-in reproduction script, so the displayed BRM source cannot drift from the executable source.

A parallel StanBlocks-native port of the EpiSewer library is at StanBlocks — EpiSewer case study; the larger multi-stream structural port is on the CDC ww-inference page.

The model has three EpiSewer modules:

1. Infection model — a renewal process. Latent incidence follows I(t) = Rt(t) · Σ_{s} g(s) · I(t-s), seeded by a per-site initial level over the first window. g is the (fixed, known) generation-interval PMF, passed as data and captured into the kernel cell. The reproduction number is a smooth latent over time (here an hsgp prior on log Rt).

2. Shedding-load model — infections are convolved with a fecal shedding-load PMF sh (also data): load(t) = Σ_{s} sh(s) · I(t-s+1).

3. Measurement model — observed log-concentration is Normal on the log load plus a load→concentration scale: log C(t) ~ Normal(log load(t) + log_scale, sigma_obs).

Multiple sites share sigma_obs / log_scale and the smooth log Rt; each site carries its own seeding. The renewal and shedding recurrences are carried-state sequential scans — which plate deliberately cannot express — so they live in @deffun Stan functions (ww_renewal, ww_shed), called loop-free from the kernel(...) cell exactly as the PMX pages call ode_rk45.

The hidden setup below evaluates those two @deffun helpers and the fixtures from the checked-in source before the declaration is rendered.

The model

The generation-interval PMF g and the shedding PMF sh are ordinary data columns referenced by name inside the kernel(...) do-block; the kernel-cell data-vector capture registers them as shared Stan data, so they are re-bindable without editing the model — no baked-in constants. A per-timepoint smooth log Rt enters the cell as a ragged(log_rt, time_site) secondary axis; each site's seeding is a (1 | site) random effect.

The declaration below is extracted verbatim from wastewater_brm_model in the reviewed source. Its kernel(...) cell threads one ragged concentration-time course per site through the renewal and shedding scans.

Wastewater renewal Rt-inference (EpiSewer / ww-inference)

function wastewater_brm_model(data = wastewater_brm_fixture())
    @brm data begin
        sigma_obs ~ Exponential(1.0)                 # measurement noise (log scale)
        log_scale ~ Normal(0.0, 1.0)                 # load -> concentration scaling
        log_rt ~ 1 + hsgp(time_x; k = 10)            # smooth log-Rt over time
        log_I0 ~ 1 + (1 | site)                      # per-site seeding
        conc ~ kernel(t_grid, dv, ragged(log_rt, time_site), log_I0) do ts, yy, logRt_i, lI0
            infections = ww_renewal(logRt_i, g, exp(lI0))   # g captured as data
            load = ww_shed(infections, sh)                   # sh captured as data
            yy ~ normal(log(load) + log_scale, sigma_obs)
            load
        end
    end
end

BRMI:
  sigma_obs ~ Exponential(1.0)
  log_scale ~ Normal(0.0, 1.0)
  time_x: data (eltype=Float64, n=84)
  log_rt ~ 1 + hsgp(time_x; k=10)
  site: data (eltype=String, n=3)
  log_I0 ~ 1 + (1 | site)
  t_grid: data (eltype=Vector{Float64}, n=3)
  dv: data (eltype=Vector{Float64}, n=3)
  time_site: data (eltype=String, n=84)
  conc ~ kernel((ts, yy, logRt_i, lI0)->begin
        #= brm-docs-example.jl:8 =#
        infections = ww_renewal(logRt_i, g, exp(lI0))
        #= brm-docs-example.jl:9 =#
        load = ww_shed(infections, sh)
        #= brm-docs-example.jl:10 =#
        yy ~ normal(log(load) + log_scale, sigma_obs)
        #= brm-docs-example.jl:11 =#
        load
    end, t_grid, dv, ragged(log_rt, time_site), log_I0)
  g: data (eltype=Float64, n=7)
  sh: data (eltype=Float64, n=10)

SBBRMI with data keys = [:PHI_hsgp_time_x, :dv, :g, :kernel_conc_log_rt_ragged, :kernel_nsub_conc, :omega2_hsgp_time_x, :rho_lower_hsgp_time_x, :sh, :t_grid, :time_x, :total_A_log_I0, :total_group_log_I0, :total_location_log_I0, :total_ng_log_I0, :total_nk_log_I0, :total_np_log_I0, :total_precision_log_I0]
configured submodels:
_brm_total_scales_configured_1 = Base.merge(BayesianRegressionModels._brm_total_scales, quote
            tau ~ (ValueFamily(brm_vector_prior_faeb6f6956d0662a))(0.0, 1.0; n = 1)
        end)
emitted @slic body:
begin
    sigma_obs ~ exponential(1.0 ./ 1.0)
    log_scale ~ normal(0.0, 1.0)
    X_log_rt = hcat(rep_vector(1.0, num_elements(time_x)))
    pop_log_rt ~ popefs(; X = X_log_rt)
    hsgp_time_x ~ _sb_hsgp(; PHI = PHI_hsgp_time_x, omega2 = omega2_hsgp_time_x, rho_lower = rho_lower_hsgp_time_x)
    log_rt = pop_log_rt + hsgp_time_x
    total_scale_log_I0 ~ _brm_total_scales_configured_1(; n = total_nk_log_I0)
    total_log_I0::matrix[total_ng_log_I0, total_nk_log_I0] ~ brm_total(total_scale_log_I0, total_A_log_I0, total_location_log_I0, total_precision_log_I0)
    population_log_I0 = brm_total_recover_rng(total_log_I0, total_scale_log_I0, total_A_log_I0, total_location_log_I0, total_precision_log_I0)
    deviation_log_I0 = brm_total_deviations(total_log_I0, total_A_log_I0 * population_log_I0)
    total_Z_log_I0 = hcat(rep_vector(1.0, num_elements(total_group_log_I0)))
    log_I0 = rows_dot_product(total_log_I0[total_group_log_I0, :], total_Z_log_I0)
    conc ~ plate(t_grid, dv, kernel_conc_log_rt_ragged, log_I0; outer = (kernel_nsub_conc,)) do ts, yy, kernel_rows_logRt_i, lI0
            #= brm-docs-example.jl:8 =#
            infections = ww_renewal(log_rt[kernel_rows_logRt_i], g, exp(lI0))
            #= brm-docs-example.jl:9 =#
            load = ww_shed(infections, sh)
            #= brm-docs-example.jl:10 =#
            yy ~ normal(log(load) + log_scale, sigma_obs)
            #= brm-docs-example.jl:11 =#
            load
        end
end

functions {
matrix hcat(vector x) {
    int n = dims(x)[1];
    return to_matrix(x, n, 1);
}
vector brm_hsgp_sqrt_spd(
    matrix omega2,
    real sigma,
    vector rho
) {
    int m = dims(omega2)[1];
    int d = dims(omega2)[2];
    if (dims(rho)[1] != d) reject("brm_hsgp_sqrt_spd: dim mismatch — `rho` dim 1 (= ", dims(rho)[1], ") does not match `d` (= ", d, "), inferred from `omega2` dim 2. `d` sizes: `omega2` dim 2 (= ", dims(omega2)[2], "), `rho` dim 1 (= ", dims(rho)[1], ").");
    vector[m] rv;
    real scale = sigma;
    for(axis in 1:d) {
        scale *= sqrt((rho[axis] * 2.5066282746310002));
    }
    for(b in 1:m) {
        real exponent = 0.0;
        for(axis in 1:d) {
            exponent += (rho[axis] * rho[axis] * omega2[b, axis]);
        }
        rv[b] = (scale * exp((-0.25 * exponent)));
    }
    return rv;
}
// value UDF brm_vector_prior_faeb6f6956d0662a_lpdf
real brm_vector_prior_faeb6f6956d0662a_lpdf(
    vector x,
    real arg_1,
    real arg_2
) {
    if((x[1] < 0.0)) {
        return negative_infinity();
    }
    return lognormal_lpdf(x[1] | arg_1, arg_2);
}
real brm_total_lpdf(
    matrix total,
    vector tau,
    matrix A,
    vector location,
    vector precision
) {
    int j = dims(total)[1];
    int k = dims(total)[2];
    int p = dims(A)[2];
    if (dims(tau)[1] != k) reject("brm_total_lpdf: dim mismatch — `tau` dim 1 (= ", dims(tau)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
    if (dims(A)[1] != k) reject("brm_total_lpdf: dim mismatch — `A` dim 1 (= ", dims(A)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
    if (dims(location)[1] != p) reject("brm_total_lpdf: dim mismatch — `location` dim 1 (= ", dims(location)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], ").");
    if (dims(precision)[1] != p) reject("brm_total_lpdf: dim mismatch — `precision` dim 1 (= ", dims(precision)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], ").");
    matrix[dims(precision)[1], dims(precision)[1]] Q = brm_total_precision(tau, A, precision, j);
    vector[dims(total)[2]] average = brm_total_mean(total);
    vector[dims(precision)[1]] beta = brm_total_conditional_mean(total, tau, A, location, precision, Q);
    vector[dims(total)[2]] residual = (average - (A * beta));
    real quadratic = 0.0;
    real lp = ((-0.5 * ((j * k) - p) * 1.8378770664093453) - (j * sum(log(tau))));
    for(c in 1:k) {
        quadratic += ((j * square(residual[c])) / square(tau[c]));
        for(g in 1:j) {
            quadratic += (square((total[g, c] - average[c])) / square(tau[c]));
        }
    }
    for(a in 1:p) {
        if((precision[a] > 0.0)) {
            lp += (0.5 * (log(precision[a]) - 1.8378770664093453));
            quadratic += (precision[a] * square((beta[a] - location[a])));
        }
    }
    return (lp - (0.5 * (log_determinant(Q) + quadratic)));
}
matrix brm_total_precision(
    vector tau,
    matrix A,
    vector precision,
    int n_groups
) {
    int k = dims(tau)[1];
    int p = dims(A)[2];
    if (dims(A)[1] != k) reject("brm_total_precision: dim mismatch — `A` dim 1 (= ", dims(A)[1], ") does not match `k` (= ", k, "), inferred from `tau` dim 1. `k` sizes: `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
    if (dims(precision)[1] != p) reject("brm_total_precision: dim mismatch — `precision` dim 1 (= ", dims(precision)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `precision` dim 1 (= ", dims(precision)[1], ").");
    matrix[dims(precision)[1], dims(precision)[1]] out = diag_matrix(precision);
    for(a in 1:p) {
        for(b in 1:p) {
            for(c in 1:k) {
                out[a, b] += ((n_groups * A[c, a] * A[c, b]) / square(tau[c]));
            }
        }
    }
    return out;
}
vector brm_total_mean(
    matrix total
) {
    int j = dims(total)[1];
    int k = dims(total)[2];
    vector[k] out = rep_vector(0.0, k);
    for(c in 1:k) {
        out[c] = (sum(total[:, c]) / j);
    }
    return out;
}
vector brm_total_conditional_mean(
    matrix total,
    vector tau,
    matrix A,
    vector location,
    vector precision,
    matrix Q
) {
    int j = dims(total)[1];
    int k = dims(total)[2];
    int p = dims(A)[2];
    if (dims(tau)[1] != k) reject("brm_total_conditional_mean: dim mismatch — `tau` dim 1 (= ", dims(tau)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
    if (dims(A)[1] != k) reject("brm_total_conditional_mean: dim mismatch — `A` dim 1 (= ", dims(A)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
    if (dims(location)[1] != p) reject("brm_total_conditional_mean: dim mismatch — `location` dim 1 (= ", dims(location)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], "), `Q` dim 1 (= ", dims(Q)[1], "), `Q` dim 2 (= ", dims(Q)[2], ").");
    if (dims(precision)[1] != p) reject("brm_total_conditional_mean: dim mismatch — `precision` dim 1 (= ", dims(precision)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], "), `Q` dim 1 (= ", dims(Q)[1], "), `Q` dim 2 (= ", dims(Q)[2], ").");
    if (dims(Q)[1] != p) reject("brm_total_conditional_mean: dim mismatch — `Q` dim 1 (= ", dims(Q)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], "), `Q` dim 1 (= ", dims(Q)[1], "), `Q` dim 2 (= ", dims(Q)[2], ").");
    if (dims(Q)[2] != p) reject("brm_total_conditional_mean: dim mismatch — `Q` dim 2 (= ", dims(Q)[2], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], "), `Q` dim 1 (= ", dims(Q)[1], "), `Q` dim 2 (= ", dims(Q)[2], ").");
    vector[dims(total)[2]] average = brm_total_mean(total);
    vector[dims(location)[1]] natural = (precision .* location);
    for(a in 1:p) {
        for(c in 1:k) {
            natural[a] += ((j * A[c, a] * average[c]) / square(tau[c]));
        }
    }
    return mdivide_left_spd(Q, natural);
}
vector brm_total_recover_rng(
    matrix total,
    vector tau,
    matrix A,
    vector location,
    vector precision
) {
    int j = dims(total)[1];
    int k = dims(total)[2];
    int p = dims(A)[2];
    if (dims(tau)[1] != k) reject("brm_total_recover_rng: dim mismatch — `tau` dim 1 (= ", dims(tau)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
    if (dims(A)[1] != k) reject("brm_total_recover_rng: dim mismatch — `A` dim 1 (= ", dims(A)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
    if (dims(location)[1] != p) reject("brm_total_recover_rng: dim mismatch — `location` dim 1 (= ", dims(location)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], ").");
    if (dims(precision)[1] != p) reject("brm_total_recover_rng: dim mismatch — `precision` dim 1 (= ", dims(precision)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], ").");
    matrix[dims(precision)[1], dims(precision)[1]] Q = brm_total_precision(tau, A, precision, j);
    vector[dims(precision)[1]] beta = brm_total_conditional_mean(total, tau, A, location, precision, Q);
    return multi_normal_rng(beta, inverse_spd(Q));
}
matrix brm_total_deviations(
    matrix total,
    vector mu
) {
    int j = dims(total)[1];
    int k = dims(total)[2];
    if (dims(mu)[1] != k) reject("brm_total_deviations: dim mismatch — `mu` dim 1 (= ", dims(mu)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `mu` dim 1 (= ", dims(mu)[1], ").");
    matrix[dims(total)[1], dims(total)[2]] out = total;
    for(c in 1:k) {
        out[:, c] = (total[:, c] - rep_vector(mu[c], j));
    }
    return out;
}
int ragged_end(array[] int ends, int i) {
    return ends[i];
}
int ragged_start(
    array[] int ends,
    int i
) {
    if((i == 1)) {
        return 1;
    } else {
        return (1 + ends[(i - 1)]);
    }
}
vector ww_renewal(
    vector logRt,
    vector g,
    real I0
) {
    int nt = dims(logRt)[1];
    int ng = dims(g)[1];
    vector[nt] I;
    for(t in 1:nt) {
        real acc = 0.0;
        for(s in 1:ng) {
            real prev = (((t - s) >= 1) ? I[(t - s)] : I0);
            acc = (acc + (g[s] * prev));
        }
        I[t] = (exp(logRt[t]) * acc);
    }
    return I;
}
vector ww_shed(
    vector I,
    vector sh
) {
    int nt = dims(I)[1];
    int nsh = dims(sh)[1];
    vector[nt] load;
    for(t in 1:nt) {
        real acc = 0.0;
        for(s in 1:nsh) {
            int idx = ((t - s) + 1);
            real contrib = ((idx >= 1) ? I[idx] : 0.0);
            acc = (acc + (sh[s] * contrib));
        }
        load[t] = acc;
    }
    return load;
}
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 time_x_n;
    vector[time_x_n] time_x;
    int omega2_hsgp_time_x_m;
    int omega2_hsgp_time_x_n;
    real rho_lower_hsgp_time_x;
    matrix[omega2_hsgp_time_x_m, omega2_hsgp_time_x_n] omega2_hsgp_time_x;
    int PHI_hsgp_time_x_m;
    int PHI_hsgp_time_x_n;
    matrix[PHI_hsgp_time_x_m, PHI_hsgp_time_x_n] PHI_hsgp_time_x;
    int total_ng_log_I0;
    int total_nk_log_I0;
    int total_A_log_I0_m;
    int total_A_log_I0_n;
    matrix[total_A_log_I0_m, total_A_log_I0_n] total_A_log_I0;
    int total_location_log_I0_n;
    vector[total_location_log_I0_n] total_location_log_I0;
    int total_precision_log_I0_n;
    vector[total_precision_log_I0_n] total_precision_log_I0;
    int total_group_log_I0_n;
    array[total_group_log_I0_n] int total_group_log_I0;
    int kernel_nsub_conc;
    int kernel_conc_log_rt_ragged_ends_n;
    int kernel_conc_log_rt_ragged_mem_n;
    tuple(array[kernel_conc_log_rt_ragged_mem_n] int, array[kernel_conc_log_rt_ragged_ends_n] int) kernel_conc_log_rt_ragged;
    int dv_mem_n;
    int dv_ends_n;
    tuple(vector[dv_mem_n], array[dv_ends_n] int) dv;
    int g_n;
    vector[g_n] g;
    int sh_n;
    vector[sh_n] sh;
}
transformed data {
    matrix[num_elements(time_x), 1] X_log_rt = hcat(rep_vector(1.0, num_elements(time_x)));
    int pop_log_rt_n_covariates = 1;
    int hsgp_time_x_n_basis = omega2_hsgp_time_x_m;
    int hsgp_time_x_n_axes = omega2_hsgp_time_x_n;
    matrix[num_elements(total_group_log_I0), 1] total_Z_log_I0 = hcat(rep_vector(1.0, num_elements(total_group_log_I0)));
    array[kernel_nsub_conc] int conc_infections__pl_len_1;
    array[kernel_nsub_conc] int conc_load__pl_len_1;
    array[kernel_nsub_conc] int conc__pl_len_1;
    for(plate_i__pl_1 in 1:kernel_nsub_conc) {
        conc_infections__pl_len_1[plate_i__pl_1] = (
            1 +
            (
                ragged_end(kernel_conc_log_rt_ragged.2, plate_i__pl_1) -
                ragged_start(kernel_conc_log_rt_ragged.2, plate_i__pl_1)
            )
        );
        conc_load__pl_len_1[plate_i__pl_1] = (
            1 +
            (
                ragged_end(kernel_conc_log_rt_ragged.2, plate_i__pl_1) -
                ragged_start(kernel_conc_log_rt_ragged.2, plate_i__pl_1)
            )
        );
        conc__pl_len_1[plate_i__pl_1] = (
            1 +
            (
                ragged_end(kernel_conc_log_rt_ragged.2, plate_i__pl_1) -
                ragged_start(kernel_conc_log_rt_ragged.2, plate_i__pl_1)
            )
        );
    }
    array[kernel_nsub_conc] int conc_infections__pl_end_1 = cumulative_sum(conc_infections__pl_len_1);
    array[kernel_nsub_conc] int conc_load__pl_end_1 = cumulative_sum(conc_load__pl_len_1);
    array[kernel_nsub_conc] int conc__pl_end_1 = cumulative_sum(conc__pl_len_1);
}
parameters {
    real<lower=0.0> sigma_obs;
    real log_scale;
    vector[pop_log_rt_n_covariates] pop_log_rt_beta_pop;
    real<lower=rho_lower_hsgp_time_x> hsgp_time_x_rho_iso;
    real<lower=0.0> hsgp_time_x_sigma;
    vector[hsgp_time_x_n_basis] hsgp_time_x_beta_raw;
    vector<lower=0.0>[1] total_scale_log_I0_tau;
    matrix[total_ng_log_I0, total_nk_log_I0] total_log_I0;
}
transformed parameters {
    vector[num_elements(time_x)] pop_log_rt = (X_log_rt * pop_log_rt_beta_pop);
    vector[hsgp_time_x_n_axes] hsgp_time_x_rho = rep_vector(hsgp_time_x_rho_iso, hsgp_time_x_n_axes);
    vector[omega2_hsgp_time_x_m] hsgp_time_x_sqrt_spd = brm_hsgp_sqrt_spd(omega2_hsgp_time_x, hsgp_time_x_sigma, hsgp_time_x_rho);
    vector[PHI_hsgp_time_x_m] hsgp_time_x = (PHI_hsgp_time_x * (hsgp_time_x_sqrt_spd .* hsgp_time_x_beta_raw));
    vector[num_elements(time_x)] log_rt = (pop_log_rt + hsgp_time_x);
    vector<lower=0.0>[1] total_scale_log_I0 = total_scale_log_I0_tau;
    vector[num_elements(total_group_log_I0)] log_I0 = rows_dot_product(total_log_I0[total_group_log_I0, :], total_Z_log_I0);
    vector[sum(conc_infections__pl_len_1)] conc_infections__pl_mem_1;
    vector[sum(conc_load__pl_len_1)] conc_load__pl_mem_1;
    for(plate_i__pl_1 in 1:kernel_nsub_conc) {
        conc_infections__pl_mem_1[
            ragged_start(conc_infections__pl_end_1, plate_i__pl_1):ragged_end(conc_infections__pl_end_1, plate_i__pl_1)
        ] = ww_renewal(
            log_rt[
                kernel_conc_log_rt_ragged.1[
                    ragged_start(kernel_conc_log_rt_ragged.2, plate_i__pl_1):ragged_end(kernel_conc_log_rt_ragged.2, plate_i__pl_1)
                ]
            ],
            g,
            exp(log_I0[plate_i__pl_1])
        );
        conc_load__pl_mem_1[
            ragged_start(conc_load__pl_end_1, plate_i__pl_1):ragged_end(conc_load__pl_end_1, plate_i__pl_1)
        ] = ww_shed(
            conc_infections__pl_mem_1[
                ragged_start(conc_infections__pl_end_1, plate_i__pl_1):ragged_end(conc_infections__pl_end_1, plate_i__pl_1)
            ],
            sh
        );
    }
}
model {
    sigma_obs ~ exponential((1.0 ./ 1.0));
    log_scale ~ normal(0.0, 1.0);
    pop_log_rt_beta_pop ~ std_normal();
    hsgp_time_x_rho_iso ~ lognormal(0.0, 1.0);
    hsgp_time_x_sigma ~ lognormal(0.0, 1.0);
    hsgp_time_x_beta_raw ~ std_normal();
    total_scale_log_I0_tau ~ brm_vector_prior_faeb6f6956d0662a(0.0, 1.0);
    total_log_I0 ~ brm_total(total_scale_log_I0, total_A_log_I0, total_location_log_I0, total_precision_log_I0);
    for(plate_i__pl_1 in 1:kernel_nsub_conc) {
        dv.1[ragged_start(dv.2, plate_i__pl_1):ragged_end(dv.2, plate_i__pl_1)] ~ normal(
            (
                log(
                    conc_load__pl_mem_1[
                        ragged_start(conc_load__pl_end_1, plate_i__pl_1):ragged_end(conc_load__pl_end_1, plate_i__pl_1)
                    ]
                ) +
                log_scale
            ),
            sigma_obs
        );
    }
}
generated quantities {
    vector[total_precision_log_I0_n] population_log_I0 = brm_total_recover_rng(
        total_log_I0,
        total_scale_log_I0,
        total_A_log_I0,
        total_location_log_I0,
        total_precision_log_I0
    );
    matrix[total_ng_log_I0, total_A_log_I0_m] deviation_log_I0 = brm_total_deviations(total_log_I0, (total_A_log_I0 * population_log_I0));
    vector[sum(conc__pl_len_1)] conc__pl_mem_1;
    vector[num_elements(dv.1)] dv_gen;
    vector[num_elements(dv.2)] dv_likelihood;
    for(plate_i__pl_1 in 1:kernel_nsub_conc) {
        dv_gen[ragged_start(dv.2, plate_i__pl_1):ragged_end(dv.2, plate_i__pl_1)] = normal_vector_rng(
            (1 + (ragged_end(dv.2, plate_i__pl_1) - ragged_start(dv.2, plate_i__pl_1))),
            (
                log(
                    conc_load__pl_mem_1[
                        ragged_start(conc_load__pl_end_1, plate_i__pl_1):ragged_end(conc_load__pl_end_1, plate_i__pl_1)
                    ]
                ) +
                log_scale
            ),
            sigma_obs
        );
        dv_likelihood[plate_i__pl_1] = normal_lpdf(dv.1[ragged_start(dv.2, plate_i__pl_1):ragged_end(dv.2, plate_i__pl_1)] | 
            (
                log(
                    conc_load__pl_mem_1[
                        ragged_start(conc_load__pl_end_1, plate_i__pl_1):ragged_end(conc_load__pl_end_1, plate_i__pl_1)
                    ]
                ) +
                log_scale
            ),
            sigma_obs
        );
        conc__pl_mem_1[
            ragged_start(conc__pl_end_1, plate_i__pl_1):ragged_end(conc__pl_end_1, plate_i__pl_1)
        ] = conc_load__pl_mem_1[
            ragged_start(conc_load__pl_end_1, plate_i__pl_1):ragged_end(conc_load__pl_end_1, plate_i__pl_1)
        ];
    }
}

Turing unsupported for this BRM example

Turing backend: direct execution requires at least one observed likelihood

The Turing pane is intentionally retained even though this structural kernel is outside the current Turing executor; its build-time construction error is part of the comparison rather than being hidden.

Provenance and scope

The reproduction lives at research/wastewater/renewal.jl alongside its README. The same file also carries a pure-StanBlocks @slic form (wastewater_model / wastewater_censored_model, the latter left-censoring non-detects below a log limit of detection) — the layer where EpiSewer's own hand-written Stan lives, and where the generation-interval / shedding PMFs are likewise data. "Faithful" refers to the EpiSewer / ww-inference-model renewal structure; the fixed epidemiological PMFs and priors here are illustrative defaults, not a claim about any specific deployment's calibrated values.

The model artifact is verified — transpile, stanc, and finite BridgeStan density/gradient — by test/wastewater_model.jl (the @slic forms) and test/kernel_capture.jl (the @brm form). It is not sampled here; the deliverable is the model.

Run the reproduction

After bootstrapping the repository's test environment:

julia --startup-file=no --project=test test/wastewater_model.jl
julia --startup-file=no --project=test test/kernel_capture.jl

Set BRM_KERNEL_RUNTIME=0 to run lowering and stanc without BridgeStan instantiation.