One-Compartment IV Infusion
Here, we write down the PK model for a one-compartment model with IV infusion and linear elimination. We also write down the statistical model for multiple error models. The Stan models are free to download using the corresponding download buttons.
1 PK Model
The data-generating model is \[C_{ij} = f(\mathbf{\theta}_i, t_{ij}) + error\] where \(error\) can be
- additive: \(C_{ij} = f(\mathbf{\theta}_i, t_{ij}) + \epsilon_{a_{ij}}\)
- proportional: \(C_{ij} = f(\mathbf{\theta}_i, t_{ij})*\left( 1 + \epsilon_{p_{ij}} \right)\)
- proportional-plus-additive: \(C_{ij} = f(\mathbf{\theta}_i, t_{ij})*\left( 1 + \epsilon_{p_{ij}} \right) + \epsilon_{a_{ij}}\)
- lognormal: \(C_{ij} = f(\mathbf{\theta}_i, t_{ij})*e^{\epsilon_{ij}}\)
where \(\mathbf{\theta}_i = \left[CL_i, \, V_{c_i}\right]^\top\) is a vector containing the individual parameters for individual \(i\)1, \(f(\cdot, \cdot)\) is the one-compartment model mean function seen here, and \(\epsilon_p, \epsilon_a,\) and \(\epsilon\) are normally distributed. The corresponding system of ODEs is:
\[\begin{align} \frac{dA_c}{dt} &= -\left(\frac{CL}{V_c}\right)A_C + rate_{in} \notag \\ \end{align} \tag{1}\]
and then \(C = \frac{A_c}{V_c}\)2.
2 Statistical Model
I’ll model the correlation matrix \(R_{\Omega}\) and standard deviations (\(\omega_p\)) of the random effects rather than the covariance matrix \(\Omega\) that is the more traditional method.
\[\begin{align} \Omega &= \begin{pmatrix} \omega^2_{CL} & \omega_{CL, V_c} \\ \omega_{CL, V_c} & \omega^2_{V_c} \\ \end{pmatrix} \\ &= \begin{pmatrix} \omega_{CL} & 0 \\ 0 & \omega_{V_c} \\ \end{pmatrix} \mathbf{R_{\Omega}} \begin{pmatrix} \omega_{CL} & 0 \\ 0 & \omega_{V_c} \\ \end{pmatrix} \\ &= \begin{pmatrix} \omega_{CL} & 0 \\ 0 & \omega_{V_c} \\ \end{pmatrix} \begin{pmatrix} 1 & \rho_{CL, V_c} \\ \rho_{CL, V_c} & 1 \\ \end{pmatrix} \begin{pmatrix} \omega_{CL} & 0 \\ 0 & \omega_{V_c} \\ \end{pmatrix} \end{align}\]
When the random effects are modeled this way, we get all of the same information as if they were modeled with a full covariance matrix, i.e., we can transform the standard deviations and correlation matrix into the covariance matrix if we want. In addition, standard deviations and correlations are more interpretable than variances and covariances. It’s also much easier to set a prior distribution on the standard deviations and the correlation matrix with an LKJ distribution than it is to set a prior on a covariance matrix with an inverse-Wishart distribution.
The full statistical model is
\[\begin{align} C_{ij} \mid \mathbf{TV}, \; \mathbf{\eta}_i, \; \mathbf{\Omega}, \; \sigma_p &\sim Normal\left( f(\mathbf{\theta}_i, t_{ij}), \; \sigma_{ij} \right) I(C_{ij} > 0) \notag \\ \mathbf{\eta}_i \; | \; \Omega &\sim Normal\left( \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} , \; \Omega\right) \notag \\ TVCL &\sim Lognormal\left(log\left(location_{TVCL}\right), scale_{TVCL}\right) \notag \\ TVVC &\sim Lognormal\left(log\left(location_{TVVC}\right), scale_{TVVC}\right) \notag \\ \omega_{CL} &\sim Half-Normal(0, scale_{\omega_{CL}}) \notag \\ \omega_{V_c} &\sim Half-Normal(0, scale_{\omega_{V_c}}) \notag \\ R_{\Omega} &\sim LKJ(df_{R_{\Omega}}) \notag \\ \sigma_p &\sim Half-Normal(0, scale_{\sigma_p}) \\ \end{align}\]
where \[\begin{align} \mathbf{TV} &= \begin{pmatrix} TVCL \\ TVVC \\ \end{pmatrix} \\ \mathbf{\theta}_i &= \begin{pmatrix} CL_i \\ V_{c_i} \\ \end{pmatrix} = \begin{pmatrix} TVCL \times e^{\eta_{CL_i}} \\ TVVC \times e^{\eta_{V_{c_i}}} \\ \end{pmatrix} \\ \mathbf{\eta}_i &= \begin{pmatrix} \eta_{CL_i} \\ \eta_{V_{c_i}} \\ \end{pmatrix} \\ \sigma_{ij} &= f(\mathbf{\theta}_i, t_{ij})\sigma_p \\ \end{align}\]
Note: The indicator for \(C_{ij} | \ldots\) indicates that we are truncating the distribution of the observed concentrations to be greater than 0.
\[\begin{align} C_{ij} \mid \mathbf{TV}, \; \mathbf{\eta}_i, \; \mathbf{\Omega}, \; \mathbf{\Sigma} &\sim Normal\left( f(\mathbf{\theta}_i, t_{ij}), \; \sigma_{ij} \right) I(C_{ij} > 0) \notag \\ \mathbf{\eta}_i \; | \; \Omega &\sim Normal\left( \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} , \; \Omega\right) \notag \\ TVCL &\sim Lognormal\left(log\left(location_{TVCL}\right), scale_{TVCL}\right) \notag \\ TVVC &\sim Lognormal\left(log\left(location_{TVVC}\right), scale_{TVVC}\right) \notag \\ \omega_{CL} &\sim Half-Normal(0, scale_{\omega_{CL}}) \notag \\ \omega_{V_c} &\sim Half-Normal(0, scale_{\omega_{V_c}}) \notag \\ R_{\Omega} &\sim LKJ(df_{R_{\Omega}}) \notag \\ \sigma_p &\sim Half-Normal(0, scale_{\sigma_p}) \\ \sigma_a &\sim Half-Normal(0, scale_{\sigma_a}) \\ R_{\Sigma} &\sim LKJ(df_{R_{\Sigma}}) \notag \\ \end{align}\]
where \[\begin{align} \mathbf{TV} &= \begin{pmatrix} TVCL \\ TVVC \\ \end{pmatrix} \\ \mathbf{\theta}_i &= \begin{pmatrix} CL_i \\ V_{c_i} \\ \end{pmatrix} = \begin{pmatrix} TVCL \times e^{\eta_{CL_i}} \\ TVVC \times e^{\eta_{V_{c_i}}} \\ \end{pmatrix} \\ \mathbf{\eta}_i &= \begin{pmatrix} \eta_{CL_i} \\ \eta_{V_{c_i}} \\ \end{pmatrix} \\ \mathbf{\Sigma} &= \begin{pmatrix} \sigma^2_{p} & \rho_{p,a}\sigma_{p}\sigma_{a} \\ \rho_{p,a}\sigma_{p}\sigma_{a} & \sigma^2_{z} \\ \end{pmatrix} \\ \sigma_{ij} &= \sqrt{f(\mathbf{\theta}_i, t_{ij})^2\sigma^2_p + \sigma^2_a + 2f(\mathbf{\theta}_i, t_{ij})\rho_{p,a}\sigma_{p}\sigma_{a}} \\ \end{align}\]
Note: The indicator for \(C_{ij} | \ldots\) indicates that we are truncating the distribution of the observed concentrations to be greater than 0.
\[\begin{align} C_{ij} \mid \mathbf{TV}, \; \mathbf{\eta}_i, \; \mathbf{\Omega}, \; \sigma &\sim Lognormal\left( log\left(f(\mathbf{\theta}_i, t_{ij})\right), \; \sigma \right) \notag \\ \mathbf{\eta}_i \; | \; \Omega &\sim Normal\left( \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} , \; \Omega\right) \notag \\ TVCL &\sim Lognormal\left(log\left(location_{TVCL}\right), scale_{TVCL}\right) \notag \\ TVVC &\sim Lognormal\left(log\left(location_{TVVC}\right), scale_{TVVC}\right) \notag \\ \omega_{CL} &\sim Half-Normal(0, scale_{\omega_{CL}}) \notag \\ \omega_{V_c} &\sim Half-Normal(0, scale_{\omega_{V_c}}) \notag \\ R_{\Omega} &\sim LKJ(df_{R_{\Omega}}) \notag \\ \sigma &\sim Half-Normal(0, scale_{\sigma}) \\ \end{align}\]
where \[\begin{align} \mathbf{TV} &= \begin{pmatrix} TVCL \\ TVVC \\ \end{pmatrix} \\ \mathbf{\theta}_i &= \begin{pmatrix} CL_i \\ V_{c_i} \\ \end{pmatrix} = \begin{pmatrix} TVCL \times e^{\eta_{CL_i}} \\ TVVC \times e^{\eta_{V_{c_i}}} \\ \end{pmatrix} \\ \mathbf{\eta}_i &= \begin{pmatrix} \eta_{CL_i} \\ \eta_{V_{c_i}} \\ \end{pmatrix} \\ \end{align}\]
3 Stan Code
3.1 Proportional Error
// IV infusion
// One-compartment PK Model
// IIV on CL and VC (full covariance matrix)
// proportional error - DV = IPRED*(1 + eps_p)
// Matrix-exponential solution using Torsten
// Implements threading for within-chain parallelization
// Deals with BLOQ values by the "CDF trick" (M4)
// Since we have a normal distribution on the error, but the DV must be > 0, it
// truncates the likelihood below at 0
// For PPC, it generates values from a normal that is truncated below at 0
functions{
array[] int sequence(int start, int end) {
array[end - start + 1] int seq;
for (n in 1:num_elements(seq)) {
seq[n] = n + start - 1;
}
return seq;
}
int num_between(int lb, int ub, array[] int y){
int n = 0;
for(i in 1:num_elements(y)){
if(y[i] >= lb && y[i] <= ub)
n = n + 1;
}
return n;
}
array[] int find_between(int lb, int ub, array[] int y) {
// vector[num_between(lb, ub, y)] result;
array[num_between(lb, ub, y)] int result;
int n = 1;
for (i in 1:num_elements(y)) {
if (y[i] >= lb && y[i] <= ub) {
result[n] = y[i];
n = n + 1;
}
}
return result;
}
vector find_between_vec(int lb, int ub, array[] int idx, vector y) {
vector[num_between(lb, ub, idx)] result;
int n = 1;
if(num_elements(idx) != num_elements(y)) reject("illegal input");
for (i in 1:rows(y)) {
if (idx[i] >= lb && idx[i] <= ub) {
result[n] = y[i];
n = n + 1;
}
}
return result;
}
real normal_lb_rng(real mu, real sigma, real lb){
real p_lb = normal_cdf(lb | mu, sigma);
real u = uniform_rng(p_lb, 1);
real y = mu + sigma * inv_Phi(u);
return y;
}
real partial_sum_lpmf(array[] int seq_subj, int start, int end,
vector dv_obs, array[] int dv_obs_id, array[] int i_obs,
array[] real amt, array[] int cmt, array[] int evid,
array[] real time, array[] real rate, array[] real ii,
array[] int addl, array[] int ss,
array[] int subj_start, array[] int subj_end,
vector CL, vector VC,
real sigma_p,
vector lloq, array[] int bloq,
int n_random, int n_subjects, int n_total,
array[] real bioav, array[] real tlag, int n_cmt){
real ptarget = 0;
int N = end - start + 1; // number of subjects in this slice
vector[n_total] dv_ipred;
matrix[n_total, n_cmt] x_ipred;
int n_obs_slice = num_between(subj_start[start], subj_end[end], i_obs);
array[n_obs_slice] int i_obs_slice = find_between(subj_start[start],
subj_end[end], i_obs);
vector[n_obs_slice] dv_obs_slice = find_between_vec(start, end,
dv_obs_id, dv_obs);
vector[n_obs_slice] ipred_slice;
vector[n_obs_slice] lloq_slice = lloq[i_obs_slice];
array[n_obs_slice] int bloq_slice = bloq[i_obs_slice];
for(n in 1:N){ // loop over subjects in this slice
int j = n + start - 1; // j is the ID of the current subject
real ke = CL[j]/VC[j];
matrix[n_cmt, n_cmt] K = rep_matrix(0, n_cmt, n_cmt);
K[1, 1] = -ke;
x_ipred[subj_start[j]:subj_end[j], ] =
pmx_solve_linode(time[subj_start[j]:subj_end[j]],
amt[subj_start[j]:subj_end[j]],
rate[subj_start[j]:subj_end[j]],
ii[subj_start[j]:subj_end[j]],
evid[subj_start[j]:subj_end[j]],
cmt[subj_start[j]:subj_end[j]],
addl[subj_start[j]:subj_end[j]],
ss[subj_start[j]:subj_end[j]],
K, bioav, tlag)';
dv_ipred[subj_start[j]:subj_end[j]] =
x_ipred[subj_start[j]:subj_end[j], 1] ./ VC[j];
}
ipred_slice = dv_ipred[i_obs_slice];
for(i in 1:n_obs_slice){
real sigma_tmp = ipred_slice[i]*sigma_p;
if(bloq_slice[i] == 1){
ptarget += log_diff_exp(normal_lcdf(lloq_slice[i] | ipred_slice[i],
sigma_tmp),
normal_lcdf(0.0 | ipred_slice[i], sigma_tmp)) -
normal_lccdf(0.0 | ipred_slice[i], sigma_tmp);
}else{
ptarget += normal_lpdf(dv_obs_slice[i] | ipred_slice[i], sigma_tmp) -
normal_lccdf(0.0 | ipred_slice[i], sigma_tmp);
}
}
return ptarget;
}
}
data{
int n_subjects;
int n_total;
int n_obs;
array[n_obs] int i_obs;
array[n_total] int ID;
array[n_total] real amt;
array[n_total] int cmt;
array[n_total] int evid;
array[n_total] real rate;
array[n_total] real ii;
array[n_total] int addl;
array[n_total] int ss;
array[n_total] real time;
vector<lower = 0>[n_total] dv;
array[n_subjects] int subj_start;
array[n_subjects] int subj_end;
vector[n_total] lloq;
array[n_total] int bloq;
real<lower = 0> location_tvcl; // Prior Location parameter for CL
real<lower = 0> location_tvvc; // Prior Location parameter for VC
real<lower = 0> scale_tvcl; // Prior Scale parameter for CL
real<lower = 0> scale_tvvc; // Prior Scale parameter for VC
real<lower = 0> scale_omega_cl; // Prior scale parameter for omega_cl
real<lower = 0> scale_omega_vc; // Prior scale parameter for omega_vc
real<lower = 0> lkj_df_omega; // Prior degrees of freedom for omega cor mat
real<lower = 0> scale_sigma_p; // Prior Scale parameter for proportional error
int<lower = 0, upper = 1> prior_only; // Want to simulate from the prior?
}
transformed data{
int grainsize = 1;
vector<lower = 0>[n_obs] dv_obs = dv[i_obs];
array[n_obs] int dv_obs_id = ID[i_obs];
vector[n_obs] lloq_obs = lloq[i_obs];
array[n_obs] int bloq_obs = bloq[i_obs];
int n_random = 2; // Number of random effects
int n_cmt = 1; // Number of states in the ODEs
array[n_random] real scale_omega = {scale_omega_cl, scale_omega_vc};
array[n_subjects] int seq_subj = sequence(1, n_subjects); // reduce_sum over subjects
array[n_cmt] real bioav = rep_array(1.0, n_cmt); // Hardcoding, but could be data or a parameter in another situation
array[n_cmt] real tlag = rep_array(0.0, n_cmt);
}
parameters{
real<lower = 0> TVCL;
real<lower = 0> TVVC;
vector<lower = 0>[n_random] omega;
cholesky_factor_corr[n_random] L;
real<lower = 0> sigma_p;
matrix[n_random, n_subjects] Z;
}
transformed parameters{
vector[n_subjects] eta_cl;
vector[n_subjects] eta_vc;
vector[n_subjects] CL;
vector[n_subjects] VC;
vector[n_subjects] KE;
{
row_vector[n_random] typical_values = to_row_vector({TVCL, TVVC});
matrix[n_subjects, n_random] eta = diag_pre_multiply(omega, L * Z)';
matrix[n_subjects, n_random] theta =
(rep_matrix(typical_values, n_subjects) .* exp(eta));
eta_cl = col(eta, 1);
eta_vc = col(eta, 2);
CL = col(theta, 1);
VC = col(theta, 2);
KE = CL ./ VC;
}
}
model{
// Priors
TVCL ~ lognormal(log(location_tvcl), scale_tvcl);
TVVC ~ lognormal(log(location_tvvc), scale_tvvc);
omega ~ normal(0, scale_omega);
L ~ lkj_corr_cholesky(lkj_df_omega);
sigma_p ~ normal(0, scale_sigma_p);
to_vector(Z) ~ std_normal();
// Likelihood
if(prior_only == 0){
target += reduce_sum(partial_sum_lupmf, seq_subj, grainsize,
dv_obs, dv_obs_id, i_obs,
amt, cmt, evid, time,
rate, ii, addl, ss, subj_start, subj_end,
CL, VC,
sigma_p,
lloq, bloq,
n_random, n_subjects, n_total,
bioav, tlag, n_cmt);
}
}
generated quantities{
real<lower = 0> sigma_sq_p = square(sigma_p);
real<lower = 0> omega_cl = omega[1];
real<lower = 0> omega_vc = omega[2];
real<lower = 0> omega_sq_cl = square(omega_cl);
real<lower = 0> omega_sq_vc = square(omega_vc);
real cor_cl_vc;
real omega_cl_vc;
vector[n_obs] ipred;
vector[n_obs] pred;
vector[n_obs] dv_ppc;
vector[n_obs] log_lik;
vector[n_obs] res;
vector[n_obs] wres;
vector[n_obs] ires;
vector[n_obs] iwres;
{
matrix[n_random, n_random] R = multiply_lower_tri_self_transpose(L);
matrix[n_random, n_random] Omega = quad_form_diag(R, omega);
vector[n_total] dv_pred;
matrix[n_total, n_cmt] x_pred;
vector[n_total] dv_ipred;
matrix[n_total, n_cmt] x_ipred;
cor_cl_vc = R[1, 2];
omega_cl_vc = Omega[1, 2];
for(j in 1:n_subjects){
real tvke = TVCL/TVVC;
matrix[n_cmt, n_cmt] K = rep_matrix(0, n_cmt, n_cmt);
matrix[n_cmt, n_cmt] K_tv = rep_matrix(0, n_cmt, n_cmt);
K[1, 1] = -KE[j];
x_ipred[subj_start[j]:subj_end[j], ] =
pmx_solve_linode(time[subj_start[j]:subj_end[j]],
amt[subj_start[j]:subj_end[j]],
rate[subj_start[j]:subj_end[j]],
ii[subj_start[j]:subj_end[j]],
evid[subj_start[j]:subj_end[j]],
cmt[subj_start[j]:subj_end[j]],
addl[subj_start[j]:subj_end[j]],
ss[subj_start[j]:subj_end[j]],
K, bioav, tlag)';
K_tv[1, 1] = -tvke;
x_pred[subj_start[j]:subj_end[j],] =
pmx_solve_linode(time[subj_start[j]:subj_end[j]],
amt[subj_start[j]:subj_end[j]],
rate[subj_start[j]:subj_end[j]],
ii[subj_start[j]:subj_end[j]],
evid[subj_start[j]:subj_end[j]],
cmt[subj_start[j]:subj_end[j]],
addl[subj_start[j]:subj_end[j]],
ss[subj_start[j]:subj_end[j]],
K_tv, bioav, tlag)';
dv_ipred[subj_start[j]:subj_end[j]] =
x_ipred[subj_start[j]:subj_end[j], 1] ./ VC[j];
dv_pred[subj_start[j]:subj_end[j]] =
x_pred[subj_start[j]:subj_end[j], 1] ./ TVVC;
}
pred = dv_pred[i_obs];
ipred = dv_ipred[i_obs];
}
res = dv_obs - pred;
ires = dv_obs - ipred;
for(i in 1:n_obs){
real ipred_tmp = ipred[i];
real sigma_tmp = ipred_tmp*sigma_p;
dv_ppc[i] = normal_lb_rng(ipred_tmp, sigma_tmp, 0.0);
if(bloq_obs[i] == 1){
// log_lik[i] = log(normal_cdf(lloq_obs[i] | ipred_tmp, sigma_tmp) -
// normal_cdf(0.0 | ipred_tmp, sigma_tmp)) -
// normal_lccdf(0.0 | ipred_tmp, sigma_tmp);
log_lik[i] = log_diff_exp(normal_lcdf(lloq_obs[i] | ipred_tmp, sigma_tmp),
normal_lcdf(0.0 | ipred_tmp, sigma_tmp)) -
normal_lccdf(0.0 | ipred_tmp, sigma_tmp);
}else{
log_lik[i] = normal_lpdf(dv_obs[i] | ipred_tmp, sigma_tmp) -
normal_lccdf(0.0 | ipred_tmp, sigma_tmp);
}
wres[i] = res[i]/sigma_tmp;
iwres[i] = ires[i]/sigma_tmp;
}
}3.2 Proportional-Plus-Additive Error
// IV infusion
// One-compartment PK Model
// IIV on CL and VC (full covariance matrix)
// proportional plus additive error - DV = IPRED*(1 + eps_p) + eps_a
// Matrix-exponential solution using Torsten
// Implements threading for within-chain parallelization
// Deals with BLOQ values by the "CDF trick" (M4)
// Since we have a normal distribution on the error, but the DV must be > 0, it
// truncates the likelihood below at 0
// For PPC, it generates values from a normal that is truncated below at 0
functions{
array[] int sequence(int start, int end) {
array[end - start + 1] int seq;
for (n in 1:num_elements(seq)) {
seq[n] = n + start - 1;
}
return seq;
}
int num_between(int lb, int ub, array[] int y){
int n = 0;
for(i in 1:num_elements(y)){
if(y[i] >= lb && y[i] <= ub)
n = n + 1;
}
return n;
}
array[] int find_between(int lb, int ub, array[] int y) {
// vector[num_between(lb, ub, y)] result;
array[num_between(lb, ub, y)] int result;
int n = 1;
for (i in 1:num_elements(y)) {
if (y[i] >= lb && y[i] <= ub) {
result[n] = y[i];
n = n + 1;
}
}
return result;
}
vector find_between_vec(int lb, int ub, array[] int idx, vector y) {
vector[num_between(lb, ub, idx)] result;
int n = 1;
if(num_elements(idx) != num_elements(y)) reject("illegal input");
for (i in 1:rows(y)) {
if (idx[i] >= lb && idx[i] <= ub) {
result[n] = y[i];
n = n + 1;
}
}
return result;
}
real normal_lb_rng(real mu, real sigma, real lb){
real p_lb = normal_cdf(lb | mu, sigma);
real u = uniform_rng(p_lb, 1);
real y = mu + sigma * inv_Phi(u);
return y;
}
real partial_sum_lpmf(array[] int seq_subj, int start, int end,
vector dv_obs, array[] int dv_obs_id, array[] int i_obs,
array[] real amt, array[] int cmt, array[] int evid,
array[] real time, array[] real rate, array[] real ii,
array[] int addl, array[] int ss,
array[] int subj_start, array[] int subj_end,
vector CL, vector VC,
real sigma_sq_p, real sigma_sq_a, real sigma_p_a,
vector lloq, array[] int bloq,
int n_random, int n_subjects, int n_total,
array[] real bioav, array[] real tlag, int n_cmt){
real ptarget = 0;
int N = end - start + 1; // number of subjects in this slice
vector[n_total] dv_ipred;
matrix[n_total, n_cmt] x_ipred;
int n_obs_slice = num_between(subj_start[start], subj_end[end], i_obs);
array[n_obs_slice] int i_obs_slice = find_between(subj_start[start],
subj_end[end], i_obs);
vector[n_obs_slice] dv_obs_slice = find_between_vec(start, end,
dv_obs_id, dv_obs);
vector[n_obs_slice] ipred_slice;
vector[n_obs_slice] lloq_slice = lloq[i_obs_slice];
array[n_obs_slice] int bloq_slice = bloq[i_obs_slice];
for(n in 1:N){ // loop over subjects in this slice
int j = n + start - 1; // j is the ID of the current subject
real ke = CL[j]/VC[j];
matrix[n_cmt, n_cmt] K = rep_matrix(0, n_cmt, n_cmt);
K[1, 1] = -ke;
x_ipred[subj_start[j]:subj_end[j], ] =
pmx_solve_linode(time[subj_start[j]:subj_end[j]],
amt[subj_start[j]:subj_end[j]],
rate[subj_start[j]:subj_end[j]],
ii[subj_start[j]:subj_end[j]],
evid[subj_start[j]:subj_end[j]],
cmt[subj_start[j]:subj_end[j]],
addl[subj_start[j]:subj_end[j]],
ss[subj_start[j]:subj_end[j]],
K, bioav, tlag)';
dv_ipred[subj_start[j]:subj_end[j]] =
x_ipred[subj_start[j]:subj_end[j], 1] ./ VC[j];
}
ipred_slice = dv_ipred[i_obs_slice];
for(i in 1:n_obs_slice){
real ipred_tmp = ipred_slice[i];
real sigma_tmp = sqrt(square(ipred_tmp) * sigma_sq_p + sigma_sq_a +
2*ipred_tmp*sigma_p_a);
if(bloq_slice[i] == 1){
ptarget += log_diff_exp(normal_lcdf(lloq_slice[i] | ipred_tmp,
sigma_tmp),
normal_lcdf(0.0 | ipred_tmp, sigma_tmp)) -
normal_lccdf(0.0 | ipred_tmp, sigma_tmp);
}else{
ptarget += normal_lpdf(dv_obs_slice[i] | ipred_tmp, sigma_tmp) -
normal_lccdf(0.0 | ipred_tmp, sigma_tmp);
}
}
return ptarget;
}
}
data{
int n_subjects;
int n_total;
int n_obs;
array[n_obs] int i_obs;
array[n_total] int ID;
array[n_total] real amt;
array[n_total] int cmt;
array[n_total] int evid;
array[n_total] real rate;
array[n_total] real ii;
array[n_total] int addl;
array[n_total] int ss;
array[n_total] real time;
vector<lower = 0>[n_total] dv;
array[n_subjects] int subj_start;
array[n_subjects] int subj_end;
vector[n_total] lloq;
array[n_total] int bloq;
real<lower = 0> location_tvcl; // Prior Location parameter for CL
real<lower = 0> location_tvvc; // Prior Location parameter for VC
real<lower = 0> scale_tvcl; // Prior Scale parameter for CL
real<lower = 0> scale_tvvc; // Prior Scale parameter for VC
real<lower = 0> scale_omega_cl; // Prior scale parameter for omega_cl
real<lower = 0> scale_omega_vc; // Prior scale parameter for omega_vc
real<lower = 0> lkj_df_omega; // Prior degrees of freedom for omega cor mat
real<lower = 0> scale_sigma_p; // Prior Scale parameter for proportional error
real<lower = 0> scale_sigma_a; // Prior Scale parameter for additive error
real<lower = 0> lkj_df_sigma; // Prior degrees of freedom for sigma cor mat
int<lower = 0, upper = 1> prior_only; // Want to simulate from the prior?
}
transformed data{
int grainsize = 1;
vector<lower = 0>[n_obs] dv_obs = dv[i_obs];
array[n_obs] int dv_obs_id = ID[i_obs];
vector[n_obs] lloq_obs = lloq[i_obs];
array[n_obs] int bloq_obs = bloq[i_obs];
int n_random = 2; // Number of random effects
int n_cmt = 1; // Number of states in the ODEs
array[n_random] real scale_omega = {scale_omega_cl, scale_omega_vc};
array[2] real scale_sigma = {scale_sigma_p, scale_sigma_a};
array[n_subjects] int seq_subj = sequence(1, n_subjects); // reduce_sum over subjects
array[n_cmt] real bioav = rep_array(1.0, n_cmt); // Hardcoding, but could be data or a parameter in another situation
array[n_cmt] real tlag = rep_array(0.0, n_cmt);
}
parameters{
real<lower = 0> TVCL;
real<lower = 0> TVVC;
vector<lower = 0>[n_random] omega;
cholesky_factor_corr[n_random] L;
vector<lower = 0>[2] sigma;
cholesky_factor_corr[2] L_Sigma;
matrix[n_random, n_subjects] Z;
}
transformed parameters{
vector[n_subjects] eta_cl;
vector[n_subjects] eta_vc;
vector[n_subjects] CL;
vector[n_subjects] VC;
vector[n_subjects] KE;
real<lower = 0> sigma_p = sigma[1];
real<lower = 0> sigma_a = sigma[2];
real<lower = 0> sigma_sq_p = square(sigma_p);
real<lower = 0> sigma_sq_a = square(sigma_a);
real cor_p_a;
real sigma_p_a;
{
row_vector[n_random] typical_values = to_row_vector({TVCL, TVVC});
matrix[n_subjects, n_random] eta = diag_pre_multiply(omega, L * Z)';
matrix[n_subjects, n_random] theta =
(rep_matrix(typical_values, n_subjects) .* exp(eta));
matrix[2, 2] R_Sigma = multiply_lower_tri_self_transpose(L_Sigma);
matrix[2, 2] Sigma = quad_form_diag(R_Sigma, sigma);
eta_cl = col(eta, 1);
eta_vc = col(eta, 2);
CL = col(theta, 1);
VC = col(theta, 2);
KE = CL ./ VC;
cor_p_a = R_Sigma[1, 2];
sigma_p_a = Sigma[1, 2];
}
}
model{
// Priors
TVCL ~ lognormal(log(location_tvcl), scale_tvcl);
TVVC ~ lognormal(log(location_tvvc), scale_tvvc);
omega ~ normal(0, scale_omega);
L ~ lkj_corr_cholesky(lkj_df_omega);
sigma ~ normal(0, scale_sigma);
L_Sigma ~ lkj_corr_cholesky(lkj_df_sigma);
to_vector(Z) ~ std_normal();
// Likelihood
if(prior_only == 0){
target += reduce_sum(partial_sum_lupmf, seq_subj, grainsize,
dv_obs, dv_obs_id, i_obs,
amt, cmt, evid, time,
rate, ii, addl, ss, subj_start, subj_end,
CL, VC,
sigma_sq_p, sigma_sq_a, sigma_p_a,
lloq, bloq,
n_random, n_subjects, n_total,
bioav, tlag, n_cmt);
}
}
generated quantities{
real<lower = 0> omega_cl = omega[1];
real<lower = 0> omega_vc = omega[2];
real<lower = 0> omega_sq_cl = square(omega_cl);
real<lower = 0> omega_sq_vc = square(omega_vc);
real cor_cl_vc;
real omega_cl_vc;
vector[n_obs] ipred;
vector[n_obs] pred;
vector[n_obs] dv_ppc;
vector[n_obs] log_lik;
vector[n_obs] res;
vector[n_obs] wres;
vector[n_obs] ires;
vector[n_obs] iwres;
{
matrix[n_random, n_random] R = multiply_lower_tri_self_transpose(L);
matrix[n_random, n_random] Omega = quad_form_diag(R, omega);
vector[n_total] dv_pred;
matrix[n_total, n_cmt] x_pred;
vector[n_total] dv_ipred;
matrix[n_total, n_cmt] x_ipred;
cor_cl_vc = R[1, 2];
omega_cl_vc = Omega[1, 2];
for(j in 1:n_subjects){
real tvke = TVCL/TVVC;
matrix[n_cmt, n_cmt] K = rep_matrix(0, n_cmt, n_cmt);
matrix[n_cmt, n_cmt] K_tv = rep_matrix(0, n_cmt, n_cmt);
K[1, 1] = -KE[j];
x_ipred[subj_start[j]:subj_end[j], ] =
pmx_solve_linode(time[subj_start[j]:subj_end[j]],
amt[subj_start[j]:subj_end[j]],
rate[subj_start[j]:subj_end[j]],
ii[subj_start[j]:subj_end[j]],
evid[subj_start[j]:subj_end[j]],
cmt[subj_start[j]:subj_end[j]],
addl[subj_start[j]:subj_end[j]],
ss[subj_start[j]:subj_end[j]],
K, bioav, tlag)';
K_tv[1, 1] = -tvke;
x_pred[subj_start[j]:subj_end[j],] =
pmx_solve_linode(time[subj_start[j]:subj_end[j]],
amt[subj_start[j]:subj_end[j]],
rate[subj_start[j]:subj_end[j]],
ii[subj_start[j]:subj_end[j]],
evid[subj_start[j]:subj_end[j]],
cmt[subj_start[j]:subj_end[j]],
addl[subj_start[j]:subj_end[j]],
ss[subj_start[j]:subj_end[j]],
K_tv, bioav, tlag)';
dv_ipred[subj_start[j]:subj_end[j]] =
x_ipred[subj_start[j]:subj_end[j], 1] ./ VC[j];
dv_pred[subj_start[j]:subj_end[j]] =
x_pred[subj_start[j]:subj_end[j], 1] ./ TVVC;
}
pred = dv_pred[i_obs];
ipred = dv_ipred[i_obs];
}
res = dv_obs - pred;
ires = dv_obs - ipred;
for(i in 1:n_obs){
real ipred_tmp = ipred[i];
real sigma_tmp = sqrt(square(ipred_tmp) * sigma_sq_p + sigma_sq_a +
2*ipred_tmp*sigma_p_a);
dv_ppc[i] = normal_lb_rng(ipred_tmp, sigma_tmp, 0.0);
if(bloq_obs[i] == 1){
// log_lik[i] = log(normal_cdf(lloq_obs[i] | ipred_tmp, sigma_tmp) -
// normal_cdf(0.0 | ipred_tmp, sigma_tmp)) -
// normal_lccdf(0.0 | ipred_tmp, sigma_tmp);
log_lik[i] = log_diff_exp(normal_lcdf(lloq_obs[i] | ipred_tmp, sigma_tmp),
normal_lcdf(0.0 | ipred_tmp, sigma_tmp)) -
normal_lccdf(0.0 | ipred_tmp, sigma_tmp);
}else{
log_lik[i] = normal_lpdf(dv_obs[i] | ipred_tmp, sigma_tmp) -
normal_lccdf(0.0 | ipred_tmp, sigma_tmp);
}
wres[i] = res[i]/sigma_tmp;
iwres[i] = ires[i]/sigma_tmp;
}
}3.3 Lognormal Error
// IV infusion
// One-compartment PK Model
// IIV on CL and VC (full covariance matrix)
// exponential error - DV = IPRED*exp(eps)
// Matrix-exponential solution using Torsten
// Implements threading for within-chain parallelization
// Deals with BLOQ values by the M3 method (M3 and M4 are equivalent with this
// error model)
functions{
array[] int sequence(int start, int end) {
array[end - start + 1] int seq;
for (n in 1:num_elements(seq)) {
seq[n] = n + start - 1;
}
return seq;
}
int num_between(int lb, int ub, array[] int y){
int n = 0;
for(i in 1:num_elements(y)){
if(y[i] >= lb && y[i] <= ub)
n = n + 1;
}
return n;
}
array[] int find_between(int lb, int ub, array[] int y) {
// vector[num_between(lb, ub, y)] result;
array[num_between(lb, ub, y)] int result;
int n = 1;
for (i in 1:num_elements(y)) {
if (y[i] >= lb && y[i] <= ub) {
result[n] = y[i];
n = n + 1;
}
}
return result;
}
vector find_between_vec(int lb, int ub, array[] int idx, vector y) {
vector[num_between(lb, ub, idx)] result;
int n = 1;
if(num_elements(idx) != num_elements(y)) reject("illegal input");
for (i in 1:rows(y)) {
if (idx[i] >= lb && idx[i] <= ub) {
result[n] = y[i];
n = n + 1;
}
}
return result;
}
real partial_sum_lpmf(array[] int seq_subj, int start, int end,
vector dv_obs, array[] int dv_obs_id, array[] int i_obs,
array[] real amt, array[] int cmt, array[] int evid,
array[] real time, array[] real rate, array[] real ii,
array[] int addl, array[] int ss,
array[] int subj_start, array[] int subj_end,
vector CL, vector VC,
real sigma,
vector lloq, array[] int bloq,
int n_random, int n_subjects, int n_total,
array[] real bioav, array[] real tlag, int n_cmt){
real ptarget = 0;
int N = end - start + 1; // number of subjects in this slice
vector[n_total] dv_ipred;
matrix[n_total, n_cmt] x_ipred;
int n_obs_slice = num_between(subj_start[start], subj_end[end], i_obs);
array[n_obs_slice] int i_obs_slice = find_between(subj_start[start],
subj_end[end], i_obs);
vector[n_obs_slice] dv_obs_slice = find_between_vec(start, end,
dv_obs_id, dv_obs);
vector[n_obs_slice] ipred_slice;
vector[n_obs_slice] lloq_slice = lloq[i_obs_slice];
array[n_obs_slice] int bloq_slice = bloq[i_obs_slice];
for(n in 1:N){ // loop over subjects in this slice
int j = n + start - 1; // j is the ID of the current subject
real ke = CL[j]/VC[j];
matrix[n_cmt, n_cmt] K = rep_matrix(0, n_cmt, n_cmt);
K[1, 1] = -ke;
x_ipred[subj_start[j]:subj_end[j], ] =
pmx_solve_linode(time[subj_start[j]:subj_end[j]],
amt[subj_start[j]:subj_end[j]],
rate[subj_start[j]:subj_end[j]],
ii[subj_start[j]:subj_end[j]],
evid[subj_start[j]:subj_end[j]],
cmt[subj_start[j]:subj_end[j]],
addl[subj_start[j]:subj_end[j]],
ss[subj_start[j]:subj_end[j]],
K, bioav, tlag)';
dv_ipred[subj_start[j]:subj_end[j]] =
x_ipred[subj_start[j]:subj_end[j], 1] ./ VC[j];
}
ipred_slice = dv_ipred[i_obs_slice];
for(i in 1:n_obs_slice){
if(bloq_slice[i] == 1){
ptarget += lognormal_lcdf(lloq_slice[i] | log(ipred_slice[i]), sigma);
}else{
ptarget += lognormal_lpdf(dv_obs_slice[i] | log(ipred_slice[i]), sigma);
}
}
return ptarget;
}
}
data{
int n_subjects;
int n_total;
int n_obs;
array[n_obs] int i_obs;
array[n_total] int ID;
array[n_total] real amt;
array[n_total] int cmt;
array[n_total] int evid;
array[n_total] real rate;
array[n_total] real ii;
array[n_total] int addl;
array[n_total] int ss;
array[n_total] real time;
vector<lower = 0>[n_total] dv;
array[n_subjects] int subj_start;
array[n_subjects] int subj_end;
vector[n_total] lloq;
array[n_total] int bloq;
real<lower = 0> location_tvcl; // Prior Location parameter for CL
real<lower = 0> location_tvvc; // Prior Location parameter for VC
real<lower = 0> scale_tvcl; // Prior Scale parameter for CL
real<lower = 0> scale_tvvc; // Prior Scale parameter for VC
real<lower = 0> scale_omega_cl; // Prior scale parameter for omega_cl
real<lower = 0> scale_omega_vc; // Prior scale parameter for omega_vc
real<lower = 0> lkj_df_omega; // Prior degrees of freedom for omega cor mat
real<lower = 0> scale_sigma; // Prior Scale parameter for exponential error
int<lower = 0, upper = 1> prior_only; // Want to simulate from the prior?
}
transformed data{
int grainsize = 1;
vector<lower = 0>[n_obs] dv_obs = dv[i_obs];
array[n_obs] int dv_obs_id = ID[i_obs];
vector[n_obs] lloq_obs = lloq[i_obs];
array[n_obs] int bloq_obs = bloq[i_obs];
int n_random = 2; // Number of random effects
int n_cmt = 1; // Number of states in the ODEs
array[n_random] real scale_omega = {scale_omega_cl, scale_omega_vc};
array[n_subjects] int seq_subj = sequence(1, n_subjects); // reduce_sum over subjects
array[n_cmt] real bioav = rep_array(1.0, n_cmt); // Hardcoding, but could be data or a parameter in another situation
array[n_cmt] real tlag = rep_array(0.0, n_cmt);
}
parameters{
real<lower = 0> TVCL;
real<lower = 0> TVVC;
vector<lower = 0>[n_random] omega;
cholesky_factor_corr[n_random] L;
real<lower = 0> sigma;
matrix[n_random, n_subjects] Z;
}
transformed parameters{
vector[n_subjects] eta_cl;
vector[n_subjects] eta_vc;
vector[n_subjects] CL;
vector[n_subjects] VC;
vector[n_subjects] KE;
{
row_vector[n_random] typical_values = to_row_vector({TVCL, TVVC});
matrix[n_subjects, n_random] eta = diag_pre_multiply(omega, L * Z)';
matrix[n_subjects, n_random] theta =
(rep_matrix(typical_values, n_subjects) .* exp(eta));
eta_cl = col(eta, 1);
eta_vc = col(eta, 2);
CL = col(theta, 1);
VC = col(theta, 2);
KE = CL ./ VC;
}
}
model{
// Priors
TVCL ~ lognormal(log(location_tvcl), scale_tvcl);
TVVC ~ lognormal(log(location_tvvc), scale_tvvc);
omega ~ normal(0, scale_omega);
L ~ lkj_corr_cholesky(lkj_df_omega);
sigma ~ normal(0, scale_sigma);
to_vector(Z) ~ std_normal();
// Likelihood
if(prior_only == 0){
target += reduce_sum(partial_sum_lupmf, seq_subj, grainsize,
dv_obs, dv_obs_id, i_obs,
amt, cmt, evid, time,
rate, ii, addl, ss, subj_start, subj_end,
CL, VC,
sigma,
lloq, bloq,
n_random, n_subjects, n_total,
bioav, tlag, n_cmt);
}
}
generated quantities{
real<lower = 0> sigma_sq = square(sigma);
real<lower = 0> omega_cl = omega[1];
real<lower = 0> omega_vc = omega[2];
real<lower = 0> omega_sq_cl = square(omega_cl);
real<lower = 0> omega_sq_vc = square(omega_vc);
real cor_cl_vc;
real omega_cl_vc;
vector[n_obs] ipred;
vector[n_obs] pred;
vector[n_obs] dv_ppc;
vector[n_obs] log_lik;
vector[n_obs] res;
vector[n_obs] wres;
vector[n_obs] ires;
vector[n_obs] iwres;
{
matrix[n_random, n_random] R = multiply_lower_tri_self_transpose(L);
matrix[n_random, n_random] Omega = quad_form_diag(R, omega);
vector[n_total] dv_pred;
matrix[n_total, n_cmt] x_pred;
vector[n_total] dv_ipred;
matrix[n_total, n_cmt] x_ipred;
cor_cl_vc = R[1, 2];
omega_cl_vc = Omega[1, 2];
for(j in 1:n_subjects){
real tvke = TVCL/TVVC;
matrix[n_cmt, n_cmt] K = rep_matrix(0, n_cmt, n_cmt);
matrix[n_cmt, n_cmt] K_tv = rep_matrix(0, n_cmt, n_cmt);
K[1, 1] = -KE[j];
x_ipred[subj_start[j]:subj_end[j], ] =
pmx_solve_linode(time[subj_start[j]:subj_end[j]],
amt[subj_start[j]:subj_end[j]],
rate[subj_start[j]:subj_end[j]],
ii[subj_start[j]:subj_end[j]],
evid[subj_start[j]:subj_end[j]],
cmt[subj_start[j]:subj_end[j]],
addl[subj_start[j]:subj_end[j]],
ss[subj_start[j]:subj_end[j]],
K, bioav, tlag)';
K_tv[1, 1] = -tvke;
x_pred[subj_start[j]:subj_end[j],] =
pmx_solve_linode(time[subj_start[j]:subj_end[j]],
amt[subj_start[j]:subj_end[j]],
rate[subj_start[j]:subj_end[j]],
ii[subj_start[j]:subj_end[j]],
evid[subj_start[j]:subj_end[j]],
cmt[subj_start[j]:subj_end[j]],
addl[subj_start[j]:subj_end[j]],
ss[subj_start[j]:subj_end[j]],
K_tv, bioav, tlag)';
dv_ipred[subj_start[j]:subj_end[j]] =
x_ipred[subj_start[j]:subj_end[j], 1] ./ VC[j];
dv_pred[subj_start[j]:subj_end[j]] =
x_pred[subj_start[j]:subj_end[j], 1] ./ TVVC;
}
pred = dv_pred[i_obs];
ipred = dv_ipred[i_obs];
}
res = log(dv_obs) - log(pred);
ires = log(dv_obs) - log(ipred);
for(i in 1:n_obs){
real log_ipred_tmp = log(ipred[i]);
dv_ppc[i] = lognormal_rng(log_ipred_tmp, sigma);
if(bloq_obs[i] == 1){
log_lik[i] = lognormal_lcdf(lloq_obs[i] | log_ipred_tmp, sigma);
}else{
log_lik[i] = lognormal_lpdf(dv_obs[i] | log_ipred_tmp, sigma);
}
wres[i] = res[i]/sigma;
iwres[i] = ires[i]/sigma;
}
}Footnotes
These parameters are the macro-parameterization of this model - clearance and central compartment volume, respectively. I have chosen to use \(V_c\) for central compartment volume rather than \(V_1\) or simply \(V\) to reduce any confusion about the numbering or notation.↩︎
For this model, the matrix-exponential is faster than the analytical solution with Torsten, so I’ve written the Stan models to use the linear ODE solution. Therefore, since this is purely IV, I’ve written the ODE with one compartment. If you want to use Torsten’s analytical solution, you would write the ODEs as \[\begin{align} \frac{dA_d}{dt} &= -K_a*A_d \notag \\ \frac{dA_c}{dt} &= K_a*A_d - \left(\frac{CL}{V_c}\right)A_C + rate_{in} \notag \\ \end{align} \tag{2}\] You’ll notice that this is written as if there is first-order absorption. This is to emphasize that the central compartment is the second compartment in the setup for Torsten’s analytical solver, and you would have to make sure \(CMT = 2\) for the dosing and observation records. Since there is no absorption, we just leave it out of our model and set it to 0 when using Torsten’s analytical solver.↩︎