implementing IV infusion rate and covariate as piecewise input #66
Description
Hello Torsten Developers!
Would like to ask if you could provide feedback on implementing IV infusion rate and WT as piece-wise input. I have browsed the internet extensively but was not successful at finding guidance on how to go about it. My starting point was the neutropeniaPopulation.stan from the example models and would like to keep using the ODE solver for this routine. My goal was to convert the code for popPK for a 2 comp IV infusion with weight as piece-wise input but keep on getting the same error output pointing towards line item real rate = rdummy[1]; which probably indicates real weight = rdummy[2]; is also not recognized as it should be.
the error message reads for all chains:
Chain 1 Unrecoverable error evaluating the log probability at the initial value.
Chain 1 Exception: Exception: array[uni, ...] index: accessing element out of range. index 1 out of range; container is empty and cannot be indexed (in '/tmp/RtmpNNQioA/model-16e2132cfc2485.stan', line 16, column 5 to column 27) (in '/tmp/RtmpNNQioA/model-16e2132cfc2485.stan', line 112, column 4 to line 122, column 57)
Warning: Chain 1 finished unexpectedly!
This is the R code am using for the run and data conversion:
get data file
xdata <- read.csv("/home/andras/Torsten/example-models/covariates_IV/nonmem_IV.csv", as.is = TRUE)
xdata<-metrumrg::as.keyed(xdata, key = c("ID", "TIME", "CMT", "EVID"))
class(xdata)
#xdata <- metrumrg::as.best(as.data.frame(xdata)) #this will put out error
#getS3method("as.best", "data.frame")
nt <- nrow(xdata)
start <- (1:nt)[!duplicated(xdata$ID)]
end <- c(start[-1] - 1, nt)
Indices of records containing observed concentrations
iObs <- with(xdata, (1:nrow(xdata))[!is.na(DV) & EVID == 0])
nObs <- length(iObs)
time_sim <- seq(0, 48, by = 1)
create data set
data <- with(xdata,
list(
nSubjects = length(unique(ID)),
nt = nt,
nObsPK = nObs,
iObsPK = iObs,
amt = AMT,
rate = RATE,
ii = II,
cmt = CMT,
evid = EVID,
addl = ADDL,
ss = SS,
start = start,
end = end,
time = TIME,
cObs = as.numeric(DV[iObs]),
weight = WT,
CLHatPrior = 3.0, # Prior median for CL
QHatPrior = 10.0, # Prior median for Q
V1HatPrior = 10.0, # Prior median for V1
V2HatPrior = 200.0, # Prior median for V2
CLHatPriorCV = 0.25, # Prior CV for CL
QHatPriorCV = 0.25, # Prior CV for Q
V1HatPriorCV = 0.25, # Prior CV for V1
V2HatPriorCV = 0.25 # Prior CV for V2
))
create initial estimates
init <- function(){
list(CLHat = exp(rnorm(1, log(3), 0.2)),
QHat = exp(rnorm(1, log(10), 0.2)),
V1Hat = exp(rnorm(1, log(35), 0.2)),
V2Hat = exp(rnorm(1, log(290), 0.2)),
sigma = abs(rnorm(1, 0.1, 0.02)),
# Between-Subject Variability (IIV)
L = diag(4), # Initialize as identity matrix for LKJ prior
omega = abs(rnorm(4, 0.2, 0.05)), # Std dev of random effects
# Standardized Individual Deviations
etaStd = matrix(rnorm(4 * data$nSubjects, 0, 1), nrow = 4, ncol = data$nSubjects))
}
#with(data, stan_rdump(ls(data), file = paste0(modelName, ".data.R")))
#inits <- init()
#with(inits, stan_rdump(ls(inits), file = paste0(modelName, ".init.R")))
str(data)
Run the model in generated quantities mode to simulate data
fit <- model$sample(
data = data,
init=init,
iter_sampling = 1000,
iter_warmup = 500,
chains = 4,
refresh = 1,
fixed_param = F
)
This is how the stan model code is like at the moment:
functions{
vector twoCptNeutModelODE(real t,
vector x,
array[] real parms,
array[] real rdummy,
array[] int idummy){
real CL = parms[1];
real Q = parms[2];
real V1 = parms[3];
real V2 = parms[4];
real k10 = CL / V1;
real k12 = Q / V1;
real k21 = Q / V2;
vector[2] dxdt;
real rate = rdummy[1];
real weight = rdummy[2];
// Central compartment receives drug directly from infusion
dxdt[1] = - (k10 + k12) * x[1] + k21 * x[2] + rate; // Include infusion rate
dxdt[2] = k12 * x[1] - k21 * x[2]; // Peripheral compartment
return dxdt;
}
}
data{
int<lower = 1> nt;
int<lower = 1> nObsPK;
array[nObsPK] int<lower = 1> iObsPK;
array[nt] real<lower = 0> amt;
array[nt] int cmt;
array[nt] int<lower = 0> evid;
array[nt] real<lower = 0> time;
array[nt] real<lower = 0> ii;
array[nt] int<lower = 0> addl;
array[nt] int<lower = 0> ss;
array[nt] real rate;
array[nt] real weight;
vector<lower = 0>[nObsPK] cObs;
// data for population model
int<lower = 1> nSubjects;
array[nSubjects] int<lower = 1> start;
array[nSubjects] int<lower = 1> end;
// data for priors
real<lower = 0> CLHatPrior;
real<lower = 0> QHatPrior;
real<lower = 0> V1HatPrior;
real<lower = 0> V2HatPrior;
real<lower = 0> CLHatPriorCV;
real<lower = 0> QHatPriorCV;
real<lower = 0> V1HatPriorCV;
real<lower = 0> V2HatPriorCV;
}
transformed data{
vector[nObsPK] logCObs = log(cObs);
int nTheta = 4; // number of ODE parameters
int nIIV = 4; // parameters with IIV
int nCmt = 2; // number of compartments
array[nCmt] real biovar;
array[nCmt] real tlag;
for (i in 1:nCmt) {
biovar[i] = 1;
tlag[i] = 0;
}
}
parameters{
real<lower = 0> CLHat;
real<lower = 0> QHat;
real<lower = 0> V1Hat;
real<lower = 0> V2Hat;
real<lower = 0> sigma;
// IIV parameters
cholesky_factor_corr[nIIV] L;
vector<lower = 0>[nIIV] omega;
matrix[nIIV, nSubjects] etaStd;
}
transformed parameters{
row_vector[nt] cHat;
row_vector[nObsPK] cHatObs;
matrix[nCmt, nt] x;
array[nTheta] real<lower = 0> parms; // The [1] indicates the parameters are constant
// variables for Matt's trick
vector<lower = 0>[nIIV] thetaHat;
matrix<lower = 0>[nSubjects, nIIV] thetaM;
// Matt's trick to use unit scale
thetaHat[1] = CLHat;
thetaHat[2] = QHat;
thetaHat[3] = V1Hat;
thetaHat[4] = V2Hat;
thetaM = (rep_matrix(thetaHat, nSubjects) .*
exp(diag_pre_multiply(omega, L * etaStd)))';
for(i in 1:nSubjects) {
parms[1] = thetaM[i, 1] * (weight[i] / 70)^0.75; // CL, but is implementation of piece-wise correct like this for weight?
parms[2] = thetaM[i, 2] * (weight[i] / 70)^0.75; // Q
parms[3] = thetaM[i, 3] * (weight[i] / 70); // V1
parms[4] = thetaM[i, 4] * (weight[i] / 70); // V2
x[start[i]:end[i]] = pmx_solve_rk45(twoCptNeutModelODE, nCmt,
time[start[i]:end[i]],
amt[start[i]:end[i]],
rate[start[i]:end[i]],
ii[start[i]:end[i]],
evid[start[i]:end[i]],
cmt[start[i]:end[i]],
addl[start[i]:end[i]],
ss[start[i]:end[i]],
parms, biovar, tlag,
1e-6, 1e-6, 1e6);
cHat[start[i]:end[i]] = x[1, start[i]:end[i]] / parms[3]; // divide by V1
}
cHatObs = cHat[iObsPK];
}
model{
// Priors
CLHat ~ lognormal(log(CLHatPrior), CLHatPriorCV);
QHat ~ lognormal(log(QHatPrior), QHatPriorCV);
V1Hat ~ lognormal(log(V1HatPrior), V1HatPriorCV);
V2Hat ~ lognormal(log(V2HatPrior), V2HatPriorCV);
sigma ~ cauchy(0, 1);
// Parameters for Matt's trick
L ~ lkj_corr_cholesky(1);
to_vector(etaStd) ~ normal(0, 1);
// observed data likelihood
logCObs ~ normal(log(cHatObs), sigma);
}
generated quantities{
matrix[nCmt, nt] xPred;
array[nTheta] real<lower = 0> parmsPred;
row_vector[nt] cHatPred;
vector<lower = 0>[nObsPK] cHatObsCond;
row_vector<lower = 0>[nObsPK] cHatObsPred;
// Variables for IIV
matrix[nIIV, nSubjects] etaStdPred;
matrix<lower = 0>[nSubjects, nIIV] thetaPredM;
corr_matrix[nIIV] rho;
rho = L * L';
for(i in 1:nSubjects) {
for(j in 1:nIIV) {
etaStdPred[j, i] = normal_rng(0, 1);
}
}
thetaPredM = (rep_matrix(thetaHat, nSubjects) .*
exp(diag_pre_multiply(omega, L * etaStdPred)))';
for(i in 1:nSubjects) {
parmsPred[1] = thetaPredM[i, 1] * (weight[i] / 70)^0.75; // CL
parmsPred[2] = thetaPredM[i, 2] * (weight[i] / 70)^0.75; // Q
parmsPred[3] = thetaPredM[i, 3] * (weight[i] / 70); // V1
parmsPred[4] = thetaPredM[i, 4] * (weight[i] / 70); // V2
xPred[start[i]:end[i]] = pmx_solve_rk45(twoCptNeutModelODE, nCmt,
time[start[i]:end[i]],
amt[start[i]:end[i]],
rate[start[i]:end[i]],
ii[start[i]:end[i]],
evid[start[i]:end[i]],
cmt[start[i]:end[i]],
addl[start[i]:end[i]],
ss[start[i]:end[i]],
parmsPred, biovar, tlag,
1e-6, 1e-6, 1e6);
cHatPred[start[i]:end[i]] = xPred[1, start[i]:end[i]] / parmsPred[3]; // divide by V1
}
// predictions for observed data records
cHatObsPred = cHatPred[iObsPK];
for(i in 1:nObsPK) {
cHatObsCond[i] = exp(normal_rng(log(fmax(machine_precision(), cHatObs[i])), sigma));
cHatObsPred[i] = exp(normal_rng(log(fmax(machine_precision(), cHatObsPred[i])), sigma));
}
}
I have attached my nonmem-like data file. Would appreciate your support to get this model running.
thanks,
Andras