Warfarin PK/PD: two-stage and joint models
This page presents two distinct contracts built from Sebastian Weber's public StanCon 2018 Warfarin program:
the faithful public two-stage workflow, where the PK posterior is reduced to fixed subject medians before the PD fit; and
a joint PK/PD model, where both likelihoods contribute to one posterior over shared latent PK effects.
Every displayed declaration uses the standard four-pane documentation view: the exact BRM authoring source, the StanBlocks model it emits, generated Stan, and the selected Turing model. The docs build reads the declarations directly from the checked-in reproduction scripts, so the displayed BRM source cannot drift from the executable source. Unsupported backends remain visible and show their current construction error.
The hidden setup below evaluates the shared public-data fixture, custom gamma2_overdisp distribution, analytical PK helpers, and turnover ODE helper before the three declarations are rendered.
Public two-stage workflow
The public reproduction is a sequential PK → PD workflow:
Fit the population PK model.
Reduce each subject's PK posterior to fixed medians for
log(tlag),log(ka), allometriclog(CL), and allometriclog(V).Fit the turnover-PD model conditional on those fixed PK values.
This is not the joint model
The PD likelihood does not update the PK posterior. The fixed PK posterior medians are inputs to the second fit, exactly as in the public two-stage source.
The checked-in reproduction includes the complete StanBlocks helper functions, custom observation family, public two-subject fixture, and executable stanc/BridgeStan acceptance harness. Read the executable source alongside its provenance and capability report. The reviewed two-stage source revision is 9fa8f77.
Stage 1: population PK
The PK stage is a one-compartment oral model with first-order absorption and a bounded lag. Clearance and volume carry fixed allometric weight effects as formula offset terms — log_cl0 ~ 1 + offset(0.75 * log_weight_ratio) + (1|cl_bsv|subject) and log_v0 ~ 1 + offset(log_weight_ratio) + … — rather than hand-wiring the weight scaling inside the kernel cell. Four independent subject effects modify lag, absorption, clearance, and volume.
The declaration below is extracted verbatim from warfarin_pk_brmi in the reviewed source. Its kernel(...) cell evaluates one ragged concentration-time course per subject.
brm-comparison
Public Warfarin population PKjulia
function warfarin_pk_brmi(data = warfarin_fixture())
return @brm data begin
sigma_pk ~ Normal(0.0, 2.0; lower=0.0)
kappa_pk ~ Gamma(0.2, 5.0)
lag_logit ~ 1 + (1 | tlag_bsv | subject)
log_ka ~ 1 + (1 | ka_bsv | subject)
# Allometric weight scaling as a formula OFFSET (fixed exponents 0.75 / 1.0),
# not hand-wired in the kernel cell: clearance ~ weight^0.75, volume ~ weight^1.
log_cl0 ~ 1 + offset(0.75 * log_weight_ratio) + (1 | cl_bsv | subject)
log_v0 ~ 1 + offset(log_weight_ratio) + (1 | v_bsv | subject)
effect(lag_logit, :) ~ Normal(0.0, 2.0)
effect(log_ka, :) ~ Normal(log(1.0), log(2.0) / 1.96)
effect(log_cl0, :) ~ Normal(log(0.1), log(10.0) / 1.96)
effect(log_v0, :) ~ Normal(log(10.0), log(10.0) / 1.96)
sd(:, tlag_bsv) ~ Normal(0.0, 0.5)
sd(:, ka_bsv) ~ Normal(0.0, 0.5)
sd(:, cl_bsv) ~ Normal(0.0, 0.5)
sd(:, v_bsv) ~ Normal(0.0, 0.5)
pk_pred ~ kernel(
pk_time, dose, pk_dv,
lag_logit, log_ka, log_cl0, log_v0,
) do times, dose_i, observed,
lag_i, lka_i, lcl_i, lv_i
# log_cl0 / log_v0 already carry the allometric offset (formula term above),
# so the cell uses them directly — no weight scaling hand-wired here.
log_tlag_i = log_inv_logit(lag_i)
prediction = exp(warfarin_pk_logconcentration(
times, log(dose_i), lka_i, lcl_i, lv_i, log_tlag_i,
)) + 1e-5
observed ~ warfarin_gamma2_overdisp(
prediction, sigma_pk, kappa_pk * 25.0,
)
prediction
end
end
endjulia
BRMI:
sigma_pk ~ Normal(0.0, 2.0; lower=0.0)
kappa_pk ~ Gamma(0.2, 5.0)
subject: data (eltype=Int64, n=2)
lag_logit ~ 1 + (1 | tlag_bsv | subject)
log_ka ~ 1 + (1 | ka_bsv | subject)
log_weight_ratio: data (eltype=Float64, n=2)
log_cl0 ~ 1 + offset((0.75 * log_weight_ratio)) + (1 | cl_bsv | subject)
log_v0 ~ 1 + offset(log_weight_ratio) + (1 | v_bsv | subject)
effect(lag_logit, :) ~ Normal(0.0, 2.0)
effect(log_ka, :) ~ Normal(log(1.0), /(log(2.0), 1.96))
effect(log_cl0, :) ~ Normal(log(0.1), /(log(10.0), 1.96))
effect(log_v0, :) ~ Normal(log(10.0), /(log(10.0), 1.96))
effect(sd, tlag_bsv) ~ Normal(0.0, 0.5)
effect(sd, ka_bsv) ~ Normal(0.0, 0.5)
effect(sd, cl_bsv) ~ Normal(0.0, 0.5)
effect(sd, v_bsv) ~ Normal(0.0, 0.5)
pk_time: data (eltype=Vector{Float64}, n=2)
dose: data (eltype=Float64, n=2)
pk_dv: data (eltype=Vector{Float64}, n=2)
pk_pred ~ kernel((times, dose_i, observed, lag_i, lka_i, lcl_i, lv_i)->begin
#= brm-docs-example.jl:29 =#
log_tlag_i = log_inv_logit(lag_i)
#= brm-docs-example.jl:30 =#
prediction = exp(warfarin_pk_logconcentration(times, log(dose_i), lka_i, lcl_i, lv_i, log_tlag_i)) + 1.0e-5
#= brm-docs-example.jl:33 =#
observed ~ warfarin_gamma2_overdisp(prediction, sigma_pk, kappa_pk * 25.0)
#= brm-docs-example.jl:36 =#
prediction
end, pk_time, dose, pk_dv, lag_logit, log_ka, log_cl0, log_v0)julia
SBBRMI with data keys = [:dose, :kernel_nsub_pk_pred, :log_weight_ratio, :n_subject, :n_terms_cl_bsv_subject, :n_terms_ka_bsv_subject, :n_terms_tlag_bsv_subject, :n_terms_v_bsv_subject, :pk_dv, :pk_time, :subject_idx]
configured submodels:
ranef_correlated_draws_generic_configured_1 = Base.merge(BayesianRegressionModels.ranef_correlated_draws_generic, quote
tau ~ normal(0.0, 0.5; n = n_terms, lower = 0.0)
end)
emitted @slic body:
begin
b_tlag_bsv_subject ~ ranef_correlated_draws_generic_configured_1(; group_idx = subject_idx, n_groups = n_subject, n_terms = n_terms_tlag_bsv_subject, lkj_eta = 1.0)
b_ka_bsv_subject ~ ranef_correlated_draws_generic_configured_1(; group_idx = subject_idx, n_groups = n_subject, n_terms = n_terms_ka_bsv_subject, lkj_eta = 1.0)
b_cl_bsv_subject ~ ranef_correlated_draws_generic_configured_1(; group_idx = subject_idx, n_groups = n_subject, n_terms = n_terms_cl_bsv_subject, lkj_eta = 1.0)
b_v_bsv_subject ~ ranef_correlated_draws_generic_configured_1(; group_idx = subject_idx, n_groups = n_subject, n_terms = n_terms_v_bsv_subject, lkj_eta = 1.0)
sigma_pk ~ normal(0.0, 2.0; lower = 0.0)
kappa_pk ~ gamma(0.2, 1.0 ./ 5.0)
X_lag_logit = hcat(rep_vector(1.0, num_elements(subject_idx)))
pop_lag_logit ~ _popefs_normal(; X = X_lag_logit, beta_loc = [0.0], beta_scale = [2.0])
r_lag_logit_tlag_bsv_subject = b_tlag_bsv_subject[subject_idx, 1]
lag_logit = pop_lag_logit + r_lag_logit_tlag_bsv_subject
X_log_ka = hcat(rep_vector(1.0, num_elements(subject_idx)))
pop_log_ka ~ _popefs_normal(; X = X_log_ka, beta_loc = [0.0], beta_scale = [0.35364652069384966])
r_log_ka_ka_bsv_subject = b_ka_bsv_subject[subject_idx, 1]
log_ka = pop_log_ka + r_log_ka_ka_bsv_subject
X_log_cl0 = hcat(rep_vector(1.0, num_elements(log_weight_ratio)))
pop_log_cl0 ~ _popefs_normal(; X = X_log_cl0, beta_loc = [-2.3025850929940455], beta_scale = [1.1747883127520642])
r_log_cl0_cl_bsv_subject = b_cl_bsv_subject[subject_idx, 1]
log_cl0 = pop_log_cl0 + 0.75 .* log_weight_ratio + r_log_cl0_cl_bsv_subject
X_log_v0 = hcat(rep_vector(1.0, num_elements(log_weight_ratio)))
pop_log_v0 ~ _popefs_normal(; X = X_log_v0, beta_loc = [2.302585092994046], beta_scale = [1.1747883127520642])
r_log_v0_v_bsv_subject = b_v_bsv_subject[subject_idx, 1]
log_v0 = pop_log_v0 + log_weight_ratio + r_log_v0_v_bsv_subject
pk_pred ~ plate(pk_time, dose, pk_dv, lag_logit, log_ka, log_cl0, log_v0; outer = (kernel_nsub_pk_pred,)) do times, dose_i, observed, lag_i, lka_i, lcl_i, lv_i
#= brm-docs-example.jl:29 =#
log_tlag_i = log_inv_logit(lag_i)
#= brm-docs-example.jl:30 =#
prediction = exp(warfarin_pk_logconcentration(times, log(dose_i), lka_i, lcl_i, lv_i, log_tlag_i)) + 1.0e-5
#= brm-docs-example.jl:33 =#
observed ~ warfarin_gamma2_overdisp(prediction, sigma_pk, kappa_pk * 25.0)
#= brm-docs-example.jl:36 =#
prediction
end
endstan
functions {
matrix hcat(vector x) {
int n = dims(x)[1];
return to_matrix(x, n, 1);
}
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 warfarin_pk_logconcentration(
vector tad,
real log_dose,
real log_ka,
real log_cl,
real log_v,
real log_tlag
) {
int n = dims(tad)[1];
real log_ke = (log_cl - log_v);
real log_delta = warfarin_log_diff_exp_abs(log_ka, log_ke);
real log_scale = (((log_dose - log_v) + log_ka) - log_delta);
vector[n] value;
for(i in 1:n) {
if((tad[i] <= exp(log_tlag))) {
value[i] = negative_infinity();
} else {
real log_tad = log((tad[i] - exp(log_tlag)));
real a = (-exp((log_ke + log_tad)));
real b = (-exp((log_ka + log_tad)));
value[i] = (log_scale + warfarin_log_diff_exp_abs(a, b));
}
}
return value;
}
real warfarin_log_diff_exp_abs(
real log_a,
real log_b
) {
return (0.5 * log_diff_exp(log_sum_exp((2.0 * log_a), (2.0 * log_b)), (log(2.0) + log_a + log_b)));
}
real warfarin_gamma2_overdisp_lpdf(
vector y,
vector mu,
real sigma,
real kappa
) {
int n = dims(y)[1];
if (dims(mu)[1] != n) reject("warfarin_gamma2_overdisp_lpdf: dim mismatch — `mu` dim 1 (= ", dims(mu)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `mu` dim 1 (= ", dims(mu)[1], ").");
real lp = 0.0;
for(i in 1:n) {
real variance = (square(sigma) + (square(mu[i]) / kappa));
real shape = (square(mu[i]) / variance);
real rate = (mu[i] / variance);
lp += gamma_lpdf(y[i] | shape, rate);
}
return lp;
}
vector warfarin_gamma2_overdisp_lpdfs(
vector y,
vector mu,
real sigma,
real kappa
) {
int n = dims(y)[1];
if (dims(mu)[1] != n) reject("warfarin_gamma2_overdisp_lpdfs: dim mismatch — `mu` dim 1 (= ", dims(mu)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `mu` dim 1 (= ", dims(mu)[1], ").");
vector[n] lp;
for(i in 1:n) {
real variance = (square(sigma) + (square(mu[i]) / kappa));
real shape = (square(mu[i]) / variance);
real rate = (mu[i] / variance);
lp[i] = gamma_lpdf(y[i] | shape, rate);
}
return lp;
}
vector warfarin_gamma2_overdisp_vector_rng(
int anontok__1,
vector mu,
real sigma,
real kappa
) {
int n = anontok__1;
if (dims(mu)[1] != n) reject("warfarin_gamma2_overdisp_rng: dim mismatch — `mu` dim 1 (= ", dims(mu)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `mu` dim 1 (= ", dims(mu)[1], ").");
vector[n] draw;
for(i in 1:n) {
draw[i] = warfarin_gamma2_overdisp_rng(mu[i], sigma, kappa);
}
return draw;
}
real warfarin_gamma2_overdisp_rng(
real mu,
real sigma,
real kappa
) {
real variance = (square(sigma) + (square(mu) / kappa));
real shape = (square(mu) / variance);
real rate = (mu / variance);
return gamma_rng(shape, rate);
}
}
data {
int n_terms_tlag_bsv_subject;
int n_subject;
int n_terms_ka_bsv_subject;
int n_terms_cl_bsv_subject;
int n_terms_v_bsv_subject;
int subject_idx_n;
array[subject_idx_n] int subject_idx;
int log_weight_ratio_n;
vector[log_weight_ratio_n] log_weight_ratio;
int kernel_nsub_pk_pred;
int pk_time_ends_n;
int pk_time_mem_n;
tuple(vector[pk_time_mem_n], array[pk_time_ends_n] int) pk_time;
int pk_dv_mem_n;
int pk_dv_ends_n;
tuple(vector[pk_dv_mem_n], array[pk_dv_ends_n] int) pk_dv;
int dose_n;
vector[dose_n] dose;
}
transformed data {
matrix[num_elements(subject_idx), 1] X_lag_logit = hcat(rep_vector(1.0, num_elements(subject_idx)));
int pop_lag_logit_n_covariates = 1;
matrix[num_elements(subject_idx), 1] X_log_ka = hcat(rep_vector(1.0, num_elements(subject_idx)));
int pop_log_ka_n_covariates = 1;
matrix[num_elements(log_weight_ratio), 1] X_log_cl0 = hcat(rep_vector(1.0, num_elements(log_weight_ratio)));
int pop_log_cl0_n_covariates = 1;
matrix[num_elements(log_weight_ratio), 1] X_log_v0 = hcat(rep_vector(1.0, num_elements(log_weight_ratio)));
int pop_log_v0_n_covariates = 1;
array[kernel_nsub_pk_pred] int pk_pred_prediction__pl_len_1;
array[kernel_nsub_pk_pred] int pk_pred__pl_len_1;
for(plate_i__pl_1 in 1:kernel_nsub_pk_pred) {
pk_pred_prediction__pl_len_1[plate_i__pl_1] = (1 + (ragged_end(pk_time.2, plate_i__pl_1) - ragged_start(pk_time.2, plate_i__pl_1)));
pk_pred__pl_len_1[plate_i__pl_1] = (1 + (ragged_end(pk_time.2, plate_i__pl_1) - ragged_start(pk_time.2, plate_i__pl_1)));
}
array[kernel_nsub_pk_pred] int pk_pred_prediction__pl_end_1 = cumulative_sum(pk_pred_prediction__pl_len_1);
array[kernel_nsub_pk_pred] int pk_pred__pl_end_1 = cumulative_sum(pk_pred__pl_len_1);
}
parameters {
cholesky_factor_corr[n_terms_tlag_bsv_subject] b_tlag_bsv_subject_L;
vector<lower=0.0>[n_terms_tlag_bsv_subject] b_tlag_bsv_subject_tau;
vector[(n_terms_tlag_bsv_subject * n_subject)] b_tlag_bsv_subject_z_flat;
cholesky_factor_corr[n_terms_ka_bsv_subject] b_ka_bsv_subject_L;
vector<lower=0.0>[n_terms_ka_bsv_subject] b_ka_bsv_subject_tau;
vector[(n_terms_ka_bsv_subject * n_subject)] b_ka_bsv_subject_z_flat;
cholesky_factor_corr[n_terms_cl_bsv_subject] b_cl_bsv_subject_L;
vector<lower=0.0>[n_terms_cl_bsv_subject] b_cl_bsv_subject_tau;
vector[(n_terms_cl_bsv_subject * n_subject)] b_cl_bsv_subject_z_flat;
cholesky_factor_corr[n_terms_v_bsv_subject] b_v_bsv_subject_L;
vector<lower=0.0>[n_terms_v_bsv_subject] b_v_bsv_subject_tau;
vector[(n_terms_v_bsv_subject * n_subject)] b_v_bsv_subject_z_flat;
real<lower=0.0> sigma_pk;
real<lower=0.0> kappa_pk;
vector[pop_lag_logit_n_covariates] pop_lag_logit_beta_pop;
vector[pop_log_ka_n_covariates] pop_log_ka_beta_pop;
vector[pop_log_cl0_n_covariates] pop_log_cl0_beta_pop;
vector[pop_log_v0_n_covariates] pop_log_v0_beta_pop;
}
transformed parameters {
matrix[n_terms_tlag_bsv_subject, n_subject] b_tlag_bsv_subject_z = to_matrix(b_tlag_bsv_subject_z_flat, n_terms_tlag_bsv_subject, n_subject);
matrix[n_subject, n_terms_tlag_bsv_subject] b_tlag_bsv_subject = ((diag_pre_multiply(b_tlag_bsv_subject_tau, b_tlag_bsv_subject_L) * b_tlag_bsv_subject_z)');
matrix[n_terms_ka_bsv_subject, n_subject] b_ka_bsv_subject_z = to_matrix(b_ka_bsv_subject_z_flat, n_terms_ka_bsv_subject, n_subject);
matrix[n_subject, n_terms_ka_bsv_subject] b_ka_bsv_subject = ((diag_pre_multiply(b_ka_bsv_subject_tau, b_ka_bsv_subject_L) * b_ka_bsv_subject_z)');
matrix[n_terms_cl_bsv_subject, n_subject] b_cl_bsv_subject_z = to_matrix(b_cl_bsv_subject_z_flat, n_terms_cl_bsv_subject, n_subject);
matrix[n_subject, n_terms_cl_bsv_subject] b_cl_bsv_subject = ((diag_pre_multiply(b_cl_bsv_subject_tau, b_cl_bsv_subject_L) * b_cl_bsv_subject_z)');
matrix[n_terms_v_bsv_subject, n_subject] b_v_bsv_subject_z = to_matrix(b_v_bsv_subject_z_flat, n_terms_v_bsv_subject, n_subject);
matrix[n_subject, n_terms_v_bsv_subject] b_v_bsv_subject = ((diag_pre_multiply(b_v_bsv_subject_tau, b_v_bsv_subject_L) * b_v_bsv_subject_z)');
vector[num_elements(subject_idx)] pop_lag_logit = (X_lag_logit * pop_lag_logit_beta_pop);
vector[subject_idx_n] r_lag_logit_tlag_bsv_subject = b_tlag_bsv_subject[subject_idx, 1];
vector[num_elements(subject_idx)] lag_logit = (pop_lag_logit + r_lag_logit_tlag_bsv_subject);
vector[num_elements(subject_idx)] pop_log_ka = (X_log_ka * pop_log_ka_beta_pop);
vector[subject_idx_n] r_log_ka_ka_bsv_subject = b_ka_bsv_subject[subject_idx, 1];
vector[num_elements(subject_idx)] log_ka = (pop_log_ka + r_log_ka_ka_bsv_subject);
vector[num_elements(log_weight_ratio)] pop_log_cl0 = (X_log_cl0 * pop_log_cl0_beta_pop);
vector[subject_idx_n] r_log_cl0_cl_bsv_subject = b_cl_bsv_subject[subject_idx, 1];
vector[num_elements(log_weight_ratio)] log_cl0 = (pop_log_cl0 + (0.75 .* log_weight_ratio) + r_log_cl0_cl_bsv_subject);
vector[num_elements(log_weight_ratio)] pop_log_v0 = (X_log_v0 * pop_log_v0_beta_pop);
vector[subject_idx_n] r_log_v0_v_bsv_subject = b_v_bsv_subject[subject_idx, 1];
vector[num_elements(log_weight_ratio)] log_v0 = (pop_log_v0 + log_weight_ratio + r_log_v0_v_bsv_subject);
real pk_pred__pl_inv1_1 = (kappa_pk * 25.0);
vector[sum(pk_pred_prediction__pl_len_1)] pk_pred_prediction__pl_mem_1;
vector[kernel_nsub_pk_pred] pk_pred_log_tlag_i;
for(plate_i__pl_1 in 1:kernel_nsub_pk_pred) {
pk_pred_log_tlag_i[plate_i__pl_1] = log_inv_logit(lag_logit[plate_i__pl_1]);
pk_pred_prediction__pl_mem_1[
ragged_start(pk_pred_prediction__pl_end_1, plate_i__pl_1):ragged_end(pk_pred_prediction__pl_end_1, plate_i__pl_1)
] = (
exp(
warfarin_pk_logconcentration(
pk_time.1[ragged_start(pk_time.2, plate_i__pl_1):ragged_end(pk_time.2, plate_i__pl_1)],
log(dose[plate_i__pl_1]),
log_ka[plate_i__pl_1],
log_cl0[plate_i__pl_1],
log_v0[plate_i__pl_1],
pk_pred_log_tlag_i[plate_i__pl_1]
)
) +
1.0e-5
);
}
}
model {
b_tlag_bsv_subject_L ~ lkj_corr_cholesky(1.0);
b_tlag_bsv_subject_tau ~ normal(0.0, 0.5);
b_tlag_bsv_subject_z_flat ~ std_normal();
b_ka_bsv_subject_L ~ lkj_corr_cholesky(1.0);
b_ka_bsv_subject_tau ~ normal(0.0, 0.5);
b_ka_bsv_subject_z_flat ~ std_normal();
b_cl_bsv_subject_L ~ lkj_corr_cholesky(1.0);
b_cl_bsv_subject_tau ~ normal(0.0, 0.5);
b_cl_bsv_subject_z_flat ~ std_normal();
b_v_bsv_subject_L ~ lkj_corr_cholesky(1.0);
b_v_bsv_subject_tau ~ normal(0.0, 0.5);
b_v_bsv_subject_z_flat ~ std_normal();
sigma_pk ~ normal(0.0, 2.0);
kappa_pk ~ gamma(0.2, (1.0 ./ 5.0));
pop_lag_logit_beta_pop ~ normal([0.0]', [2.0]');
pop_log_ka_beta_pop ~ normal([0.0]', [0.35364652069384966]');
pop_log_cl0_beta_pop ~ normal([-2.3025850929940455]', [1.1747883127520642]');
pop_log_v0_beta_pop ~ normal([2.302585092994046]', [1.1747883127520642]');
for(plate_i__pl_1 in 1:kernel_nsub_pk_pred) {
pk_dv.1[ragged_start(pk_dv.2, plate_i__pl_1):ragged_end(pk_dv.2, plate_i__pl_1)] ~ warfarin_gamma2_overdisp(
pk_pred_prediction__pl_mem_1[
ragged_start(pk_pred_prediction__pl_end_1, plate_i__pl_1):ragged_end(pk_pred_prediction__pl_end_1, plate_i__pl_1)
],
sigma_pk,
pk_pred__pl_inv1_1
);
}
}
generated quantities {
vector[sum(pk_pred__pl_len_1)] pk_pred__pl_mem_1;
vector[num_elements(pk_dv.1)] pk_dv_gen;
vector[num_elements(pk_dv.2)] pk_dv_likelihood;
for(plate_i__pl_1 in 1:kernel_nsub_pk_pred) {
pk_dv_gen[ragged_start(pk_dv.2, plate_i__pl_1):ragged_end(pk_dv.2, plate_i__pl_1)] = warfarin_gamma2_overdisp_vector_rng(
(1 + (ragged_end(pk_dv.2, plate_i__pl_1) - ragged_start(pk_dv.2, plate_i__pl_1))),
pk_pred_prediction__pl_mem_1[
ragged_start(pk_pred_prediction__pl_end_1, plate_i__pl_1):ragged_end(pk_pred_prediction__pl_end_1, plate_i__pl_1)
],
sigma_pk,
pk_pred__pl_inv1_1
);
pk_dv_likelihood[plate_i__pl_1] = warfarin_gamma2_overdisp_lpdf(pk_dv.1[ragged_start(pk_dv.2, plate_i__pl_1):ragged_end(pk_dv.2, plate_i__pl_1)] |
pk_pred_prediction__pl_mem_1[
ragged_start(pk_pred_prediction__pl_end_1, plate_i__pl_1):ragged_end(pk_pred_prediction__pl_end_1, plate_i__pl_1)
],
sigma_pk,
pk_pred__pl_inv1_1
);
pk_pred__pl_mem_1[
ragged_start(pk_pred__pl_end_1, plate_i__pl_1):ragged_end(pk_pred__pl_end_1, plate_i__pl_1)
] = pk_pred_prediction__pl_mem_1[
ragged_start(pk_pred_prediction__pl_end_1, plate_i__pl_1):ragged_end(pk_pred_prediction__pl_end_1, plate_i__pl_1)
];
}
}julia
Turing unsupported for this BRM example
Turing backend: direct execution requires at least one observed likelihoodThe model preserves the source priors, tlagMax = 1, the 0.75 clearance and 1.0 volume allometric exponents, and the PK overdispersion multiplier 5² = 25.
The fixed-median handoff
The second fit receives four values per subject:
pk_log_tlagpk_log_kapk_log_clpk_log_v
They are posterior medians from the PK fit, already transformed and—with clearance and volume—already including the weight effects. They are ordinary fixed data columns in the PD model, not sampled parameters shared between the two likelihoods. The small executable fixture uses the corresponding medians from the public PD model's Stan data dump.
Stage 2: turnover PD
The PD stage solves a turnover ODE conditional on those fixed PK medians. Three independent subject effects modify baseline response, inverse turnover rate, and EC50.
brm-comparison
Public Warfarin turnover PDjulia
function warfarin_pd_brmi(data = warfarin_fixture())
return @brm data begin
sigma_pd ~ Normal(0.0, 10.0; lower=0.0)
kappa_pd ~ Gamma(0.2, 5.0)
log_r0 ~ 1 + (1 | r0_bsv | subject)
log_inv_kout ~ 1 + (1 | kout_bsv | subject)
log_ec50 ~ 1 + (1 | ec50_bsv | subject)
effect(log_r0, :) ~ Normal(log(80.0), log(10.0) / 1.96)
effect(log_inv_kout, :) ~ Normal(log(30.0), log(10.0) / 1.96)
effect(log_ec50, :) ~ Normal(log(2.5), log(10.0) / 1.96)
sd(:, r0_bsv) ~ Normal(0.0, 0.5)
sd(:, kout_bsv) ~ Normal(0.0, 0.5)
sd(:, ec50_bsv) ~ Normal(0.0, 0.5)
pd_pred ~ kernel(
pd_time, dose, pd_dv,
pk_log_tlag, pk_log_ka, pk_log_cl, pk_log_v,
log_r0, log_inv_kout, log_ec50,
) do times, dose_i, observed,
log_tlag_i, log_ka_i, log_cl_i, log_v_i,
lr0_i, linvkout_i, lec50_i
prediction = warfarin_turnover_prediction(
times, log(dose_i), log_ka_i, log_cl_i, log_v_i,
log_tlag_i, lr0_i, linvkout_i, lec50_i,
)
observed ~ warfarin_gamma2_overdisp(
prediction, sigma_pd, kappa_pd * 625.0,
)
prediction
end
end
endjulia
BRMI:
sigma_pd ~ Normal(0.0, 10.0; lower=0.0)
kappa_pd ~ Gamma(0.2, 5.0)
subject: data (eltype=Int64, n=2)
log_r0 ~ 1 + (1 | r0_bsv | subject)
log_inv_kout ~ 1 + (1 | kout_bsv | subject)
log_ec50 ~ 1 + (1 | ec50_bsv | subject)
effect(log_r0, :) ~ Normal(log(80.0), /(log(10.0), 1.96))
effect(log_inv_kout, :) ~ Normal(log(30.0), /(log(10.0), 1.96))
effect(log_ec50, :) ~ Normal(log(2.5), /(log(10.0), 1.96))
effect(sd, r0_bsv) ~ Normal(0.0, 0.5)
effect(sd, kout_bsv) ~ Normal(0.0, 0.5)
effect(sd, ec50_bsv) ~ Normal(0.0, 0.5)
pd_time: data (eltype=Vector{Float64}, n=2)
dose: data (eltype=Float64, n=2)
pd_dv: data (eltype=Vector{Float64}, n=2)
pk_log_tlag: data (eltype=Float64, n=2)
pk_log_ka: data (eltype=Float64, n=2)
pk_log_cl: data (eltype=Float64, n=2)
pk_log_v: data (eltype=Float64, n=2)
pd_pred ~ kernel((times, dose_i, observed, log_tlag_i, log_ka_i, log_cl_i, log_v_i, lr0_i, linvkout_i, lec50_i)->begin
#= brm-docs-example.jl:24 =#
prediction = warfarin_turnover_prediction(times, log(dose_i), log_ka_i, log_cl_i, log_v_i, log_tlag_i, lr0_i, linvkout_i, lec50_i)
#= brm-docs-example.jl:28 =#
observed ~ warfarin_gamma2_overdisp(prediction, sigma_pd, kappa_pd * 625.0)
#= brm-docs-example.jl:31 =#
prediction
end, pd_time, dose, pd_dv, pk_log_tlag, pk_log_ka, pk_log_cl, pk_log_v, log_r0, log_inv_kout, log_ec50)julia
SBBRMI with data keys = [:dose, :kernel_nsub_pd_pred, :pd_dv, :pd_time, :pk_log_cl, :pk_log_ka, :pk_log_tlag, :pk_log_v, :subject, :total_A_log_ec50, :total_A_log_inv_kout, :total_A_log_r0, :total_group_log_ec50, :total_group_log_inv_kout, :total_group_log_r0, :total_location_log_ec50, :total_location_log_inv_kout, :total_location_log_r0, :total_ng_log_ec50, :total_ng_log_inv_kout, :total_ng_log_r0, :total_nk_log_ec50, :total_nk_log_inv_kout, :total_nk_log_r0, :total_np_log_ec50, :total_np_log_inv_kout, :total_np_log_r0, :total_precision_log_ec50, :total_precision_log_inv_kout, :total_precision_log_r0]
configured submodels:
_brm_total_scales_configured_1 = Base.merge(BayesianRegressionModels._brm_total_scales, quote
tau ~ (ValueFamily(brm_vector_prior_ca4b8a1c1bc116d6))(0.0, 0.5; n = 1)
end)
emitted @slic body:
begin
sigma_pd ~ normal(0.0, 10.0; lower = 0.0)
kappa_pd ~ gamma(0.2, 1.0 ./ 5.0)
total_scale_log_r0 ~ _brm_total_scales_configured_1(; n = total_nk_log_r0)
total_log_r0::matrix[total_ng_log_r0, total_nk_log_r0] ~ brm_total(total_scale_log_r0, total_A_log_r0, total_location_log_r0, total_precision_log_r0)
population_log_r0 = brm_total_recover_rng(total_log_r0, total_scale_log_r0, total_A_log_r0, total_location_log_r0, total_precision_log_r0)
deviation_log_r0 = brm_total_deviations(total_log_r0, total_A_log_r0 * population_log_r0)
total_Z_log_r0 = hcat(rep_vector(1.0, num_elements(total_group_log_r0)))
log_r0 = rows_dot_product(total_log_r0[total_group_log_r0, :], total_Z_log_r0)
total_scale_log_inv_kout ~ _brm_total_scales_configured_1(; n = total_nk_log_inv_kout)
total_log_inv_kout::matrix[total_ng_log_inv_kout, total_nk_log_inv_kout] ~ brm_total(total_scale_log_inv_kout, total_A_log_inv_kout, total_location_log_inv_kout, total_precision_log_inv_kout)
population_log_inv_kout = brm_total_recover_rng(total_log_inv_kout, total_scale_log_inv_kout, total_A_log_inv_kout, total_location_log_inv_kout, total_precision_log_inv_kout)
deviation_log_inv_kout = brm_total_deviations(total_log_inv_kout, total_A_log_inv_kout * population_log_inv_kout)
total_Z_log_inv_kout = hcat(rep_vector(1.0, num_elements(total_group_log_inv_kout)))
log_inv_kout = rows_dot_product(total_log_inv_kout[total_group_log_inv_kout, :], total_Z_log_inv_kout)
total_scale_log_ec50 ~ _brm_total_scales_configured_1(; n = total_nk_log_ec50)
total_log_ec50::matrix[total_ng_log_ec50, total_nk_log_ec50] ~ brm_total(total_scale_log_ec50, total_A_log_ec50, total_location_log_ec50, total_precision_log_ec50)
population_log_ec50 = brm_total_recover_rng(total_log_ec50, total_scale_log_ec50, total_A_log_ec50, total_location_log_ec50, total_precision_log_ec50)
deviation_log_ec50 = brm_total_deviations(total_log_ec50, total_A_log_ec50 * population_log_ec50)
total_Z_log_ec50 = hcat(rep_vector(1.0, num_elements(total_group_log_ec50)))
log_ec50 = rows_dot_product(total_log_ec50[total_group_log_ec50, :], total_Z_log_ec50)
pd_pred ~ plate(pd_time, dose, pd_dv, pk_log_tlag, pk_log_ka, pk_log_cl, pk_log_v, log_r0, log_inv_kout, log_ec50; outer = (kernel_nsub_pd_pred,)) do times, dose_i, observed, log_tlag_i, log_ka_i, log_cl_i, log_v_i, lr0_i, linvkout_i, lec50_i
#= brm-docs-example.jl:24 =#
prediction = warfarin_turnover_prediction(times, log(dose_i), log_ka_i, log_cl_i, log_v_i, log_tlag_i, lr0_i, linvkout_i, lec50_i)
#= brm-docs-example.jl:28 =#
observed ~ warfarin_gamma2_overdisp(prediction, sigma_pd, kappa_pd * 625.0)
#= brm-docs-example.jl:31 =#
prediction
end
endstan
functions {
// value UDF brm_vector_prior_ca4b8a1c1bc116d6_lpdf
real brm_vector_prior_ca4b8a1c1bc116d6_lpdf(
vector x,
real arg_1,
real arg_2
) {
if((x[1] < 0.0)) {
return negative_infinity();
}
return normal_lpdf(x[1] | arg_1, arg_2);
}
real brm_total_lpdf(
matrix total,
vector tau,
matrix A,
vector location,
vector precision
) {
int j = dims(total)[1];
int k = dims(total)[2];
int p = dims(A)[2];
if (dims(tau)[1] != k) reject("brm_total_lpdf: dim mismatch — `tau` dim 1 (= ", dims(tau)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
if (dims(A)[1] != k) reject("brm_total_lpdf: dim mismatch — `A` dim 1 (= ", dims(A)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
if (dims(location)[1] != p) reject("brm_total_lpdf: dim mismatch — `location` dim 1 (= ", dims(location)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], ").");
if (dims(precision)[1] != p) reject("brm_total_lpdf: dim mismatch — `precision` dim 1 (= ", dims(precision)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], ").");
matrix[dims(precision)[1], dims(precision)[1]] Q = brm_total_precision(tau, A, precision, j);
vector[dims(total)[2]] average = brm_total_mean(total);
vector[dims(precision)[1]] beta = brm_total_conditional_mean(total, tau, A, location, precision, Q);
vector[dims(total)[2]] residual = (average - (A * beta));
real quadratic = 0.0;
real lp = ((-0.5 * ((j * k) - p) * 1.8378770664093453) - (j * sum(log(tau))));
for(c in 1:k) {
quadratic += ((j * square(residual[c])) / square(tau[c]));
for(g in 1:j) {
quadratic += (square((total[g, c] - average[c])) / square(tau[c]));
}
}
for(a in 1:p) {
if((precision[a] > 0.0)) {
lp += (0.5 * (log(precision[a]) - 1.8378770664093453));
quadratic += (precision[a] * square((beta[a] - location[a])));
}
}
return (lp - (0.5 * (log_determinant(Q) + quadratic)));
}
matrix brm_total_precision(
vector tau,
matrix A,
vector precision,
int n_groups
) {
int k = dims(tau)[1];
int p = dims(A)[2];
if (dims(A)[1] != k) reject("brm_total_precision: dim mismatch — `A` dim 1 (= ", dims(A)[1], ") does not match `k` (= ", k, "), inferred from `tau` dim 1. `k` sizes: `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
if (dims(precision)[1] != p) reject("brm_total_precision: dim mismatch — `precision` dim 1 (= ", dims(precision)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `precision` dim 1 (= ", dims(precision)[1], ").");
matrix[dims(precision)[1], dims(precision)[1]] out = diag_matrix(precision);
for(a in 1:p) {
for(b in 1:p) {
for(c in 1:k) {
out[a, b] += ((n_groups * A[c, a] * A[c, b]) / square(tau[c]));
}
}
}
return out;
}
vector brm_total_mean(
matrix total
) {
int j = dims(total)[1];
int k = dims(total)[2];
vector[k] out = rep_vector(0.0, k);
for(c in 1:k) {
out[c] = (sum(total[:, c]) / j);
}
return out;
}
vector brm_total_conditional_mean(
matrix total,
vector tau,
matrix A,
vector location,
vector precision,
matrix Q
) {
int j = dims(total)[1];
int k = dims(total)[2];
int p = dims(A)[2];
if (dims(tau)[1] != k) reject("brm_total_conditional_mean: dim mismatch — `tau` dim 1 (= ", dims(tau)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
if (dims(A)[1] != k) reject("brm_total_conditional_mean: dim mismatch — `A` dim 1 (= ", dims(A)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
if (dims(location)[1] != p) reject("brm_total_conditional_mean: dim mismatch — `location` dim 1 (= ", dims(location)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], "), `Q` dim 1 (= ", dims(Q)[1], "), `Q` dim 2 (= ", dims(Q)[2], ").");
if (dims(precision)[1] != p) reject("brm_total_conditional_mean: dim mismatch — `precision` dim 1 (= ", dims(precision)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], "), `Q` dim 1 (= ", dims(Q)[1], "), `Q` dim 2 (= ", dims(Q)[2], ").");
if (dims(Q)[1] != p) reject("brm_total_conditional_mean: dim mismatch — `Q` dim 1 (= ", dims(Q)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], "), `Q` dim 1 (= ", dims(Q)[1], "), `Q` dim 2 (= ", dims(Q)[2], ").");
if (dims(Q)[2] != p) reject("brm_total_conditional_mean: dim mismatch — `Q` dim 2 (= ", dims(Q)[2], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], "), `Q` dim 1 (= ", dims(Q)[1], "), `Q` dim 2 (= ", dims(Q)[2], ").");
vector[dims(total)[2]] average = brm_total_mean(total);
vector[dims(location)[1]] natural = (precision .* location);
for(a in 1:p) {
for(c in 1:k) {
natural[a] += ((j * A[c, a] * average[c]) / square(tau[c]));
}
}
return mdivide_left_spd(Q, natural);
}
vector brm_total_recover_rng(
matrix total,
vector tau,
matrix A,
vector location,
vector precision
) {
int j = dims(total)[1];
int k = dims(total)[2];
int p = dims(A)[2];
if (dims(tau)[1] != k) reject("brm_total_recover_rng: dim mismatch — `tau` dim 1 (= ", dims(tau)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
if (dims(A)[1] != k) reject("brm_total_recover_rng: dim mismatch — `A` dim 1 (= ", dims(A)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `tau` dim 1 (= ", dims(tau)[1], "), `A` dim 1 (= ", dims(A)[1], ").");
if (dims(location)[1] != p) reject("brm_total_recover_rng: dim mismatch — `location` dim 1 (= ", dims(location)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], ").");
if (dims(precision)[1] != p) reject("brm_total_recover_rng: dim mismatch — `precision` dim 1 (= ", dims(precision)[1], ") does not match `p` (= ", p, "), inferred from `A` dim 2. `p` sizes: `A` dim 2 (= ", dims(A)[2], "), `location` dim 1 (= ", dims(location)[1], "), `precision` dim 1 (= ", dims(precision)[1], ").");
matrix[dims(precision)[1], dims(precision)[1]] Q = brm_total_precision(tau, A, precision, j);
vector[dims(precision)[1]] beta = brm_total_conditional_mean(total, tau, A, location, precision, Q);
return multi_normal_rng(beta, inverse_spd(Q));
}
matrix brm_total_deviations(
matrix total,
vector mu
) {
int j = dims(total)[1];
int k = dims(total)[2];
if (dims(mu)[1] != k) reject("brm_total_deviations: dim mismatch — `mu` dim 1 (= ", dims(mu)[1], ") does not match `k` (= ", k, "), inferred from `total` dim 2. `k` sizes: `total` dim 2 (= ", dims(total)[2], "), `mu` dim 1 (= ", dims(mu)[1], ").");
matrix[dims(total)[1], dims(total)[2]] out = total;
for(c in 1:k) {
out[:, c] = (total[:, c] - rep_vector(mu[c], j));
}
return out;
}
matrix hcat(vector x) {
int n = dims(x)[1];
return to_matrix(x, n, 1);
}
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 warfarin_turnover_prediction(
vector times,
real log_dose,
real log_ka,
real log_cl,
real log_v,
real log_tlag,
real log_r0,
real log_inv_kout,
real log_ec50
) {
int n = dims(times)[1];
vector[1] initial_state = rep_vector(exp(log_r0), 1);
array[dims(times)[1]] vector[1] trajectory = ode_rk45_tol(
warfarin_turnover_rhs,
initial_state,
-0.0001,
to_array_1d(times),
1.0e-5,
0.001,
500,
log_dose,
log_ka,
log_cl,
log_v,
log_tlag,
log_r0,
log_inv_kout,
log_ec50
);
return to_vector(trajectory[:, 1]);
}
vector warfarin_turnover_rhs(
real time,
vector state,
real log_dose,
real log_ka,
real log_cl,
real log_v,
real log_tlag,
real log_r0,
real log_inv_kout,
real log_ec50
) {
int m = dims(state)[1];
real log_conc = warfarin_pk_logconcentration_one(time, log_dose, log_ka, log_cl, log_v, log_tlag);
real log_kout = (-log_inv_kout);
real log_kin = (log_r0 + log_kout);
real log_inhibition = log_inv_logit((log_conc - log_ec50));
vector[m] derivative;
derivative[1] = (exp((log_kin + log1m_exp(log_inhibition))) - (state[1] * exp(log_kout)));
return derivative;
}
real warfarin_pk_logconcentration_one(
real time,
real log_dose,
real log_ka,
real log_cl,
real log_v,
real log_tlag
) {
if((time < exp(log_tlag))) {
return -25.0;
}
real log_ke = (log_cl - log_v);
real log_delta = warfarin_log_diff_exp_abs(log_ka, log_ke);
real log_scale = (((log_dose - log_v) + log_ka) - log_delta);
real log_tad = log((time - exp(log_tlag)));
real a = (-exp((log_ke + log_tad)));
real b = (-exp((log_ka + log_tad)));
return (log_scale + warfarin_log_diff_exp_abs(a, b));
}
real warfarin_log_diff_exp_abs(
real log_a,
real log_b
) {
return (0.5 * log_diff_exp(log_sum_exp((2.0 * log_a), (2.0 * log_b)), (log(2.0) + log_a + log_b)));
}
real warfarin_gamma2_overdisp_lpdf(
vector y,
vector mu,
real sigma,
real kappa
) {
int n = dims(y)[1];
if (dims(mu)[1] != n) reject("warfarin_gamma2_overdisp_lpdf: dim mismatch — `mu` dim 1 (= ", dims(mu)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `mu` dim 1 (= ", dims(mu)[1], ").");
real lp = 0.0;
for(i in 1:n) {
real variance = (square(sigma) + (square(mu[i]) / kappa));
real shape = (square(mu[i]) / variance);
real rate = (mu[i] / variance);
lp += gamma_lpdf(y[i] | shape, rate);
}
return lp;
}
vector warfarin_gamma2_overdisp_lpdfs(
vector y,
vector mu,
real sigma,
real kappa
) {
int n = dims(y)[1];
if (dims(mu)[1] != n) reject("warfarin_gamma2_overdisp_lpdfs: dim mismatch — `mu` dim 1 (= ", dims(mu)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `mu` dim 1 (= ", dims(mu)[1], ").");
vector[n] lp;
for(i in 1:n) {
real variance = (square(sigma) + (square(mu[i]) / kappa));
real shape = (square(mu[i]) / variance);
real rate = (mu[i] / variance);
lp[i] = gamma_lpdf(y[i] | shape, rate);
}
return lp;
}
vector warfarin_gamma2_overdisp_vector_rng(
int anontok__1,
vector mu,
real sigma,
real kappa
) {
int n = anontok__1;
if (dims(mu)[1] != n) reject("warfarin_gamma2_overdisp_rng: dim mismatch — `mu` dim 1 (= ", dims(mu)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `mu` dim 1 (= ", dims(mu)[1], ").");
vector[n] draw;
for(i in 1:n) {
draw[i] = warfarin_gamma2_overdisp_rng(mu[i], sigma, kappa);
}
return draw;
}
real warfarin_gamma2_overdisp_rng(
real mu,
real sigma,
real kappa
) {
real variance = (square(sigma) + (square(mu) / kappa));
real shape = (square(mu) / variance);
real rate = (mu / variance);
return gamma_rng(shape, rate);
}
}
data {
int total_ng_log_r0;
int total_nk_log_r0;
int total_A_log_r0_m;
int total_A_log_r0_n;
matrix[total_A_log_r0_m, total_A_log_r0_n] total_A_log_r0;
int total_location_log_r0_n;
vector[total_location_log_r0_n] total_location_log_r0;
int total_precision_log_r0_n;
vector[total_precision_log_r0_n] total_precision_log_r0;
int total_group_log_r0_n;
array[total_group_log_r0_n] int total_group_log_r0;
int total_ng_log_inv_kout;
int total_nk_log_inv_kout;
int total_A_log_inv_kout_m;
int total_A_log_inv_kout_n;
matrix[total_A_log_inv_kout_m, total_A_log_inv_kout_n] total_A_log_inv_kout;
int total_location_log_inv_kout_n;
vector[total_location_log_inv_kout_n] total_location_log_inv_kout;
int total_precision_log_inv_kout_n;
vector[total_precision_log_inv_kout_n] total_precision_log_inv_kout;
int total_group_log_inv_kout_n;
array[total_group_log_inv_kout_n] int total_group_log_inv_kout;
int total_ng_log_ec50;
int total_nk_log_ec50;
int total_A_log_ec50_m;
int total_A_log_ec50_n;
matrix[total_A_log_ec50_m, total_A_log_ec50_n] total_A_log_ec50;
int total_location_log_ec50_n;
vector[total_location_log_ec50_n] total_location_log_ec50;
int total_precision_log_ec50_n;
vector[total_precision_log_ec50_n] total_precision_log_ec50;
int total_group_log_ec50_n;
array[total_group_log_ec50_n] int total_group_log_ec50;
int kernel_nsub_pd_pred;
int pd_time_ends_n;
int pd_time_mem_n;
tuple(vector[pd_time_mem_n], array[pd_time_ends_n] int) pd_time;
int pd_dv_mem_n;
int pd_dv_ends_n;
tuple(vector[pd_dv_mem_n], array[pd_dv_ends_n] int) pd_dv;
int dose_n;
vector[dose_n] dose;
int pk_log_ka_n;
vector[pk_log_ka_n] pk_log_ka;
int pk_log_cl_n;
vector[pk_log_cl_n] pk_log_cl;
int pk_log_v_n;
vector[pk_log_v_n] pk_log_v;
int pk_log_tlag_n;
vector[pk_log_tlag_n] pk_log_tlag;
}
transformed data {
matrix[num_elements(total_group_log_r0), 1] total_Z_log_r0 = hcat(rep_vector(1.0, num_elements(total_group_log_r0)));
matrix[num_elements(total_group_log_inv_kout), 1] total_Z_log_inv_kout = hcat(rep_vector(1.0, num_elements(total_group_log_inv_kout)));
matrix[num_elements(total_group_log_ec50), 1] total_Z_log_ec50 = hcat(rep_vector(1.0, num_elements(total_group_log_ec50)));
array[kernel_nsub_pd_pred] int pd_pred__pl_len_1;
array[kernel_nsub_pd_pred] int pd_pred_prediction__pl_len_1;
for(plate_i__pl_1 in 1:kernel_nsub_pd_pred) {
pd_pred__pl_len_1[plate_i__pl_1] = (1 + (ragged_end(pd_time.2, plate_i__pl_1) - ragged_start(pd_time.2, plate_i__pl_1)));
pd_pred_prediction__pl_len_1[plate_i__pl_1] = (1 + (ragged_end(pd_time.2, plate_i__pl_1) - ragged_start(pd_time.2, plate_i__pl_1)));
}
array[kernel_nsub_pd_pred] int pd_pred__pl_end_1 = cumulative_sum(pd_pred__pl_len_1);
array[kernel_nsub_pd_pred] int pd_pred_prediction__pl_end_1 = cumulative_sum(pd_pred_prediction__pl_len_1);
}
parameters {
real<lower=0.0> sigma_pd;
real<lower=0.0> kappa_pd;
vector<lower=0.0>[1] total_scale_log_r0_tau;
matrix[total_ng_log_r0, total_nk_log_r0] total_log_r0;
vector<lower=0.0>[1] total_scale_log_inv_kout_tau;
matrix[total_ng_log_inv_kout, total_nk_log_inv_kout] total_log_inv_kout;
vector<lower=0.0>[1] total_scale_log_ec50_tau;
matrix[total_ng_log_ec50, total_nk_log_ec50] total_log_ec50;
}
transformed parameters {
vector<lower=0.0>[1] total_scale_log_r0 = total_scale_log_r0_tau;
vector[num_elements(total_group_log_r0)] log_r0 = rows_dot_product(total_log_r0[total_group_log_r0, :], total_Z_log_r0);
vector<lower=0.0>[1] total_scale_log_inv_kout = total_scale_log_inv_kout_tau;
vector[num_elements(total_group_log_inv_kout)] log_inv_kout = rows_dot_product(total_log_inv_kout[total_group_log_inv_kout, :], total_Z_log_inv_kout);
vector<lower=0.0>[1] total_scale_log_ec50 = total_scale_log_ec50_tau;
vector[num_elements(total_group_log_ec50)] log_ec50 = rows_dot_product(total_log_ec50[total_group_log_ec50, :], total_Z_log_ec50);
real pd_pred__pl_inv1_1 = (kappa_pd * 625.0);
vector[sum(pd_pred_prediction__pl_len_1)] pd_pred_prediction__pl_mem_1;
for(plate_i__pl_1 in 1:kernel_nsub_pd_pred) {
pd_pred_prediction__pl_mem_1[
ragged_start(pd_pred_prediction__pl_end_1, plate_i__pl_1):ragged_end(pd_pred_prediction__pl_end_1, plate_i__pl_1)
] = warfarin_turnover_prediction(
pd_time.1[ragged_start(pd_time.2, plate_i__pl_1):ragged_end(pd_time.2, plate_i__pl_1)],
log(dose[plate_i__pl_1]),
pk_log_ka[plate_i__pl_1],
pk_log_cl[plate_i__pl_1],
pk_log_v[plate_i__pl_1],
pk_log_tlag[plate_i__pl_1],
log_r0[plate_i__pl_1],
log_inv_kout[plate_i__pl_1],
log_ec50[plate_i__pl_1]
);
}
}
model {
sigma_pd ~ normal(0.0, 10.0);
kappa_pd ~ gamma(0.2, (1.0 ./ 5.0));
total_scale_log_r0_tau ~ brm_vector_prior_ca4b8a1c1bc116d6(0.0, 0.5);
total_log_r0 ~ brm_total(total_scale_log_r0, total_A_log_r0, total_location_log_r0, total_precision_log_r0);
total_scale_log_inv_kout_tau ~ brm_vector_prior_ca4b8a1c1bc116d6(0.0, 0.5);
total_log_inv_kout ~ brm_total(
total_scale_log_inv_kout,
total_A_log_inv_kout,
total_location_log_inv_kout,
total_precision_log_inv_kout
);
total_scale_log_ec50_tau ~ brm_vector_prior_ca4b8a1c1bc116d6(0.0, 0.5);
total_log_ec50 ~ brm_total(total_scale_log_ec50, total_A_log_ec50, total_location_log_ec50, total_precision_log_ec50);
for(plate_i__pl_1 in 1:kernel_nsub_pd_pred) {
pd_dv.1[ragged_start(pd_dv.2, plate_i__pl_1):ragged_end(pd_dv.2, plate_i__pl_1)] ~ warfarin_gamma2_overdisp(
pd_pred_prediction__pl_mem_1[
ragged_start(pd_pred_prediction__pl_end_1, plate_i__pl_1):ragged_end(pd_pred_prediction__pl_end_1, plate_i__pl_1)
],
sigma_pd,
pd_pred__pl_inv1_1
);
}
}
generated quantities {
vector[total_precision_log_r0_n] population_log_r0 = brm_total_recover_rng(
total_log_r0,
total_scale_log_r0,
total_A_log_r0,
total_location_log_r0,
total_precision_log_r0
);
matrix[total_ng_log_r0, total_A_log_r0_m] deviation_log_r0 = brm_total_deviations(total_log_r0, (total_A_log_r0 * population_log_r0));
vector[total_precision_log_inv_kout_n] population_log_inv_kout = brm_total_recover_rng(
total_log_inv_kout,
total_scale_log_inv_kout,
total_A_log_inv_kout,
total_location_log_inv_kout,
total_precision_log_inv_kout
);
matrix[total_ng_log_inv_kout, total_A_log_inv_kout_m] deviation_log_inv_kout = brm_total_deviations(total_log_inv_kout, (total_A_log_inv_kout * population_log_inv_kout));
vector[total_precision_log_ec50_n] population_log_ec50 = brm_total_recover_rng(
total_log_ec50,
total_scale_log_ec50,
total_A_log_ec50,
total_location_log_ec50,
total_precision_log_ec50
);
matrix[total_ng_log_ec50, total_A_log_ec50_m] deviation_log_ec50 = brm_total_deviations(total_log_ec50, (total_A_log_ec50 * population_log_ec50));
vector[sum(pd_pred__pl_len_1)] pd_pred__pl_mem_1;
vector[num_elements(pd_dv.1)] pd_dv_gen;
vector[num_elements(pd_dv.2)] pd_dv_likelihood;
for(plate_i__pl_1 in 1:kernel_nsub_pd_pred) {
pd_dv_gen[ragged_start(pd_dv.2, plate_i__pl_1):ragged_end(pd_dv.2, plate_i__pl_1)] = warfarin_gamma2_overdisp_vector_rng(
(1 + (ragged_end(pd_dv.2, plate_i__pl_1) - ragged_start(pd_dv.2, plate_i__pl_1))),
pd_pred_prediction__pl_mem_1[
ragged_start(pd_pred_prediction__pl_end_1, plate_i__pl_1):ragged_end(pd_pred_prediction__pl_end_1, plate_i__pl_1)
],
sigma_pd,
pd_pred__pl_inv1_1
);
pd_dv_likelihood[plate_i__pl_1] = warfarin_gamma2_overdisp_lpdf(pd_dv.1[ragged_start(pd_dv.2, plate_i__pl_1):ragged_end(pd_dv.2, plate_i__pl_1)] |
pd_pred_prediction__pl_mem_1[
ragged_start(pd_pred_prediction__pl_end_1, plate_i__pl_1):ragged_end(pd_pred_prediction__pl_end_1, plate_i__pl_1)
],
sigma_pd,
pd_pred__pl_inv1_1
);
pd_pred__pl_mem_1[
ragged_start(pd_pred__pl_end_1, plate_i__pl_1):ragged_end(pd_pred__pl_end_1, plate_i__pl_1)
] = pd_pred_prediction__pl_mem_1[
ragged_start(pd_pred_prediction__pl_end_1, plate_i__pl_1):ragged_end(pd_pred_prediction__pl_end_1, plate_i__pl_1)
];
}
}julia
Turing unsupported for this BRM example
Turing backend: direct execution requires at least one observed likelihoodThe helper called by the kernel uses RK45 with t₀ = -1e-4, relative tolerance 1e-5, absolute tolerance 1e-3, and at most 500 steps. The PD overdispersion multiplier is 25² = 625.
Both stages retain the source observation variance sigma² + mu² / kappa. The reviewed executable passed lowering, stanc, and finite BridgeStan density/gradient checks for both models.
Joint PK/PD: one posterior
The joint model combines the two public likelihoods in a single BRM declaration. It samples each subject's latent lag, absorption, clearance, and volume effects once and feeds those same quantities into both the PK concentration kernel and the concentration-driven PD turnover ODE.
That shared latent path creates bidirectional updating: PK uncertainty propagates into the PD model, while PD observations can update the posterior for the PK quantities. There is no fixed-median handoff and the joint model does not consume pk_log_tlag, pk_log_ka, pk_log_cl, or pk_log_v.
The PK and PD observations still use independent residual distributions. They are measured on different time grids and in different units; adding residual covariance would require a new alignment and measurement model, not merely a shared posterior. The coupling here is through the shared latent PK effects.
The joint declaration and acceptance harness were reviewed at canonical commit f64cf6d291bf40565a6e73299623ea61ead34aa3. See the executable joint source and its joint-model contract.
brm-comparison
Joint Warfarin PK/PDjulia
function warfarin_joint_brmi(data = warfarin_fixture())
return @brm data begin
sigma_pk ~ Normal(0.0, 2.0; lower=0.0)
kappa_pk ~ Gamma(0.2, 5.0)
sigma_pd ~ Normal(0.0, 10.0; lower=0.0)
kappa_pd ~ Gamma(0.2, 5.0)
lag_logit ~ 1 + (1 | tlag_bsv | subject)
log_ka ~ 1 + (1 | ka_bsv | subject)
# Allometric weight scaling as a formula OFFSET (fixed exponents 0.75 / 1.0),
# shared by both the PK and PD cells below — not hand-wired in either cell.
log_cl0 ~ 1 + offset(0.75 * log_weight_ratio) + (1 | cl_bsv | subject)
log_v0 ~ 1 + offset(log_weight_ratio) + (1 | v_bsv | subject)
log_r0 ~ 1 + (1 | r0_bsv | subject)
log_inv_kout ~ 1 + (1 | kout_bsv | subject)
log_ec50 ~ 1 + (1 | ec50_bsv | subject)
effect(lag_logit, :) ~ Normal(0.0, 2.0)
effect(log_ka, :) ~ Normal(log(1.0), log(2.0) / 1.96)
effect(log_cl0, :) ~ Normal(log(0.1), log(10.0) / 1.96)
effect(log_v0, :) ~ Normal(log(10.0), log(10.0) / 1.96)
effect(log_r0, :) ~ Normal(log(80.0), log(10.0) / 1.96)
effect(log_inv_kout, :) ~ Normal(log(30.0), log(10.0) / 1.96)
effect(log_ec50, :) ~ Normal(log(2.5), log(10.0) / 1.96)
sd(:, tlag_bsv) ~ Normal(0.0, 0.5)
sd(:, ka_bsv) ~ Normal(0.0, 0.5)
sd(:, cl_bsv) ~ Normal(0.0, 0.5)
sd(:, v_bsv) ~ Normal(0.0, 0.5)
sd(:, r0_bsv) ~ Normal(0.0, 0.5)
sd(:, kout_bsv) ~ Normal(0.0, 0.5)
sd(:, ec50_bsv) ~ Normal(0.0, 0.5)
pk_pred ~ kernel(
pk_time, dose, pk_dv,
lag_logit, log_ka, log_cl0, log_v0,
) do times, dose_i, observed,
lag_i, lka_i, lcl_i, lv_i
log_tlag_i = log_inv_logit(lag_i)
prediction = exp(warfarin_pk_logconcentration(
times, log(dose_i), lka_i, lcl_i, lv_i, log_tlag_i,
)) + 1e-5
pk_pointwise_loglik = warfarin_gamma2_overdisp_lpdfs(
observed, prediction, sigma_pk, kappa_pk * 25.0,
)
observed ~ warfarin_gamma2_overdisp(
prediction, sigma_pk, kappa_pk * 25.0,
)
prediction
end
pd_pred ~ kernel(
pd_time, dose, pd_dv,
lag_logit, log_ka, log_cl0, log_v0,
log_r0, log_inv_kout, log_ec50,
) do times, dose_i, observed,
lag_i, lka_i, lcl_i, lv_i,
lr0_i, linvkout_i, lec50_i
log_tlag_i = log_inv_logit(lag_i)
prediction = warfarin_turnover_prediction(
times, log(dose_i), lka_i, lcl_i, lv_i,
log_tlag_i, lr0_i, linvkout_i, lec50_i,
)
pd_pointwise_loglik = warfarin_gamma2_overdisp_lpdfs(
observed, prediction, sigma_pd, kappa_pd * 625.0,
)
observed ~ warfarin_gamma2_overdisp(
prediction, sigma_pd, kappa_pd * 625.0,
)
prediction
end
end
endjulia
BRMI:
sigma_pk ~ Normal(0.0, 2.0; lower=0.0)
kappa_pk ~ Gamma(0.2, 5.0)
sigma_pd ~ Normal(0.0, 10.0; lower=0.0)
kappa_pd ~ Gamma(0.2, 5.0)
subject: data (eltype=Int64, n=2)
lag_logit ~ 1 + (1 | tlag_bsv | subject)
log_ka ~ 1 + (1 | ka_bsv | subject)
log_weight_ratio: data (eltype=Float64, n=2)
log_cl0 ~ 1 + offset((0.75 * log_weight_ratio)) + (1 | cl_bsv | subject)
log_v0 ~ 1 + offset(log_weight_ratio) + (1 | v_bsv | subject)
log_r0 ~ 1 + (1 | r0_bsv | subject)
log_inv_kout ~ 1 + (1 | kout_bsv | subject)
log_ec50 ~ 1 + (1 | ec50_bsv | subject)
effect(lag_logit, :) ~ Normal(0.0, 2.0)
effect(log_ka, :) ~ Normal(log(1.0), /(log(2.0), 1.96))
effect(log_cl0, :) ~ Normal(log(0.1), /(log(10.0), 1.96))
effect(log_v0, :) ~ Normal(log(10.0), /(log(10.0), 1.96))
effect(log_r0, :) ~ Normal(log(80.0), /(log(10.0), 1.96))
effect(log_inv_kout, :) ~ Normal(log(30.0), /(log(10.0), 1.96))
effect(log_ec50, :) ~ Normal(log(2.5), /(log(10.0), 1.96))
effect(sd, tlag_bsv) ~ Normal(0.0, 0.5)
effect(sd, ka_bsv) ~ Normal(0.0, 0.5)
effect(sd, cl_bsv) ~ Normal(0.0, 0.5)
effect(sd, v_bsv) ~ Normal(0.0, 0.5)
effect(sd, r0_bsv) ~ Normal(0.0, 0.5)
effect(sd, kout_bsv) ~ Normal(0.0, 0.5)
effect(sd, ec50_bsv) ~ Normal(0.0, 0.5)
pk_time: data (eltype=Vector{Float64}, n=2)
dose: data (eltype=Float64, n=2)
pk_dv: data (eltype=Vector{Float64}, n=2)
pk_pred ~ kernel((times, dose_i, observed, lag_i, lka_i, lcl_i, lv_i)->begin
#= brm-docs-example.jl:38 =#
log_tlag_i = log_inv_logit(lag_i)
#= brm-docs-example.jl:39 =#
prediction = exp(warfarin_pk_logconcentration(times, log(dose_i), lka_i, lcl_i, lv_i, log_tlag_i)) + 1.0e-5
#= brm-docs-example.jl:42 =#
pk_pointwise_loglik = warfarin_gamma2_overdisp_lpdfs(observed, prediction, sigma_pk, kappa_pk * 25.0)
#= brm-docs-example.jl:45 =#
observed ~ warfarin_gamma2_overdisp(prediction, sigma_pk, kappa_pk * 25.0)
#= brm-docs-example.jl:48 =#
prediction
end, pk_time, dose, pk_dv, lag_logit, log_ka, log_cl0, log_v0)
pd_time: data (eltype=Vector{Float64}, n=2)
pd_dv: data (eltype=Vector{Float64}, n=2)
pd_pred ~ kernel((times, dose_i, observed, lag_i, lka_i, lcl_i, lv_i, lr0_i, linvkout_i, lec50_i)->begin
#= brm-docs-example.jl:58 =#
log_tlag_i = log_inv_logit(lag_i)
#= brm-docs-example.jl:59 =#
prediction = warfarin_turnover_prediction(times, log(dose_i), lka_i, lcl_i, lv_i, log_tlag_i, lr0_i, linvkout_i, lec50_i)
#= brm-docs-example.jl:63 =#
pd_pointwise_loglik = warfarin_gamma2_overdisp_lpdfs(observed, prediction, sigma_pd, kappa_pd * 625.0)
#= brm-docs-example.jl:66 =#
observed ~ warfarin_gamma2_overdisp(prediction, sigma_pd, kappa_pd * 625.0)
#= brm-docs-example.jl:69 =#
prediction
end, pd_time, dose, pd_dv, lag_logit, log_ka, log_cl0, log_v0, log_r0, log_inv_kout, log_ec50)julia
SBBRMI with data keys = [:dose, :kernel_nsub_pd_pred, :kernel_nsub_pk_pred, :log_weight_ratio, :n_subject, :n_terms_cl_bsv_subject, :n_terms_ec50_bsv_subject, :n_terms_ka_bsv_subject, :n_terms_kout_bsv_subject, :n_terms_r0_bsv_subject, :n_terms_tlag_bsv_subject, :n_terms_v_bsv_subject, :pd_dv, :pd_time, :pk_dv, :pk_time, :subject_idx]
configured submodels:
ranef_correlated_draws_generic_configured_1 = Base.merge(BayesianRegressionModels.ranef_correlated_draws_generic, quote
tau ~ normal(0.0, 0.5; n = n_terms, lower = 0.0)
end)
emitted @slic body:
begin
b_tlag_bsv_subject ~ ranef_correlated_draws_generic_configured_1(; group_idx = subject_idx, n_groups = n_subject, n_terms = n_terms_tlag_bsv_subject, lkj_eta = 1.0)
b_ka_bsv_subject ~ ranef_correlated_draws_generic_configured_1(; group_idx = subject_idx, n_groups = n_subject, n_terms = n_terms_ka_bsv_subject, lkj_eta = 1.0)
b_cl_bsv_subject ~ ranef_correlated_draws_generic_configured_1(; group_idx = subject_idx, n_groups = n_subject, n_terms = n_terms_cl_bsv_subject, lkj_eta = 1.0)
b_v_bsv_subject ~ ranef_correlated_draws_generic_configured_1(; group_idx = subject_idx, n_groups = n_subject, n_terms = n_terms_v_bsv_subject, lkj_eta = 1.0)
b_r0_bsv_subject ~ ranef_correlated_draws_generic_configured_1(; group_idx = subject_idx, n_groups = n_subject, n_terms = n_terms_r0_bsv_subject, lkj_eta = 1.0)
b_kout_bsv_subject ~ ranef_correlated_draws_generic_configured_1(; group_idx = subject_idx, n_groups = n_subject, n_terms = n_terms_kout_bsv_subject, lkj_eta = 1.0)
b_ec50_bsv_subject ~ ranef_correlated_draws_generic_configured_1(; group_idx = subject_idx, n_groups = n_subject, n_terms = n_terms_ec50_bsv_subject, lkj_eta = 1.0)
sigma_pk ~ normal(0.0, 2.0; lower = 0.0)
kappa_pk ~ gamma(0.2, 1.0 ./ 5.0)
sigma_pd ~ normal(0.0, 10.0; lower = 0.0)
kappa_pd ~ gamma(0.2, 1.0 ./ 5.0)
X_lag_logit = hcat(rep_vector(1.0, num_elements(subject_idx)))
pop_lag_logit ~ _popefs_normal(; X = X_lag_logit, beta_loc = [0.0], beta_scale = [2.0])
r_lag_logit_tlag_bsv_subject = b_tlag_bsv_subject[subject_idx, 1]
lag_logit = pop_lag_logit + r_lag_logit_tlag_bsv_subject
X_log_ka = hcat(rep_vector(1.0, num_elements(subject_idx)))
pop_log_ka ~ _popefs_normal(; X = X_log_ka, beta_loc = [0.0], beta_scale = [0.35364652069384966])
r_log_ka_ka_bsv_subject = b_ka_bsv_subject[subject_idx, 1]
log_ka = pop_log_ka + r_log_ka_ka_bsv_subject
X_log_cl0 = hcat(rep_vector(1.0, num_elements(log_weight_ratio)))
pop_log_cl0 ~ _popefs_normal(; X = X_log_cl0, beta_loc = [-2.3025850929940455], beta_scale = [1.1747883127520642])
r_log_cl0_cl_bsv_subject = b_cl_bsv_subject[subject_idx, 1]
log_cl0 = pop_log_cl0 + 0.75 .* log_weight_ratio + r_log_cl0_cl_bsv_subject
X_log_v0 = hcat(rep_vector(1.0, num_elements(log_weight_ratio)))
pop_log_v0 ~ _popefs_normal(; X = X_log_v0, beta_loc = [2.302585092994046], beta_scale = [1.1747883127520642])
r_log_v0_v_bsv_subject = b_v_bsv_subject[subject_idx, 1]
log_v0 = pop_log_v0 + log_weight_ratio + r_log_v0_v_bsv_subject
X_log_r0 = hcat(rep_vector(1.0, num_elements(subject_idx)))
pop_log_r0 ~ _popefs_normal(; X = X_log_r0, beta_loc = [4.382026634673881], beta_scale = [1.1747883127520642])
r_log_r0_r0_bsv_subject = b_r0_bsv_subject[subject_idx, 1]
log_r0 = pop_log_r0 + r_log_r0_r0_bsv_subject
X_log_inv_kout = hcat(rep_vector(1.0, num_elements(subject_idx)))
pop_log_inv_kout ~ _popefs_normal(; X = X_log_inv_kout, beta_loc = [3.4011973816621555], beta_scale = [1.1747883127520642])
r_log_inv_kout_kout_bsv_subject = b_kout_bsv_subject[subject_idx, 1]
log_inv_kout = pop_log_inv_kout + r_log_inv_kout_kout_bsv_subject
X_log_ec50 = hcat(rep_vector(1.0, num_elements(subject_idx)))
pop_log_ec50 ~ _popefs_normal(; X = X_log_ec50, beta_loc = [0.9162907318741551], beta_scale = [1.1747883127520642])
r_log_ec50_ec50_bsv_subject = b_ec50_bsv_subject[subject_idx, 1]
log_ec50 = pop_log_ec50 + r_log_ec50_ec50_bsv_subject
pk_pred ~ plate(pk_time, dose, pk_dv, lag_logit, log_ka, log_cl0, log_v0; outer = (kernel_nsub_pk_pred,)) do times, dose_i, observed, lag_i, lka_i, lcl_i, lv_i
#= brm-docs-example.jl:38 =#
log_tlag_i = log_inv_logit(lag_i)
#= brm-docs-example.jl:39 =#
prediction = exp(warfarin_pk_logconcentration(times, log(dose_i), lka_i, lcl_i, lv_i, log_tlag_i)) + 1.0e-5
#= brm-docs-example.jl:42 =#
pk_pointwise_loglik = warfarin_gamma2_overdisp_lpdfs(observed, prediction, sigma_pk, kappa_pk * 25.0)
#= brm-docs-example.jl:45 =#
observed ~ warfarin_gamma2_overdisp(prediction, sigma_pk, kappa_pk * 25.0)
#= brm-docs-example.jl:48 =#
prediction
end
pd_pred ~ plate(pd_time, dose, pd_dv, lag_logit, log_ka, log_cl0, log_v0, log_r0, log_inv_kout, log_ec50; outer = (kernel_nsub_pd_pred,)) do times, dose_i, observed, lag_i, lka_i, lcl_i, lv_i, lr0_i, linvkout_i, lec50_i
#= brm-docs-example.jl:58 =#
log_tlag_i = log_inv_logit(lag_i)
#= brm-docs-example.jl:59 =#
prediction = warfarin_turnover_prediction(times, log(dose_i), lka_i, lcl_i, lv_i, log_tlag_i, lr0_i, linvkout_i, lec50_i)
#= brm-docs-example.jl:63 =#
pd_pointwise_loglik = warfarin_gamma2_overdisp_lpdfs(observed, prediction, sigma_pd, kappa_pd * 625.0)
#= brm-docs-example.jl:66 =#
observed ~ warfarin_gamma2_overdisp(prediction, sigma_pd, kappa_pd * 625.0)
#= brm-docs-example.jl:69 =#
prediction
end
endstan
functions {
matrix hcat(vector x) {
int n = dims(x)[1];
return to_matrix(x, n, 1);
}
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 warfarin_pk_logconcentration(
vector tad,
real log_dose,
real log_ka,
real log_cl,
real log_v,
real log_tlag
) {
int n = dims(tad)[1];
real log_ke = (log_cl - log_v);
real log_delta = warfarin_log_diff_exp_abs(log_ka, log_ke);
real log_scale = (((log_dose - log_v) + log_ka) - log_delta);
vector[n] value;
for(i in 1:n) {
if((tad[i] <= exp(log_tlag))) {
value[i] = negative_infinity();
} else {
real log_tad = log((tad[i] - exp(log_tlag)));
real a = (-exp((log_ke + log_tad)));
real b = (-exp((log_ka + log_tad)));
value[i] = (log_scale + warfarin_log_diff_exp_abs(a, b));
}
}
return value;
}
real warfarin_log_diff_exp_abs(
real log_a,
real log_b
) {
return (0.5 * log_diff_exp(log_sum_exp((2.0 * log_a), (2.0 * log_b)), (log(2.0) + log_a + log_b)));
}
vector warfarin_gamma2_overdisp_lpdfs(
vector y,
vector mu,
real sigma,
real kappa
) {
int n = dims(y)[1];
if (dims(mu)[1] != n) reject("warfarin_gamma2_overdisp_lpdfs: dim mismatch — `mu` dim 1 (= ", dims(mu)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `mu` dim 1 (= ", dims(mu)[1], ").");
vector[n] lp;
for(i in 1:n) {
real variance = (square(sigma) + (square(mu[i]) / kappa));
real shape = (square(mu[i]) / variance);
real rate = (mu[i] / variance);
lp[i] = gamma_lpdf(y[i] | shape, rate);
}
return lp;
}
vector warfarin_gamma2_overdisp_vector_rng(
int anontok__1,
vector mu,
real sigma,
real kappa
) {
int n = anontok__1;
if (dims(mu)[1] != n) reject("warfarin_gamma2_overdisp_rng: dim mismatch — `mu` dim 1 (= ", dims(mu)[1], ") does not match `n` (= ", n, "), inferred from `anontok__1` dim 1. `n` sizes: `anontok__1` dim 1 (= ", anontok__1, "), `mu` dim 1 (= ", dims(mu)[1], ").");
vector[n] draw;
for(i in 1:n) {
draw[i] = warfarin_gamma2_overdisp_rng(mu[i], sigma, kappa);
}
return draw;
}
real warfarin_gamma2_overdisp_rng(
real mu,
real sigma,
real kappa
) {
real variance = (square(sigma) + (square(mu) / kappa));
real shape = (square(mu) / variance);
real rate = (mu / variance);
return gamma_rng(shape, rate);
}
real warfarin_gamma2_overdisp_lpdf(
vector y,
vector mu,
real sigma,
real kappa
) {
int n = dims(y)[1];
if (dims(mu)[1] != n) reject("warfarin_gamma2_overdisp_lpdf: dim mismatch — `mu` dim 1 (= ", dims(mu)[1], ") does not match `n` (= ", n, "), inferred from `y` dim 1. `n` sizes: `y` dim 1 (= ", dims(y)[1], "), `mu` dim 1 (= ", dims(mu)[1], ").");
real lp = 0.0;
for(i in 1:n) {
real variance = (square(sigma) + (square(mu[i]) / kappa));
real shape = (square(mu[i]) / variance);
real rate = (mu[i] / variance);
lp += gamma_lpdf(y[i] | shape, rate);
}
return lp;
}
vector warfarin_turnover_prediction(
vector times,
real log_dose,
real log_ka,
real log_cl,
real log_v,
real log_tlag,
real log_r0,
real log_inv_kout,
real log_ec50
) {
int n = dims(times)[1];
vector[1] initial_state = rep_vector(exp(log_r0), 1);
array[dims(times)[1]] vector[1] trajectory = ode_rk45_tol(
warfarin_turnover_rhs,
initial_state,
-0.0001,
to_array_1d(times),
1.0e-5,
0.001,
500,
log_dose,
log_ka,
log_cl,
log_v,
log_tlag,
log_r0,
log_inv_kout,
log_ec50
);
return to_vector(trajectory[:, 1]);
}
vector warfarin_turnover_rhs(
real time,
vector state,
real log_dose,
real log_ka,
real log_cl,
real log_v,
real log_tlag,
real log_r0,
real log_inv_kout,
real log_ec50
) {
int m = dims(state)[1];
real log_conc = warfarin_pk_logconcentration_one(time, log_dose, log_ka, log_cl, log_v, log_tlag);
real log_kout = (-log_inv_kout);
real log_kin = (log_r0 + log_kout);
real log_inhibition = log_inv_logit((log_conc - log_ec50));
vector[m] derivative;
derivative[1] = (exp((log_kin + log1m_exp(log_inhibition))) - (state[1] * exp(log_kout)));
return derivative;
}
real warfarin_pk_logconcentration_one(
real time,
real log_dose,
real log_ka,
real log_cl,
real log_v,
real log_tlag
) {
if((time < exp(log_tlag))) {
return -25.0;
}
real log_ke = (log_cl - log_v);
real log_delta = warfarin_log_diff_exp_abs(log_ka, log_ke);
real log_scale = (((log_dose - log_v) + log_ka) - log_delta);
real log_tad = log((time - exp(log_tlag)));
real a = (-exp((log_ke + log_tad)));
real b = (-exp((log_ka + log_tad)));
return (log_scale + warfarin_log_diff_exp_abs(a, b));
}
}
data {
int n_terms_tlag_bsv_subject;
int n_subject;
int n_terms_ka_bsv_subject;
int n_terms_cl_bsv_subject;
int n_terms_v_bsv_subject;
int n_terms_r0_bsv_subject;
int n_terms_kout_bsv_subject;
int n_terms_ec50_bsv_subject;
int subject_idx_n;
array[subject_idx_n] int subject_idx;
int log_weight_ratio_n;
vector[log_weight_ratio_n] log_weight_ratio;
int kernel_nsub_pk_pred;
int pk_time_ends_n;
int pk_time_mem_n;
tuple(vector[pk_time_mem_n], array[pk_time_ends_n] int) pk_time;
int pk_dv_mem_n;
int pk_dv_ends_n;
tuple(vector[pk_dv_mem_n], array[pk_dv_ends_n] int) pk_dv;
int dose_n;
vector[dose_n] dose;
int kernel_nsub_pd_pred;
int pd_time_ends_n;
int pd_time_mem_n;
tuple(vector[pd_time_mem_n], array[pd_time_ends_n] int) pd_time;
int pd_dv_mem_n;
int pd_dv_ends_n;
tuple(vector[pd_dv_mem_n], array[pd_dv_ends_n] int) pd_dv;
}
transformed data {
matrix[num_elements(subject_idx), 1] X_lag_logit = hcat(rep_vector(1.0, num_elements(subject_idx)));
int pop_lag_logit_n_covariates = 1;
matrix[num_elements(subject_idx), 1] X_log_ka = hcat(rep_vector(1.0, num_elements(subject_idx)));
int pop_log_ka_n_covariates = 1;
matrix[num_elements(log_weight_ratio), 1] X_log_cl0 = hcat(rep_vector(1.0, num_elements(log_weight_ratio)));
int pop_log_cl0_n_covariates = 1;
matrix[num_elements(log_weight_ratio), 1] X_log_v0 = hcat(rep_vector(1.0, num_elements(log_weight_ratio)));
int pop_log_v0_n_covariates = 1;
matrix[num_elements(subject_idx), 1] X_log_r0 = hcat(rep_vector(1.0, num_elements(subject_idx)));
int pop_log_r0_n_covariates = 1;
matrix[num_elements(subject_idx), 1] X_log_inv_kout = hcat(rep_vector(1.0, num_elements(subject_idx)));
int pop_log_inv_kout_n_covariates = 1;
matrix[num_elements(subject_idx), 1] X_log_ec50 = hcat(rep_vector(1.0, num_elements(subject_idx)));
int pop_log_ec50_n_covariates = 1;
array[kernel_nsub_pk_pred] int pk_pred_prediction__pl_len_1;
array[kernel_nsub_pk_pred] int pk_pred_pk_pointwise_loglik__pl_len_1;
array[kernel_nsub_pk_pred] int pk_pred__pl_len_1;
for(plate_i__pl_1 in 1:kernel_nsub_pk_pred) {
pk_pred_prediction__pl_len_1[plate_i__pl_1] = (1 + (ragged_end(pk_time.2, plate_i__pl_1) - ragged_start(pk_time.2, plate_i__pl_1)));
pk_pred_pk_pointwise_loglik__pl_len_1[plate_i__pl_1] = (1 + (ragged_end(pk_time.2, plate_i__pl_1) - ragged_start(pk_time.2, plate_i__pl_1)));
pk_pred__pl_len_1[plate_i__pl_1] = (1 + (ragged_end(pk_time.2, plate_i__pl_1) - ragged_start(pk_time.2, plate_i__pl_1)));
}
array[kernel_nsub_pk_pred] int pk_pred_prediction__pl_end_1 = cumulative_sum(pk_pred_prediction__pl_len_1);
array[kernel_nsub_pk_pred] int pk_pred_pk_pointwise_loglik__pl_end_1 = cumulative_sum(pk_pred_pk_pointwise_loglik__pl_len_1);
array[kernel_nsub_pk_pred] int pk_pred__pl_end_1 = cumulative_sum(pk_pred__pl_len_1);
array[kernel_nsub_pd_pred] int pd_pred__pl_len_2;
array[kernel_nsub_pd_pred] int pd_pred_pd_pointwise_loglik__pl_len_2;
array[kernel_nsub_pd_pred] int pd_pred_prediction__pl_len_2;
for(plate_i__pl_2 in 1:kernel_nsub_pd_pred) {
pd_pred__pl_len_2[plate_i__pl_2] = (1 + (ragged_end(pd_time.2, plate_i__pl_2) - ragged_start(pd_time.2, plate_i__pl_2)));
pd_pred_pd_pointwise_loglik__pl_len_2[plate_i__pl_2] = (1 + (ragged_end(pd_time.2, plate_i__pl_2) - ragged_start(pd_time.2, plate_i__pl_2)));
pd_pred_prediction__pl_len_2[plate_i__pl_2] = (1 + (ragged_end(pd_time.2, plate_i__pl_2) - ragged_start(pd_time.2, plate_i__pl_2)));
}
array[kernel_nsub_pd_pred] int pd_pred__pl_end_2 = cumulative_sum(pd_pred__pl_len_2);
array[kernel_nsub_pd_pred] int pd_pred_pd_pointwise_loglik__pl_end_2 = cumulative_sum(pd_pred_pd_pointwise_loglik__pl_len_2);
array[kernel_nsub_pd_pred] int pd_pred_prediction__pl_end_2 = cumulative_sum(pd_pred_prediction__pl_len_2);
}
parameters {
cholesky_factor_corr[n_terms_tlag_bsv_subject] b_tlag_bsv_subject_L;
vector<lower=0.0>[n_terms_tlag_bsv_subject] b_tlag_bsv_subject_tau;
vector[(n_terms_tlag_bsv_subject * n_subject)] b_tlag_bsv_subject_z_flat;
cholesky_factor_corr[n_terms_ka_bsv_subject] b_ka_bsv_subject_L;
vector<lower=0.0>[n_terms_ka_bsv_subject] b_ka_bsv_subject_tau;
vector[(n_terms_ka_bsv_subject * n_subject)] b_ka_bsv_subject_z_flat;
cholesky_factor_corr[n_terms_cl_bsv_subject] b_cl_bsv_subject_L;
vector<lower=0.0>[n_terms_cl_bsv_subject] b_cl_bsv_subject_tau;
vector[(n_terms_cl_bsv_subject * n_subject)] b_cl_bsv_subject_z_flat;
cholesky_factor_corr[n_terms_v_bsv_subject] b_v_bsv_subject_L;
vector<lower=0.0>[n_terms_v_bsv_subject] b_v_bsv_subject_tau;
vector[(n_terms_v_bsv_subject * n_subject)] b_v_bsv_subject_z_flat;
cholesky_factor_corr[n_terms_r0_bsv_subject] b_r0_bsv_subject_L;
vector<lower=0.0>[n_terms_r0_bsv_subject] b_r0_bsv_subject_tau;
vector[(n_terms_r0_bsv_subject * n_subject)] b_r0_bsv_subject_z_flat;
cholesky_factor_corr[n_terms_kout_bsv_subject] b_kout_bsv_subject_L;
vector<lower=0.0>[n_terms_kout_bsv_subject] b_kout_bsv_subject_tau;
vector[(n_terms_kout_bsv_subject * n_subject)] b_kout_bsv_subject_z_flat;
cholesky_factor_corr[n_terms_ec50_bsv_subject] b_ec50_bsv_subject_L;
vector<lower=0.0>[n_terms_ec50_bsv_subject] b_ec50_bsv_subject_tau;
vector[(n_terms_ec50_bsv_subject * n_subject)] b_ec50_bsv_subject_z_flat;
real<lower=0.0> sigma_pk;
real<lower=0.0> kappa_pk;
real<lower=0.0> sigma_pd;
real<lower=0.0> kappa_pd;
vector[pop_lag_logit_n_covariates] pop_lag_logit_beta_pop;
vector[pop_log_ka_n_covariates] pop_log_ka_beta_pop;
vector[pop_log_cl0_n_covariates] pop_log_cl0_beta_pop;
vector[pop_log_v0_n_covariates] pop_log_v0_beta_pop;
vector[pop_log_r0_n_covariates] pop_log_r0_beta_pop;
vector[pop_log_inv_kout_n_covariates] pop_log_inv_kout_beta_pop;
vector[pop_log_ec50_n_covariates] pop_log_ec50_beta_pop;
}
transformed parameters {
matrix[n_terms_tlag_bsv_subject, n_subject] b_tlag_bsv_subject_z = to_matrix(b_tlag_bsv_subject_z_flat, n_terms_tlag_bsv_subject, n_subject);
matrix[n_subject, n_terms_tlag_bsv_subject] b_tlag_bsv_subject = ((diag_pre_multiply(b_tlag_bsv_subject_tau, b_tlag_bsv_subject_L) * b_tlag_bsv_subject_z)');
matrix[n_terms_ka_bsv_subject, n_subject] b_ka_bsv_subject_z = to_matrix(b_ka_bsv_subject_z_flat, n_terms_ka_bsv_subject, n_subject);
matrix[n_subject, n_terms_ka_bsv_subject] b_ka_bsv_subject = ((diag_pre_multiply(b_ka_bsv_subject_tau, b_ka_bsv_subject_L) * b_ka_bsv_subject_z)');
matrix[n_terms_cl_bsv_subject, n_subject] b_cl_bsv_subject_z = to_matrix(b_cl_bsv_subject_z_flat, n_terms_cl_bsv_subject, n_subject);
matrix[n_subject, n_terms_cl_bsv_subject] b_cl_bsv_subject = ((diag_pre_multiply(b_cl_bsv_subject_tau, b_cl_bsv_subject_L) * b_cl_bsv_subject_z)');
matrix[n_terms_v_bsv_subject, n_subject] b_v_bsv_subject_z = to_matrix(b_v_bsv_subject_z_flat, n_terms_v_bsv_subject, n_subject);
matrix[n_subject, n_terms_v_bsv_subject] b_v_bsv_subject = ((diag_pre_multiply(b_v_bsv_subject_tau, b_v_bsv_subject_L) * b_v_bsv_subject_z)');
matrix[n_terms_r0_bsv_subject, n_subject] b_r0_bsv_subject_z = to_matrix(b_r0_bsv_subject_z_flat, n_terms_r0_bsv_subject, n_subject);
matrix[n_subject, n_terms_r0_bsv_subject] b_r0_bsv_subject = ((diag_pre_multiply(b_r0_bsv_subject_tau, b_r0_bsv_subject_L) * b_r0_bsv_subject_z)');
matrix[n_terms_kout_bsv_subject, n_subject] b_kout_bsv_subject_z = to_matrix(b_kout_bsv_subject_z_flat, n_terms_kout_bsv_subject, n_subject);
matrix[n_subject, n_terms_kout_bsv_subject] b_kout_bsv_subject = ((diag_pre_multiply(b_kout_bsv_subject_tau, b_kout_bsv_subject_L) * b_kout_bsv_subject_z)');
matrix[n_terms_ec50_bsv_subject, n_subject] b_ec50_bsv_subject_z = to_matrix(b_ec50_bsv_subject_z_flat, n_terms_ec50_bsv_subject, n_subject);
matrix[n_subject, n_terms_ec50_bsv_subject] b_ec50_bsv_subject = ((diag_pre_multiply(b_ec50_bsv_subject_tau, b_ec50_bsv_subject_L) * b_ec50_bsv_subject_z)');
vector[num_elements(subject_idx)] pop_lag_logit = (X_lag_logit * pop_lag_logit_beta_pop);
vector[subject_idx_n] r_lag_logit_tlag_bsv_subject = b_tlag_bsv_subject[subject_idx, 1];
vector[num_elements(subject_idx)] lag_logit = (pop_lag_logit + r_lag_logit_tlag_bsv_subject);
vector[num_elements(subject_idx)] pop_log_ka = (X_log_ka * pop_log_ka_beta_pop);
vector[subject_idx_n] r_log_ka_ka_bsv_subject = b_ka_bsv_subject[subject_idx, 1];
vector[num_elements(subject_idx)] log_ka = (pop_log_ka + r_log_ka_ka_bsv_subject);
vector[num_elements(log_weight_ratio)] pop_log_cl0 = (X_log_cl0 * pop_log_cl0_beta_pop);
vector[subject_idx_n] r_log_cl0_cl_bsv_subject = b_cl_bsv_subject[subject_idx, 1];
vector[num_elements(log_weight_ratio)] log_cl0 = (pop_log_cl0 + (0.75 .* log_weight_ratio) + r_log_cl0_cl_bsv_subject);
vector[num_elements(log_weight_ratio)] pop_log_v0 = (X_log_v0 * pop_log_v0_beta_pop);
vector[subject_idx_n] r_log_v0_v_bsv_subject = b_v_bsv_subject[subject_idx, 1];
vector[num_elements(log_weight_ratio)] log_v0 = (pop_log_v0 + log_weight_ratio + r_log_v0_v_bsv_subject);
vector[num_elements(subject_idx)] pop_log_r0 = (X_log_r0 * pop_log_r0_beta_pop);
vector[subject_idx_n] r_log_r0_r0_bsv_subject = b_r0_bsv_subject[subject_idx, 1];
vector[num_elements(subject_idx)] log_r0 = (pop_log_r0 + r_log_r0_r0_bsv_subject);
vector[num_elements(subject_idx)] pop_log_inv_kout = (X_log_inv_kout * pop_log_inv_kout_beta_pop);
vector[subject_idx_n] r_log_inv_kout_kout_bsv_subject = b_kout_bsv_subject[subject_idx, 1];
vector[num_elements(subject_idx)] log_inv_kout = (pop_log_inv_kout + r_log_inv_kout_kout_bsv_subject);
vector[num_elements(subject_idx)] pop_log_ec50 = (X_log_ec50 * pop_log_ec50_beta_pop);
vector[subject_idx_n] r_log_ec50_ec50_bsv_subject = b_ec50_bsv_subject[subject_idx, 1];
vector[num_elements(subject_idx)] log_ec50 = (pop_log_ec50 + r_log_ec50_ec50_bsv_subject);
real pk_pred__pl_inv1_1 = (kappa_pk * 25.0);
vector[sum(pk_pred_prediction__pl_len_1)] pk_pred_prediction__pl_mem_1;
vector[kernel_nsub_pk_pred] pk_pred_log_tlag_i;
for(plate_i__pl_1 in 1:kernel_nsub_pk_pred) {
pk_pred_log_tlag_i[plate_i__pl_1] = log_inv_logit(lag_logit[plate_i__pl_1]);
pk_pred_prediction__pl_mem_1[
ragged_start(pk_pred_prediction__pl_end_1, plate_i__pl_1):ragged_end(pk_pred_prediction__pl_end_1, plate_i__pl_1)
] = (
exp(
warfarin_pk_logconcentration(
pk_time.1[ragged_start(pk_time.2, plate_i__pl_1):ragged_end(pk_time.2, plate_i__pl_1)],
log(dose[plate_i__pl_1]),
log_ka[plate_i__pl_1],
log_cl0[plate_i__pl_1],
log_v0[plate_i__pl_1],
pk_pred_log_tlag_i[plate_i__pl_1]
)
) +
1.0e-5
);
}
real pd_pred__pl_inv1_2 = (kappa_pd * 625.0);
vector[sum(pd_pred_prediction__pl_len_2)] pd_pred_prediction__pl_mem_2;
vector[kernel_nsub_pd_pred] pd_pred_log_tlag_i;
for(plate_i__pl_2 in 1:kernel_nsub_pd_pred) {
pd_pred_log_tlag_i[plate_i__pl_2] = log_inv_logit(lag_logit[plate_i__pl_2]);
pd_pred_prediction__pl_mem_2[
ragged_start(pd_pred_prediction__pl_end_2, plate_i__pl_2):ragged_end(pd_pred_prediction__pl_end_2, plate_i__pl_2)
] = warfarin_turnover_prediction(
pd_time.1[ragged_start(pd_time.2, plate_i__pl_2):ragged_end(pd_time.2, plate_i__pl_2)],
log(dose[plate_i__pl_2]),
log_ka[plate_i__pl_2],
log_cl0[plate_i__pl_2],
log_v0[plate_i__pl_2],
pd_pred_log_tlag_i[plate_i__pl_2],
log_r0[plate_i__pl_2],
log_inv_kout[plate_i__pl_2],
log_ec50[plate_i__pl_2]
);
}
}
model {
b_tlag_bsv_subject_L ~ lkj_corr_cholesky(1.0);
b_tlag_bsv_subject_tau ~ normal(0.0, 0.5);
b_tlag_bsv_subject_z_flat ~ std_normal();
b_ka_bsv_subject_L ~ lkj_corr_cholesky(1.0);
b_ka_bsv_subject_tau ~ normal(0.0, 0.5);
b_ka_bsv_subject_z_flat ~ std_normal();
b_cl_bsv_subject_L ~ lkj_corr_cholesky(1.0);
b_cl_bsv_subject_tau ~ normal(0.0, 0.5);
b_cl_bsv_subject_z_flat ~ std_normal();
b_v_bsv_subject_L ~ lkj_corr_cholesky(1.0);
b_v_bsv_subject_tau ~ normal(0.0, 0.5);
b_v_bsv_subject_z_flat ~ std_normal();
b_r0_bsv_subject_L ~ lkj_corr_cholesky(1.0);
b_r0_bsv_subject_tau ~ normal(0.0, 0.5);
b_r0_bsv_subject_z_flat ~ std_normal();
b_kout_bsv_subject_L ~ lkj_corr_cholesky(1.0);
b_kout_bsv_subject_tau ~ normal(0.0, 0.5);
b_kout_bsv_subject_z_flat ~ std_normal();
b_ec50_bsv_subject_L ~ lkj_corr_cholesky(1.0);
b_ec50_bsv_subject_tau ~ normal(0.0, 0.5);
b_ec50_bsv_subject_z_flat ~ std_normal();
sigma_pk ~ normal(0.0, 2.0);
kappa_pk ~ gamma(0.2, (1.0 ./ 5.0));
sigma_pd ~ normal(0.0, 10.0);
kappa_pd ~ gamma(0.2, (1.0 ./ 5.0));
pop_lag_logit_beta_pop ~ normal([0.0]', [2.0]');
pop_log_ka_beta_pop ~ normal([0.0]', [0.35364652069384966]');
pop_log_cl0_beta_pop ~ normal([-2.3025850929940455]', [1.1747883127520642]');
pop_log_v0_beta_pop ~ normal([2.302585092994046]', [1.1747883127520642]');
pop_log_r0_beta_pop ~ normal([4.382026634673881]', [1.1747883127520642]');
pop_log_inv_kout_beta_pop ~ normal([3.4011973816621555]', [1.1747883127520642]');
pop_log_ec50_beta_pop ~ normal([0.9162907318741551]', [1.1747883127520642]');
for(plate_i__pl_1 in 1:kernel_nsub_pk_pred) {
pk_dv.1[ragged_start(pk_dv.2, plate_i__pl_1):ragged_end(pk_dv.2, plate_i__pl_1)] ~ warfarin_gamma2_overdisp(
pk_pred_prediction__pl_mem_1[
ragged_start(pk_pred_prediction__pl_end_1, plate_i__pl_1):ragged_end(pk_pred_prediction__pl_end_1, plate_i__pl_1)
],
sigma_pk,
pk_pred__pl_inv1_1
);
}
for(plate_i__pl_2 in 1:kernel_nsub_pd_pred) {
pd_dv.1[ragged_start(pd_dv.2, plate_i__pl_2):ragged_end(pd_dv.2, plate_i__pl_2)] ~ warfarin_gamma2_overdisp(
pd_pred_prediction__pl_mem_2[
ragged_start(pd_pred_prediction__pl_end_2, plate_i__pl_2):ragged_end(pd_pred_prediction__pl_end_2, plate_i__pl_2)
],
sigma_pd,
pd_pred__pl_inv1_2
);
}
}
generated quantities {
vector[sum(pk_pred_pk_pointwise_loglik__pl_len_1)] pk_pred_pk_pointwise_loglik__pl_mem_1;
vector[sum(pk_pred__pl_len_1)] pk_pred__pl_mem_1;
vector[num_elements(pk_dv.1)] pk_dv_gen;
vector[num_elements(pk_dv.2)] pk_dv_likelihood;
for(plate_i__pl_1 in 1:kernel_nsub_pk_pred) {
pk_pred_pk_pointwise_loglik__pl_mem_1[
ragged_start(pk_pred_pk_pointwise_loglik__pl_end_1, plate_i__pl_1):ragged_end(pk_pred_pk_pointwise_loglik__pl_end_1, plate_i__pl_1)
] = warfarin_gamma2_overdisp_lpdfs(
pk_dv.1[ragged_start(pk_dv.2, plate_i__pl_1):ragged_end(pk_dv.2, plate_i__pl_1)],
pk_pred_prediction__pl_mem_1[
ragged_start(pk_pred_prediction__pl_end_1, plate_i__pl_1):ragged_end(pk_pred_prediction__pl_end_1, plate_i__pl_1)
],
sigma_pk,
pk_pred__pl_inv1_1
);
pk_dv_gen[ragged_start(pk_dv.2, plate_i__pl_1):ragged_end(pk_dv.2, plate_i__pl_1)] = warfarin_gamma2_overdisp_vector_rng(
(1 + (ragged_end(pk_dv.2, plate_i__pl_1) - ragged_start(pk_dv.2, plate_i__pl_1))),
pk_pred_prediction__pl_mem_1[
ragged_start(pk_pred_prediction__pl_end_1, plate_i__pl_1):ragged_end(pk_pred_prediction__pl_end_1, plate_i__pl_1)
],
sigma_pk,
pk_pred__pl_inv1_1
);
pk_dv_likelihood[plate_i__pl_1] = warfarin_gamma2_overdisp_lpdf(pk_dv.1[ragged_start(pk_dv.2, plate_i__pl_1):ragged_end(pk_dv.2, plate_i__pl_1)] |
pk_pred_prediction__pl_mem_1[
ragged_start(pk_pred_prediction__pl_end_1, plate_i__pl_1):ragged_end(pk_pred_prediction__pl_end_1, plate_i__pl_1)
],
sigma_pk,
pk_pred__pl_inv1_1
);
pk_pred__pl_mem_1[
ragged_start(pk_pred__pl_end_1, plate_i__pl_1):ragged_end(pk_pred__pl_end_1, plate_i__pl_1)
] = pk_pred_prediction__pl_mem_1[
ragged_start(pk_pred_prediction__pl_end_1, plate_i__pl_1):ragged_end(pk_pred_prediction__pl_end_1, plate_i__pl_1)
];
}
vector[sum(pd_pred__pl_len_2)] pd_pred__pl_mem_2;
vector[sum(pd_pred_pd_pointwise_loglik__pl_len_2)] pd_pred_pd_pointwise_loglik__pl_mem_2;
vector[num_elements(pd_dv.1)] pd_dv_gen;
vector[num_elements(pd_dv.2)] pd_dv_likelihood;
for(plate_i__pl_2 in 1:kernel_nsub_pd_pred) {
pd_pred_pd_pointwise_loglik__pl_mem_2[
ragged_start(pd_pred_pd_pointwise_loglik__pl_end_2, plate_i__pl_2):ragged_end(pd_pred_pd_pointwise_loglik__pl_end_2, plate_i__pl_2)
] = warfarin_gamma2_overdisp_lpdfs(
pd_dv.1[ragged_start(pd_dv.2, plate_i__pl_2):ragged_end(pd_dv.2, plate_i__pl_2)],
pd_pred_prediction__pl_mem_2[
ragged_start(pd_pred_prediction__pl_end_2, plate_i__pl_2):ragged_end(pd_pred_prediction__pl_end_2, plate_i__pl_2)
],
sigma_pd,
pd_pred__pl_inv1_2
);
pd_dv_gen[ragged_start(pd_dv.2, plate_i__pl_2):ragged_end(pd_dv.2, plate_i__pl_2)] = warfarin_gamma2_overdisp_vector_rng(
(1 + (ragged_end(pd_dv.2, plate_i__pl_2) - ragged_start(pd_dv.2, plate_i__pl_2))),
pd_pred_prediction__pl_mem_2[
ragged_start(pd_pred_prediction__pl_end_2, plate_i__pl_2):ragged_end(pd_pred_prediction__pl_end_2, plate_i__pl_2)
],
sigma_pd,
pd_pred__pl_inv1_2
);
pd_dv_likelihood[plate_i__pl_2] = warfarin_gamma2_overdisp_lpdf(pd_dv.1[ragged_start(pd_dv.2, plate_i__pl_2):ragged_end(pd_dv.2, plate_i__pl_2)] |
pd_pred_prediction__pl_mem_2[
ragged_start(pd_pred_prediction__pl_end_2, plate_i__pl_2):ragged_end(pd_pred_prediction__pl_end_2, plate_i__pl_2)
],
sigma_pd,
pd_pred__pl_inv1_2
);
pd_pred__pl_mem_2[
ragged_start(pd_pred__pl_end_2, plate_i__pl_2):ragged_end(pd_pred__pl_end_2, plate_i__pl_2)
] = pd_pred_prediction__pl_mem_2[
ragged_start(pd_pred_prediction__pl_end_2, plate_i__pl_2):ragged_end(pd_pred_prediction__pl_end_2, plate_i__pl_2)
];
}
}julia
Turing unsupported for this BRM example
Turing backend: direct execution requires at least one observed likelihoodThe 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 boundary
The public brms issue mentions a Warfarin PK/PD model but supplies no equations, source, data, or citation. It therefore cannot establish identity with an unavailable private model.
The public two-stage sections faithfully reproduce Weber's StanCon 2018 Warfarin model directory, including its separate PK and PD programs. “Faithful” refers to that public two-stage program; it is not a claim about unspecified private source. The joint section is explicitly a new model built from those public likelihoods, not a claim that the public program was joint.
Run the complete reproductions
After bootstrapping the repository's test environment, run the public two-stage checks with:
sh
BRM_WARFARIN_RUNTIME=1 julia --startup-file=no --project=test \
research/warfarin/reproduce.jlRun the joint checks with:
sh
BRM_WARFARIN_JOINT_RUNTIME=1 julia --startup-file=no --project=test \
research/warfarin/joint_reproduce.jlSet the corresponding runtime variable to 0 to run lowering and stanc without BridgeStan instantiation.