Hi, this is my first post on here and am newbie to stan, so please let me know if its not detailed enough, or way too complicated than it needs to be, as I’m using a modified version of a template i found on metrum tutorial ‘introduction to bayesian modelling using stan’.
I’m attempting to fit a multiple dose 1 cmt- PK model and here is my stan model file;
functions{
vector oneCptModel1(real dt, vector init, real amt, int cmt, int evid,
real CL, real V){
vector[1] x;
real alpha;
alpha = CL / V;
x = rep_vector(0.0, 1);
x[1] = init[1] * exp(-alpha * dt);
if (evid == 1) x[cmt] = x[cmt] + amt;
return x;
}
matrix oneCptModel(real[] time, real[] amt, int[] cmt, int[] evid,
real CL, real V){
vector[1] init;
real dt;
real t0;
matrix[size(time), 1] result;
int nt;
nt = size(time);
init = rep_vector(0, 1);
t0 = time[1];
for(i in 1:nt){
dt = time[i] - t0;
init = oneCptModel1(dt, init, amt[i], cmt[i], evid[i],
CL, V);
result[i, 1] = init[1];
t0 = time[i];
}
return result;
}
}
data{
int nSubjects;
int nt;
int nObs;
int iObs[nObs];
real amt[nt];
int cmt[nt];
int evid[nt];
int start[nSubjects];
int end[nSubjects];
real time[nt];
vector[nObs] cObs;
int weight[nSubjects];
}
transformed data{
vector[nObs] logCObs;
int nti[nSubjects];
int nRandom;
logCObs = log(cObs);
nRandom = 2;
}
parameters{
real<lower = 0.001> CLHat;
real<lower = 0.001> VHat;
vector<lower = 0>[nRandom] omega;
real<lower = 0> sigma;
vector[nRandom] coeff[nSubjects];
vector<lower = 0> [nSubjects] eta1;
vector<lower = 0> [nSubjects] eta2;
}
transformed parameters{
vector<lower = 0.001>[nRandom] thetaHat;
vector[nRandom] ltv;
cov_matrix[nRandom] Omega;
real<lower = 0> CL[nSubjects];
real<lower = 0> V[nSubjects];
vector [nt] cHat;
vector [nObs] cHatObs;
matrix [nt, 1] x;
vector[nRandom] logtheta[nSubjects];
thetaHat[1] = CLHat;
thetaHat[2] = VHat;
ltv[1] = log(thetaHat[1]);
ltv[2] = log(thetaHat[2]);
Omega = diag_matrix(omega);
for(j in 1:nSubjects){
logtheta[j,1] = (ltv[1] + eta1[j]);
logtheta[j,2] = (ltv[2] + eta2[j]);
CL[j] = exp(logtheta[j, 1])*(weight[j]/1000);
V[j] = exp(logtheta[j, 2])*(weight[j]/1000);
x[start[j]:end[j],] = oneCptModel(time[start[j]:end[j]],
amt[start[j]:end[j]],
cmt[start[j]:end[j]],
evid[start[j]:end[j]],
CL[j], V[j]);
cHat[start[j]:end[j]] = x[start[j]:end[j], 1] ./ V[j];
}
cHatObs = cHat[iObs];
}
model{
CLHat ~ uniform(1, 20);
VHat ~ uniform(1, 10);
for (j in 1:nSubjects){
eta1[j] ~ normal(0, omega[1]);
eta2[j] ~ normal(0, omega[2]);
}
omega[1] ~ uniform(0.1, 1);
omega[2] ~ uniform(0.1, 1);
sigma ~ uniform(0.1, 10);
## Inter-individual variability
// logtheta ~ multi_normal(log(thetaHat), Omega);
logtheta[1] ~ normal(ltv[1],omega[1]);
logtheta[2] ~ normal(ltv[2],omega[2]);
// for (i in 1:nObs)
// if (evid[i] == 0){
// logCObs[i] ~ normal(log(cHatObs[i]), sigma);
// }
for (i in 1: nObs){
logCObs[i] ~ normal(log(cHatObs[i]), sigma);
}
}
generated quantities{
vector[nRandom] logthetaPred[nSubjects];
vector<lower = 0>[nt] cHatPred;
real <lower=0>cObsCond[nt];
real <lower=0>cObsPred[nt];
real<lower = 0> CLPred[nSubjects];
real<lower = 0> VPred[nSubjects];
real log_lik[nObs];
matrix[nt, 2] xPred;
for(j in 1:nSubjects){
logthetaPred[j] = multi_normal_rng(log(thetaHat), Omega);
CLPred[j] = exp(logthetaPred[j, 1]);
VPred[j] = exp(logthetaPred[j, 2]);
xPred[start[j]:end[j],] = oneCptModel(time[start[j]:end[j]],
amt[start[j]:end[j]],
cmt[start[j]:end[j]],
evid[start[j]:end[j]],
CLPred[j], VPred[j]);
cHatPred[start[j]:end[j]] = xPred[start[j]:end[j], 2] ./ V[j];
}
for(i in 1:nt){
if(time[i] == 0){
cObsCond[i] = 0;
cObsPred[i] = 0;
}else{
cObsCond[i] = exp(normal_rng(log(cHat[i]), sigma));
cObsPred[i] = exp(normal_rng(log(cHatPred[i]), sigma));
}
}
}
This is my script file, with initial values:
NeoGent <- read.csv(file=“nGent_data7a-163.csv”, header = TRUE, sep = “,”)
library(data.table)
library(dplyr)
NeoGent <- as.data.table(NeoGent)
NeoGent$CMT <- 1
NeoGentEVID10 <- as.data.frame(NeoGent[NeoGent$EVID != “2”]) # training data set, incl 1+0
NeoGentEVID10 <- NeoGentEVID10%>%
mutate(DV = ifelse(EVID, NA, DV))%>%
arrange(ID, TIME, desc(EVID))
nt <- nrow(NeoGentEVID10)
start <- (1:nt)[!duplicated(NeoGentEVID10$ID)]
end <- c(start[-1] - 1, nt)
iObs <- with(NeoGentEVID10, (1:nrow(NeoGentEVID10))[!is.na(DV) & EVID == 0])
nObs <- length(iObs)
nSubjects = length(unique(NeoGentEVID10$ID))
attach(NeoGentEVID10)
data <- with(NeoGentEVID10,
list(
nSubjects = nSubjects,
nt = nt,
nObs = nObs,
iObs = iObs,
amt = AMT,
cmt = CMT,
evid = EVID,
start = start,
end = end,
time = TIME,
cObs = DV[iObs],
weight = WT[!duplicated(ID)]
))
init <- function(){
list(CLHat = 0.046,
VHat = 0.73,
omega = c(0.0576,0.0729),
sigma = runif(1, 0.09, 2),
eta1 = rep(0.5,nSubjects),
eta2 = rep(0.5,nSubjects)
#coeff = matrix(rep(c(0.002,0.0004),ea = nObs),nrow = nObs)
)
}
parametersToPlot <- c(“CLHat”, “VHat”,
“sigma”, “omega”)
otherRVs <- c(“cObsCond”, “cObsPred”, “CL”, “V”) ##, “log_lik”)
parameters <- c(parametersToPlot, otherRVs)
nChains <- 4
nPost <- 1000 ## Number of post-burn-in samples per chain after thinning
nBurn <- 1000 ## Number of burn-in samples per chain after thinning
nThin <- 1
nIter <- (nPost + nBurn) * nThin
nBurnin <- nBurn * nThin
fit <- stan(‘grasela.stan’,
data = data,
pars = parameters,
iter = nIter,
warmup = nBurnin,
thin = nThin,
init = init,
chains = nChains)
This is the error message i get:
Rejecting initial value:
Error evaluating the log probability at the initial value.
Exception: normal_lpdf: Location parameter is nan, but must be finite! (in ‘modelea8d5760c71b_grasela’ at line 137)
Initialization between (-2, 2) failed after 100 attempts.
Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] “Error in sampler$call_sampler(args_list[[i]]) : Initialization failed.”
FYI, line 136-137 corresponds to the line:
for (i in 1: nObs){
logCObs[i] ~ normal(log(cHatObs[i]), sigma);
I have tried to play around with the bounds and the parameters of the distributions, but no luck. Any help will be greatly appreciated!!