I would add up variances and not standard deviations…
1 Like
@wds15 You are obviously correct, it totally should have been real sigma = sqrt(square(sigma_absolute) + square(sigma_relative * mu[n])). This was based on my advice so thanks for noting and sorry for the confusion - assuming my older code is definitely OK and giving it as example is probably not a good policy going forward :-). (and I definitely hope @dbarneche checks what I propose, I definitely make a lot of mistakes).
Other than that looks roughly good, some more minor things I noticed:
Those are very wide priors, especially for sigma_relative, I think you can safely assume that a priori the relative error is not more than 1 (i.e. the error is not 100%), giving you something like sigma_relative ~ normal(0, 0.5); or even tighter. Similarly, knowing the error of your measuring device could give you better prior for the sigma_absolute.
This is also a bit weird - the constraint will only make the sampler throw errors at you when it is invalidated, and break the sampler if non-negligible posterior mass is outside the constraint. Are you sure Eff cannot physically be larger than 1? (I’ve discussed that in a previous post here, but not sure if you changed the model or you’ve seen a mistake in my reasoning).
1 Like
Thanks again @martinmodrak and @wds15! Even after adding the corrections above (i.e. adding variances rather than standard deviations, tighter priors on sigma), the peak is still slightly underestimated when compared to the model without an absolute error term only.
The crux of the problem is the efficiency equation. I came up with this efficiency reasoning myself borrowing some concepts from ecological theory (my field of research), so while I’d say it makes lots of sense to me, I admit I could be missing something.
If we start with concentration \mu = 0 at time T_0, and the dynamics are dictated by this single dose, first order absorption model, the only parameter in the equation that conceptually drives the drug absorption (as far as I understand!) is K_a. K_e and C_l are parameters related to elimination processes.
So, if the conceptual interpretation above is correct, K_a T sets the maximum possible concentration at any time point (even in the hypothetical absence of elimination processes). So, whatever the concentration value at the peak P is (or any concentration at any time point), it has to be lower than K_a T, even if mathematically higher values of Q (originally D/C_l) allow for E > 1.
Please let me know if you disagree. I’m not a pharmacologist, so I am going to reach out to more experienced folks to confirm this reasoning is correct!
I thought so for a while, but that is NOT correct:
library(tidyverse)
D = 10
Ka = 1
Cl = 1
Ke = 0.1
tibble(t = seq(0,10,length.out = 100)) %>%
mutate(mu = D * Ka * (exp(-Ke * t) - exp(-Ka * t)) / (Cl * (Ka - Ke))) %>%
ggplot(aes(t, mu)) + geom_line() + geom_abline(slope = Ka, color = "blue", linetype = "dashed")
The dashed line is K_a T, the curve is \mu(T)
1 Like
Hi @martinmodrak,
Thanks for this. As mentioned before, I appreciate that mathematically you can get cases where E > 1. In your plot above, for instance, you only see that because D was very high relative to C_l. But if you were to set D = 1 for example, you’d get E < 1 over the entire range of t.
I wonder two things:
1 - My proposed equation is correct, and all I am missing are the defined natural constraints of this equation that would set the values of E to always be [0,1]. I haven’t been able to find some in-depth discussion on the web about what all these parameters mean, and how (or if) they constrain one another.
2 - If there’s another efficiency equation for E that explicitly limits its value to be [0,1] regardless of the parameter space. In my mind, even if the E formula I proposed is wrong, the question is still valid, i.e. what was the efficiency of the organism in absorbing the drug up to the peak? Efficiencies by definition are never > 1. I’m not a mathematician though, so I wouldn’t know if there’s a method that allows me to derive such an equation that is independent of the parameter space. If you have ideas though, I’d love to hear them!
Thanks again!
Turns out I was using wrong formula (withoute K_e in the numerator). With the correct one, you can make E arbitrarily large even for fixed D/C_l, say:
D = 1
Ka = 1
Cl = 0.5
Ke = 10
T_p = (log(Ka) - log(Ke)) / (Ka - Ke)
P_rel = Ke * Ka * (exp(-Ke*T_p) - exp(-Ka*T_p)) / (Ka - Ke)
P = P_rel * D / Cl
tibble(t = seq(0,1,length.out = 100)) %>%
mutate(mu = D *Ke * Ka * (exp(-Ke * t) - exp(-Ka * t)) / (Cl * (Ka - Ke))) %>%
ggplot(aes(t, mu)) + geom_line() + geom_abline(slope = Ka, color = "blue", linetype = "dashed")
And note that when modelling you generally shouldn’t put hard constraints based on domain knowledge (e.g that D < C_l), this should rather be represented as a soft constraint in prior. Hard constraints generally belong to unphysical/mathematically impossible situations, which E > 1 isn’t.
Hi @martinmodrak,
I wanted to follow up on this problem, but I’m not sure if I should open a new thread. I’m assuming this next question could be more of a generalised question in stan.
I have successfully implemented the model above in a hierarchical fashion, with patient as a random effect on all parameters. However, let’s assume though that my patients were split within two treatment levels, placebo vs. the real drug. What I would like to do now is expand the code below to test whether efficiency, E (Eff in the code below), is substantially different between treatment levels (placebo vs. real drug), like a simple anova in R. Eff ~ treatment - 1. How can I do this given that E is declared inside the transformed parameters block rather than being a variable informed in the data block? Unfortunately I can’t (I think) include this as a nested random effect either (i.e. patient nested in treatment) because I only have two levels of treatment. Any help would be extremely appreciated!
data {
int<lower=1> N; // number of observations
vector[N] Y; // response variable
// covariate vectors
vector[N] Time_days;
// data for group-level effects of ID 1
int<lower=1> N_1; // number of patients
int<lower=1> M_1; // number of variables with random effects (3)
int<lower=1> J_1[N]; // vector of patient ID
}
parameters {
real<lower=0> Ka; // population-level effects
real logitKe; // population-level effects
real<lower=0> Q; // population-level effects
real<lower=0> sigma; // residual SD
vector<lower=0>[M_1] sd_1; // group-level standard deviations
matrix[M_1, N_1] z_1; // unscaled group-level effects
// cholesky factor of correlation matrix
cholesky_factor_corr[M_1] L_1;
}
transformed parameters {
// save average Ke; elimination rate, needs to be positive and lower than Ka
real<lower=0> Ke = Ka / (1 + exp(-logitKe));
// group-level effects
matrix[N_1, M_1] r_1 = (diag_pre_multiply(sd_1, L_1) * z_1)';
vector<lower=0>[N_1] r_1_Ka = r_1[, 1] + Ka;
vector<lower=0>[N_1] r_1_Ke = r_1_Ka ./ (1 + exp(-(r_1[, 2] + logitKe))); // elimination rate, needs to be positive and lower than Ka
vector<lower=0>[N_1] r_1_Q = r_1[, 3] + Q;
vector<lower=0,upper=1>[N_1] Eff;
for(i in 1:N_1) {
Eff[i] = r_1_Q[i] * (r_1_Ka[i] - r_1_Ke[i]) * (r_1_Ka[i] / r_1_Ke[i])^(r_1_Ka[i] / (r_1_Ke[i] - r_1_Ka[i])) / log(r_1_Ka[i] / r_1_Ke[i]); // constrain E mathematically based on Q, Ke, and Ka
}
}
model {
vector[N] mu;
mu = r_1_Q[J_1] .* r_1_Ka[J_1] .* r_1_Ke[J_1] .* (exp(-r_1_Ke[J_1] .* Time_days) - exp(-r_1_Ka[J_1] .* Time_days)) ./ (r_1_Ka[J_1] - r_1_Ke[J_1]);
// priors including all constants
Ka ~ normal(0.4, 2);
logitKe ~ normal(0, 5);
Q ~ normal(5, 2);
sigma ~ cauchy(0, 5);
sd_1 ~ cauchy(0, 5);
to_vector(z_1) ~ normal(0, 1);
L_1 ~ lkj_corr_cholesky(3);
// likelihood including all constants
Y ~ normal(mu, sigma);
}
So far what I have managed to do is to use the posterior samples of Eff after running the model, and for each iteration I calculated the mean value of Eff per treatment level (i.e. across patients within each treatment level), in order to have a posterior distribution of Eff per treatment level. But I’m wondering whether that’s throwing away useful uncertainty info?
Thanks for your help!
You can’t really treat derived quantities (such as Eff) as dependent on group. What you can do is make the actual parameters (Ka, Ke) depend on group. Than you can use the posterior samples of the parameters per group to estimate efficiency per group. Does that make sense?
1 Like
Yep, that’s what I imaged. Thanks!
Just one final update on this—sorry it took me almost a year to post a proper closure on the topic.
@martinmodrak your insights here were much appreciated. It turns out I had some of the parameter units wrong, and the efficiency formula now does yield values bounded between 0 and 1, so there’s no need to constrain a transformed parameter, as you emphasised.
The units are:
K_a and K_e in 1 / minute
D in mg
We can calculate the mass of administered drug, m(t), that was absorbed into the body since time t = 0 as
m(t) = D (1 - e^{-K_at})
So, when t is large, m_a(t) \equiv D when a drug is 100% absorbable. So the efficiency at the peak time, t_p, becomes easily calculated, i.e. the mass of drug that is still in circulation, m_c(t_p), relative to m(t_p). Because
m_c(t_p) = \frac{DK_a\left(e^{-K_et_p} - e^{-K_at_p}\right)}{K_a - K_e}, then
\epsilon(t_p) = m_c(t_p)/m(t_p) becomes independent of D, and is bounded between 0 and 1 no matter what the combinations of K_a and K_e are.
So now \epsilon can either be calculated inside the generated quantities block or outside in R.
4 Likes