OK, I think I’m almost there! I have been able to implement the code in stan with a gaussian distribution instead of gamma, with some key modifications:
1 – I have substituted D / C_l for a constant parameter Q, just for simplicity.
2 – In the transformed parameters block, I have included the efficiency E as a bounded parameter [0,1] based on Q, K_a, and K_e (assuming they’re all real):
E = \frac{P}{K_a T_p} = \frac{Q (Ka - K_e) (K_a/K_e)^{K_a/(K_e - K_a)}}{log(K_a) - log(K_e)}
3 – I fitted K_e in its logit form, K_e = \frac{K_a}{1 + e^{-\textrm{logit}(K_e)}}, to ensure that K_a > K_e;
4 – The formula is now fitted using T_p informed from data, substituting K_a - K_e in the original formulation with K_a - K_e = \frac{log(K_a) - log(K_e)}{T_p}:
\mu(T,T_p) = \frac{\frac{Q K_a^2}{1 + e^{-\textrm{logit}(K_e)}} (e^{-\frac{K_a T}{1 + e^{-\textrm{logit}(K_e)}} } - e^{-K_a T})}{\frac{log(K_a) - log(K_e)}{T_p}}
The only thing missing now is the error handling; the errors towards the tails on the data are most definitely underestimated, as pointed out by @martinmodrak – please see figure at the end. I am no longer fitting the function on the log scale, so I’m assuming it’s simpler to implement the absolute + relative error? I would really appreciate your help here implementing the relative error in the stan code below. It confuses me a bit because sigma here is real, and I do not understand how I can make operations in stan between a vector (which would be the case of Y relative to \mu) and a real. Thanks so much for all your guidance!
stan code:
data {
int<lower=1> N; // number of observations
vector[N] Y; // concentration
int<lower=1> K_Ka; // number of population-level effects
matrix[N, K_Ka] X_Ka; // population-level design matrix
int<lower=1> K_logitKe; // number of population-level effects
matrix[N, K_logitKe] X_logitKe; // population-level design matrix
int<lower=1> K_Q; // number of population-level effects
matrix[N, K_Q] X_Q; // population-level design matrix
vector[N] Time_days; // Time
vector[N] Time_peak_days; // Time to peak
}
parameters {
vector<lower=0>[K_Ka] Ka; // uptake rate
vector[K_logitKe] logitKe; // elimination rate (logit); Ke = (Ka / (1 + exp(-logitKe))); ensures that Ka > Ke
vector<lower=0>[K_Q] Q; // constant (originally D / Cl)
real<lower=0> sigma; // residual SD
}
transformed parameters {
vector<lower=0,upper=1>[1] Eff; // efficiency
Eff[1] = Q[1] * (Ka[1] - (Ka[1] / (1 + exp(-logitKe[1])))) * (Ka[1] / (Ka[1] / (1 + exp(-logitKe[1]))))^(Ka[1] / ((Ka[1] / (1 + exp(-logitKe[1]))) - Ka[1])) / log(Ka[1] / (Ka[1] / (1 + exp(-logitKe[1])))); // constrain E mathematically based on Q, logitKe, and Ka
}
model {
vector[N] nlp_Ka = X_Ka * Ka;
vector[N] nlp_logitKe = X_logitKe * logitKe;
vector[N] nlp_Q = X_Q * Q;
vector[N] mu;
for (n in 1:N) {
mu[n] = nlp_Q[n] * nlp_Ka[n] * (nlp_Ka[n] / (1 + exp(-nlp_logitKe[n]))) * (exp(-(nlp_Ka[n] / (1 + exp(-nlp_logitKe[n]))) * Time_days[n]) - exp(-nlp_Ka[n] * Time_days[n])) / (log(nlp_Ka[n] / (nlp_Ka[n] / (1 + exp(-nlp_logitKe[n])))) / Time_peak_days[n]);
}
// priors
target += normal_lpdf(Ka | 0, 10)
- 1 * normal_lccdf(0 | 0, 10);
target += normal_lpdf(logitKe | 0, 10)
- 1 * normal_lccdf(0 | 0, 10);
target += normal_lpdf(Q | 0, 10)
- 1 * normal_lccdf(0 | 0, 10);
target += normal_lpdf(Eff | 0.5, 0.5)
- 1 * log_diff_exp(normal_lcdf(1 | 0.5, 0.5), normal_lcdf(0 | 0.5, 0.5));
target += student_t_lpdf(sigma | 3, 0, 10)
- 1 * student_t_lccdf(0 | 3, 0, 10);
// likelihood
target += normal_lpdf(Y | mu, sigma);
}
R code:
library(rstan)
rstan::rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
data <- read.csv(text =
'Drug_concentration_mg.ml,Dose_mg,Time_days
0.0006773291,13.57695,0.01
0.4894254798,13.57695,0.88
1.0282572998,13.57695,1.24
1.6008287008,13.57695,1.77
1.5079260525,13.57695,2.10
1.4089264576,13.57695,2.85
1.1498402861,13.57695,3.93
0.8579193906,13.57695,4.95
0.5450494898,13.57695,7.05
0.2563178348,13.57695,11.97
0.1634733689,13.57695,15.98
0.0706932426,13.57695,23.74
0.0774802343,13.57695,28.96
0.0364683985,13.57695,53.01',
header = TRUE, stringsAsFactors = FALSE)
createParMat <- function (data) {
mat <- matrix(rep(1, nrow(data)))
colnames(mat) <- 'Intercept'
attr(mat, 'assign') <- 0
mat
}
brmData <- list('N' = nrow(data),
'Y' = data$Drug_concentration_mg.ml,
'Time_days' = data$Time_days,
'Time_peak_days' = rep(data$Time_days[which.max(data$Drug_concentration_mg.ml)], nrow(data)),
'K_Ka' = 1,
'X_Ka' = createParMat(data),
'K_logitKe' = 1,
'X_logitKe' = createParMat(data),
'K_Q' = 1,
'X_Q' = createParMat(data))
set.seed(10)
model <- rstan::stan('tmp_stancode_gaussian_logitKe_Tp_E.stan', data = brmData, chains = 3, cores = 3, iter = 1.5e4, warmup = 7.5e3, control = list(adapt_delta = 0.99, max_treedepth = 15))
# add time to peak, and peak, then plot the data
pars <- as.data.frame(t(summary(model, pars = c('Ka[1]', 'logitKe[1]', 'Q[1]'))$summary[, 'mean', drop = FALSE]))
names(pars) <- gsub('[1]', '', names(pars), fixed = TRUE)
pars$Ke <- pars$Ka / (1 + exp(-pars$logitKe))
T_p <- log(pars$Ka / pars$Ke) / (pars$Ka - pars$Ke) # time to peak in days
P <- pars$Q * pars$Ka * pars$Ke * (exp(-pars$Ke * T_p) - exp(-pars$Ka * T_p)) / (pars$Ka - pars$Ke) # peak in mg/ml
E <- P / (pars$Ka * T_p)
xseq <- seq(min(data$Time_days), max(data$Time_days), length.out = 200)
posts <- extract(model)
predictions <- matrix(0, nrow(posts[['Ka']]), length(xseq))
for (i in 1:nrow(predictions)) {
Ke <- posts[['Ka']][i] / (1 + exp(-posts[['logitKe']][i]))
predictions[i, ] <- posts[['Q']][i] * Ke * posts[['Ka']][i] * (exp(-Ke * xseq) - exp(-posts[['Ka']][i] * xseq)) / (posts[['Ka']][i] - Ke)
}
quantiles <- apply(predictions, 2, quantile, probs = c(0.025, 0.975))
modelPrediction <- data.frame('mean' = apply(predictions, 2, mean), 'Q2.5' = quantiles['2.5%', ], 'Q97.5' = quantiles['97.5%', ])
transparentCol <- rgb(255, 99, 71, max = 255, alpha = 125)
plot(data$Time_days, data$Drug_concentration_mg.ml, ylim = range(c(data$Drug_concentration_mg.ml, unlist(modelPrediction))), xlab = 'Time (days)', ylab = 'Drug concentration (mg / ml)', las = 1, pch = 16, col = 'tomato')
polygon(c(xseq, rev(xseq), xseq[1]), c(modelPrediction$Q2.5, rev(modelPrediction$Q97.5), modelPrediction$Q2.5[1]), border = NA, col = transparentCol) # notice negative values on posterior predictions
lines(xseq, modelPrediction$mean, lty = 2, lwd = 1.5, col = 'tomato')
