I don’t know is this is the right place but while trying to implement multi dose IV I happened to look at torsten codes and was confused by IV rate implementation. It applies both to 1cmp and 2cmp models.
In the Pred1_one function there is a statement which updates concentration in central compartment:
if ((init[1] != 0) || (rate[1] != 0)) {
a[0] = 1;
pred(0, 1) += PolyExp(x=dt, dose=init[1], rate=0, xinf=0, tau=0, ss=false, a, alpha, n=1) +
PolyExp(dt, 0, rate[1], dt, 0, false, a, alpha, 1);
}
When amt (and init) has units of mass and rate has units of mass/time, the pred functions of Torsten all calculate mass, not concentration. The user is responsible for converting mass to concentration in the Stan model code, Thus the units are consistent.
Regarding a mail list for Torsten, there is none yet. My guess is there’ll be one once we expect higher traffic, and v1.0 gets released. I’d say it’s appropriate to post issues on github, github.com/metrumresearchgroup/stan. If you post on the Stan forum, we’ll respond + other Stan users and developers can step in, which is if course a plus.
Thanks. The next item I am struggling with is the parameters that change
over the time. The signature of one cpt model requires that theta (the
parameter next to biovar) has to be real[] or real[,]. While the attached
code runs for real[] case it doesn’t for real[,] case (commented code). Any
idea where is the problem?
I think your intention is to only populate theta with parameters for individual j, but the loop index corresponds to the entire data set since N_rest appears to be the total number of event times for all individuals. Maybe you want something like “for(i in 1:nt[j])” for the loop index where nt[j] is the number of event times for the jth individual, i.e.,
for(i in 1:nSubjects)
ntj[i] = end_rest[j] - start_rest[j] + 1;
Since you’re probably dealing with a ragged array, you will also need to specify the same index range when you pass theta as an argument, i.e.,
Finally declare theta with a size corresponding to the largest number of event times for a single individual, i.e.,
real<lower = 0> theta[max(ntj), 3];
Thanks. I have corrected this. There is still relatively small mismatch
between torsten and my implementation of one compartmental model. Both are
attached in case you are interested. Torsten model gives much more
divergencies even if I use Matt’s trick and sd of estimates is much larger.
Thanks for these. Are you able to share the data? Otherwise it’s hard for me to diagnose the discrepancy. The fact the sd estimates are much larger makes me suspect the PK model calculations are different event given the same parameter values. Otherwise I would expect them to converge to the same posterior.
One possibility for mismatch is the way IV is handled. In my
implementation whenever IV occurs I add two dosing events: one with real
dose and real rate and another 1h later with 0 dose and 0 rate. Since it
is possible that some labs will occur within this 1h window I sort based
on time. I don’t know how this is handled in torsten.
If you can’t share the data set, would you share an excerpt of the data showing how the dosing is specified? I would really like to understand if there is a bug I need to fix or just a misunderstanding about how to specify dosing.
Torsten uses NONMEM conventions for dose specifications. If I understand correctly, you are trying to specify a constant rate input with a 1 hour duration. Is that correct? If so, that is specified with a single dosing event record in which rate = the input rate, amt = rate * (input duration – 1 h in this case), evid = 1, cmt = compartment receiving input (2 for PKModelOneCpt if it is IV input).