Multiple studies; All patients in each study receive 1 or 2 or 3 drugs.; Ex: stdy1=drugA,DrugB,stdy2=drugC,DrugD,DrugE. Dosing single at time = 0. Multiple concentration measurements at time =0 to 30 days. Tried 2compartment linear PK model with crossed random effect between subjects and drugs. Model fails with error below. I tried multiple times with different initial values and prior distributions. Prediction of concentration (cHat) is nan. Can you please suggest me possible reason or example code or reference?. Thanks in advance. >
R code:
1. Get unique subjects and drugs
> unique_subjects <- conc %>% distinct(ID) %>% pull(ID)
> n_subjects_stan <- length(unique_subjects)
> unique_drugs <- conc %>% distinct(COMPN) %>% pull(COMPN)
> n_unique_drugs_stan <- length(unique_drugs)
2. Determine the number of drug administrations
> unique_admin <- conc %>% distinct(ID, COMPN, DOSEA, ROUTE, ROUTEN)
> n_drugs_stan <- nrow(unique_admin)
3. Create arrays for drug dose, subject index, route, and drug ID for each administration
> drug_dose_stan <- unique_admin %>% pull(DOSEA)
> subject_index_stan <- unique_admin %>% pull(ID) %>% match(unique_subjects)
> route_stan <- unique_admin %>% pull(ROUTEN) # Assuming ROUTEN is the numerical route (0 or 1)
> drug_id_stan <- unique_admin %>% pull(COMPN) %>% match(unique_drugs)
4. Get time and concentration observations
> time_stan <- conc %>% pull(TIME)
> cObs_stan <- conc %>% pull(DVOR)
> n_obs_stan <- length(time_stan)
5. Determine start and end indices for each drug administration
> indexed_data <- conc %>%
+ group_by(ID, COMPN) %>%
+ mutate(obs_index = row_number()) %>%
+ ungroup()
Merge obs_index into unique_admin
> unique_admin_with_index <- unique_admin %>%
+ left_join(indexed_data, by = c("ID" = "ID", "COMPN" = "COMPN"))
> start_indices <- unique_admin_with_index %>%
+ group_by(ID, COMPN, DOSEA.x, `ROUTE.x`, ROUTEN.x) %>%
+ summarise(start = min(obs_index), .groups = "drop") %>%
+ pull(start)
> end_indices <- unique_admin_with_index %>%
+ group_by(ID, COMPN, DOSEA.x, `ROUTE.x`, ROUTEN.x) %>%
+ summarise(end = max(obs_index), .groups = "drop") %>%
+ pull(end)
> # --- Prepare the data list for Stan ---
> data_list <- list(
+ nSubjects = n_subjects_stan,
+ nObs = n_obs_stan,
+ nDrugs = n_drugs_stan,
+ nUniqueDrugs = n_unique_drugs_stan,
+ drug_dose = drug_dose_stan,
+ start = start_indices,
+ end = end_indices,
+ time = time_stan,
+ cObs = cObs_stan,
+ subject_index = subject_index_stan,
+ route = route_stan,
+ drug_id = drug_id_stan
+ )
> stan_model_code <- "
+ functions{
+ real twoCptModelIV1(real time, real dose, real CL, real Q, real V1, real V2){
+ real k10 = CL / V1;
+ real k12 = Q / V1;
+ real k21 = Q / V2;
+ real ksum = k10 + k12 + k21;
+ real kprod = k10 * k21;
+ real sqrt_term = sqrt(ksum * ksum - 4.0 * kprod);
+ real alpha = (ksum + sqrt_term) / 2.0;
+ real beta = (ksum - sqrt_term) / 2.0;
+ real A = (dose / V1) * (alpha - k21) / (alpha - beta);
+ real B = (dose / V1) * (k21 - beta) / (alpha - beta);
+ return A * exp(-alpha * time) + B * exp(-beta * time);
+ }
+
+ real twoCptModelSC1(real time, real dose, real CL, real Q, real V1, real V2, real ka){
+ real k10;
+ real k12;
+ real k21;
+ real ksum;
+ vector[3] a;
+ vector[3] alpha;
+
+ k10 = CL / V1;
+ k12 = Q / V1;
+ k21 = Q / V2;
+ ksum = k10 + k12 + k21;
+ alpha[1] = (ksum + sqrt(ksum * ksum - 4.0 * k10 * k21))/2.0;
+ alpha[2] = (ksum - sqrt(ksum * ksum - 4.0 * k10 * k21))/2.0;
+ alpha[3] = ka;
+
+ a[1] = ka * (k21 - alpha[1]) / ((ka - alpha[1]) * (alpha[2] - alpha[1]));
+ a[2] = ka * (k21 - alpha[2]) / ((ka - alpha[2]) * (alpha[1] - alpha[2]));
+ a[3] = -(a[1] + a[2]);
+
+ return (dose / V1) * sum(a .* exp(-alpha * time));
+ }
+
+ vector twoCptModel(real[] time, real dose, real CL, real Q, real V1, real V2, real ka, int route){
+ vector[size(time)] conc;
+ if (route == 0) { // IV bolus
+ for(i in 1:size(time))
+ conc[i] = twoCptModelIV1(time[i], dose, CL, Q, V1, V2);
+ } else { // SC
+ for(i in 1:size(time))
+ conc[i] = twoCptModelSC1(time[i], dose, CL, Q, V1, V2, ka);
+ }
+ return conc;
+ }
+ }
+
+ data{
+ int nSubjects;
+ int nObs;
+ int nDrugs; // Total number of drug administrations
+ int nUniqueDrugs; // Number of unique drugs
+ real drug_dose[nDrugs];
+ int start[nDrugs];
+ int end[nDrugs];
+ real time[nObs];
+ vector[nObs] cObs;
+ int subject_index[nDrugs];
+ int route[nDrugs];
+ int drug_id[nDrugs]; // ID of the drug administered (1 to nUniqueDrugs)
+ }
+
+ transformed data{
+ vector[nObs] logCObs = log(cObs);
+ int nRandom = 5;
+ }
+
+ parameters{
+ real<lower = 0> CLHat;
+ real<lower = 0> QHat;
+ real<lower = 0> V1Hat;
+ real<lower = 0> V2Hat;
+ real<lower = 0> kaHat;
+ corr_matrix[nRandom] rho;
+ vector<lower = 0>[nRandom] omega;
+ real<lower = 0> sigma;
+ vector[nRandom] logtheta[nSubjects, nUniqueDrugs]; // Random effects by subject and drug
+ }
+
+ transformed parameters{
+ vector<lower = 0>[nRandom] thetaHat;
+ cov_matrix[nRandom] Omega;
+ real<lower = 0> CL[nDrugs]; // Index by drug administration
+ real<lower = 0> Q[nDrugs];
+ real<lower = 0> V1[nDrugs];
+ real<lower = 0> V2[nDrugs];
+ real<lower = 0> ka[nDrugs];
+ vector<lower = 0>[nObs] cHat;
+
+ thetaHat[1] = CLHat;
+ thetaHat[2] = QHat;
+ thetaHat[3] = V1Hat;
+ thetaHat[4] = V2Hat;
+ thetaHat[5] = kaHat;
+
+ Omega = quad_form_diag(rho, omega);
+
+ for(d in 1:nDrugs){
+ int s = subject_index[d];
+ int drug = drug_id[d];
+ CL[d] = exp(logtheta[s, drug, 1]);
+ Q[d] = exp(logtheta[s, drug, 2]);
+ V1[d] = exp(logtheta[s, drug, 3]);
+ V2[d] = exp(logtheta[s, drug, 4]);
+ ka[d] = exp(logtheta[s, drug, 5]);
+ }
+
+ for(d in 1:nDrugs){
+ int subject_id = subject_index[d];
+ int r = route[d];
+ real current_dose = drug_dose[d];
+ if (r == 0) { // IV bolus
+ cHat[start[d]:end[d]] = twoCptModel(time[start[d]:end[d]],
+ current_dose, CL[d], Q[d], V1[d], V2[d], 1.0, r);
+ } else { // SC
+ cHat[start[d]:end[d]] = twoCptModel(time[start[d]:end[d]],
+ current_dose, CL[d], Q[d], V1[d], V2[d], ka[d], r);
+ }
+ }
+ }
+
+ model{
+ // Priors based on the scale of the initial values
+ CLHat ~ lognormal(log(0.6), 1); // Mean around exp(log(0.6)) ~= 0.6, SD allows for variability
+ QHat ~ lognormal(log(0.3), 1); // Mean around exp(log(0.3)) ~= 0.3
+ V1Hat ~ lognormal(log(7), 1); // Mean around exp(log(7)) ~= 11
+ V2Hat ~ lognormal(log(10), 1); // Mean around exp(log(10)) ~= 22
+ kaHat ~ lognormal(log(0.7), 1); // Mean around exp(log(0.7)) ~= 0.7
+ sigma ~ lognormal(log(0.2), 1); // Prior for the residual standard deviation
+
+ omega ~ cauchy(0, 1);
+ rho ~ lkj_corr(1);
+
+ // Inter-individual and inter-drug variability
+ for(s in 1:nSubjects){
+ for(drug in 1:nUniqueDrugs){
+ logtheta[s, drug] ~ multi_normal(log(thetaHat), Omega);
+ }
+ }
+
+ logCObs ~ normal(log(cHat), sigma);
+ }
+
+ generated quantities{
+ vector[nRandom] logthetaPred[nSubjects, nUniqueDrugs];
+ vector<lower = 0>[nObs] cHatPred;
+ real<lower = 0> cObsCond[nObs];
+ real<lower = 0> cObsPred[nObs];
+ real<lower = 0> CLPred[nDrugs];
+ real<lower = 0> QPred[nDrugs];
+ real<lower = 0> V1Pred[nDrugs];
+ real<lower = 0> V2Pred[nDrugs];
+ real<lower = 0> kaPred[nDrugs];
+
+ for(j in 1:nSubjects){
+ for(drug in 1:nUniqueDrugs){
+ logthetaPred[j, drug] = multi_normal_rng(log(thetaHat), Omega);
+ }
+ }
+
+ for(d in 1:nDrugs){
+ int s = subject_index[d];
+ int drug = drug_id[d];
+ CLPred[d] = exp(logthetaPred[s, drug, 1]);
+ QPred[d] = exp(logthetaPred[s, drug, 2]);
+ V1Pred[d] = exp(logthetaPred[s, drug, 3]);
+ V2Pred[d] = exp(logthetaPred[s, drug, 4]);
+ kaPred[d] = exp(logthetaPred[s, drug, 5]);
+ }
+
+ for(d in 1:nDrugs){
+ int subject_id = subject_index[d];
+ int r = route[d];
+ real current_dose = drug_dose[d];
+ if (r == 0) { // IV bolus
+ cHatPred[start[d]:end[d]] = twoCptModel(time[start[d]:end[d]],
+ current_dose, CLPred[d], QPred[d], V1Pred[d], V2Pred[d], 1.0, r);
+ } else { // SC
+ cHatPred[start[d]:end[d]] = twoCptModel(time[start[d]:end[d]],
+ current_dose, CLPred[d], QPred[d], V1Pred[d], V2Pred[d], kaPred[d], r);
+ }
+ }
+
+ for(i in 1:nObs){
+ cObsCond[i] = exp(normal_rng(log(cHat[i]), sigma));
+ cObsPred[i] = exp(normal_rng(log(cHatPred[i]), sigma));
+ }
+ }
+ "
R code
> init_values <- list(
+ list(CLHat = 1, QHat = 0.6, V1Hat = 10, V2Hat = 15, kaHat = 1, sigma = 0.5),
+ list(CLHat = 0.2, QHat = 0.1, V1Hat = 3, V2Hat = 5, kaHat = 0.3, sigma = 0.1),
+ list(CLHat = 0.6, QHat = 0.3, V1Hat = 5, V2Hat = 8, kaHat = 0.5, sigma = 0.2),
+ list(CLHat = 0.8, QHat = 0.4, V1Hat = 8, V2Hat = 12, kaHat = 0.8, sigma = 0.3)
+ )
Compile and fit the Stan model
> compiled_model <- cmdstan_model(write_stan_file(stan_model_code))
>
> fit <- compiled_model$sample(
+ data = data_list,
+ chains = 4,
+ parallel_chains = 4,
+ iter_warmup = 500,
+ iter_sampling = 1000,
+ init = init_values
+
+ )
Running MCMC with 4 parallel chains...
Chain 1 Rejecting initial value:
Chain 1 Error evaluating the log probability at the initial value.
Chain 1 Exception: validate transformed params: cHat[13] is nan, but must be greater than or equal to 0 (in '/tmp/RtmpEFaOHX/model-40ad6a6916d2.stan' at line 93)
Chain 1 Exception: validate transformed params: cHat[13] is nan, but must be greater than or equal to 0 (in '/tmp/RtmpEFaOHX/model-40ad6a6916d2.stan' at line 93)
[Edit: escaped code, but didn’t single-space it.]