stan and PK modeling

[Linas Mockus]

Hi,

Is it possible to implement such error model in stan:
y_i=model_i+model_i*err1+err2 where err1~N(0,sigma1) and err2~N(0,sigma2)

NONMEM allows such error model:
$ERROR
IPRED=F
Y=F+F*ERR(1)+ERR(2)

Thanks a bunch for any suggestions,
Linas

[Bob Carpenter]

You can just do this directly in Stan. First
note that the usual linear regression

y[i] = alpha * x[i] + beta + epsilon[i];
epsilon[i] ~ normal(0,sigma);

translates into Stan as follows:

y[i] ~ normal(alpha * x[i] + beta, sigma);

The epsilon[i] term is marginalized out, but you
can recover it as a generated quantity as
(y[i] - (alpha * x[i] + beta)) if you need the residuals
for something.

You can do the same thing with your system, which
I'm writing fully indexed as:

y[i] = model[i] + (err1[i] * model[i]) + err2[i];

err1[i] ~ normal(0,sigma1);
err2[i] ~ normal(0,sigma2);

This becomes

y[i] ~ normal(model[i] + (err1[i] * model[i]), sigma2);
err1[i] ~ normal(0, sigma1);

At least I think that works out to what you want. With
essentially two error terms and one that varies by item, I'm
curious as to how it'll fit.

What priors were you thinking of using for the sigma?

- Bob

> --
> You received this message because you are subscribed to the Google Groups "stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Linas Mockus]

Thanks, Bob.  I was planning to use Cauchy (0,5) for sigmas.  Based on examples in stan I understood that it is current best practice.

Does it make sense to write the model as (no index on err1):

  y[i] ~ normal(model[i] + (err1 * model[i]), sigma2);
  err1 ~ normal(0, sigma1);

Linas

[Bob Carpenter]

On Mar 26, 2014, at 1:54 AM, Linas Mockus <linasm...@gmail.com> wrote:

> Thanks, Bob. I was planning to use Cauchy (0,5) for sigmas. Based on examples in stan I understood that it is current best practice.
>

It depends on what you know about the scale. If you
expect unit normal noise, then it's overly broad, but if
you expect noise on the scale of 100, it's overly narrow.

> Does it make sense to write the model as (no index on err1):
> y[i] ~ normal(model[i] + (err1 * model[i]), sigma2);
> err1 ~ normal(0, sigma1);
>

Yes, you could do that. Then it's not really noise so much as you
want (1 + err1) as a multiplier of model[i]. Why would you want this
rather than just rescaling model directly?

- Bob

[Linas Mockus]

On 3/26/2014 12:57 PM, Bob Carpenter wrote:
>> Thanks, Bob. I was planning to use Cauchy (0,5) for sigmas. Based on examples in stan I understood that it is current best practice.
>>
> It depends on what you know about the scale. If you
> expect unit normal noise, then it's overly broad, but if
> you expect noise on the scale of 100, it's overly narrow.

I think the noise should be unit normal. Perhaps it makes sense to try
sigma ~ Uniform (0, 1).

>> Does it make sense to write the model as (no index on err1):
>> y[i] ~ normal(model[i] + (err1 * model[i]), sigma2);
>> err1 ~ normal(0, sigma1);
>>
> Yes, you could do that. Then it's not really noise so much as you
> want (1 + err1) as a multiplier of model[i]. Why would you want this
> rather than just rescaling model directly?

With regards to err1 with index I was concerned that it may not converge.

Yet another question. PK models often have lags. It means that nothing
happens until some time t0. For example, model(t)=0 when t < t0;
model(t)=exp(t-t0)-1 when t>=t0. As such there is C0 continuity but
gradient is discontinuous at t=t0. Obviously I want t0 to be a
parameter to estimate.

Linas

[Bob Carpenter]

On Mar 26, 2014, at 8:14 PM, Linas Mockus <linasm...@gmail.com> wrote:

> On 3/26/2014 12:57 PM, Bob Carpenter wrote:
>>> Thanks, Bob. I was planning to use Cauchy (0,5) for sigmas. Based on examples in stan I understood that it is current best practice.
>>>
>> It depends on what you know about the scale. If you
>> expect unit normal noise, then it's overly broad, but if
>> you expect noise on the scale of 100, it's overly narrow.
> I think the noise should be unit normal. Perhaps it makes sense to try sigma ~ Uniform (0, 1).
>>> Does it make sense to write the model as (no index on err1):
>>> y[i] ~ normal(model[i] + (err1 * model[i]), sigma2);
>>> err1 ~ normal(0, sigma1);
>>>
>> Yes, you could do that. Then it's not really noise so much as you
>> want (1 + err1) as a multiplier of model[i]. Why would you want this
>> rather than just rescaling model directly?
> With regards to err1 with index I was concerned that it may not converge.
>

That's a valid concern because you'd get a redundant parameterization.
I just don't think it's doing anything without an index other than
rescaling a bit in a way that's also going to cause a near multiplicative
invariance between all of the model[i] and err1.

> Yet another question. PK models often have lags. It means that nothing happens until some time t0. For example, model(t)=0 when t < t0; model(t)=exp(t-t0)-1 when t>=t0. As such there is C0 continuity but gradient is discontinuous at t=t0. Obviously I want t0 to be a parameter to estimate.

I don't think model(t) as written is discontinuous at t0.

But in general, you can give particular data points different models.
For instance, in a time series, the first point can get a diffuse
prior, the next point one based on just a single previous point and
so on up to the general loop.

- Bob

[Linas Mockus]

>> Yet another question. PK models often have lags. It means that nothing happens until some time t0. For example, model(t)=0 when t < t0; model(t)=exp(t-t0)-1 when t>=t0. As such there is C0 continuity but gradient is discontinuous at t=t0. Obviously I want t0 to be a parameter to estimate.
> I don't think model(t) as written is discontinuous at t0.
>
> But in general, you can give particular data points different models.
> For instance, in a time series, the first point can get a diffuse
> prior, the next point one based on just a single previous point and
> so on up to the general loop.

Thanks, Bob. Is there any reference I could read more on giving

particular data points different models.

Linas

[Bob Carpenter]

I'm not sure what you mean by a reference --- my understanding
is that it's common practice in time-series models.

If you have a time series of data points y[1],...,y[N],
and you want to in general model y[n] as a function of
y[n-1], y[n-2], ..., y[n-K], then you have to have a different
model for y[1],...y[K], because you don't have the relevant
predictors.

- Bob

[Linas Mockus]

I found the solution for 2-compartmental model in WinBUGS:
        for (t in 1 : N) {
            model_c <- B1 * exp (-lambda1 * time [t]) +
                B2 * exp (-lambda2 * time [t]) +
                B3 * exp (-ka * time [t]);
            sigma <- sqrt (sigma_m * sigma_m * model_c * model_c +
                sigma_a * sigma_a);
            concentration [t] ~ normal (model_c, sigma);
        }

However for some combination of parameters model_c is negative as well as concentration.  When using explicit likelihood function I assign very low number to the likelihood whenever model_c<0 - don't know if it is a good practice but basically I tell the solver not to explore infeasible set of  parameters.  How to do this in stan?

Also, I generate predictive posterior in generated quantities block:
        for (t in 1 : N_pred) {
            model_c <- B1 * exp (-lambda1 * time_pred [t]) +
                B2 * exp (-lambda2 * time_pred [t]) +
                B3 * exp (-ka * time_pred [t]);
            sigma <- sqrt (sigma_m * sigma_m * model_c * model_c +
                sigma_a * sigma_a);
            concentration_pred [t] <- normal_rng (model_c, sigma);
        }
For some reason traceplot gives for stanfit object:
Error in mcmc.list(x) : Arguments must be mcmc objects

I appreciate any insight on those two items.

Linas

[Bob Carpenter]

On Mar 29, 2014, at 9:45 PM, Linas Mockus <linasm...@gmail.com> wrote:

> I found the solution for 2-compartmental model in WinBUGS:
> for (t in 1 : N) {
> model_c <- B1 * exp (-lambda1 * time [t]) +
> B2 * exp (-lambda2 * time [t]) +
> B3 * exp (-ka * time [t]);
> sigma <- sqrt (sigma_m * sigma_m * model_c * model_c +
> sigma_a * sigma_a);
> concentration [t] ~ normal (model_c, sigma);
> }
>
> However for some combination of parameters model_c is negative as well as concentration. When using explicit likelihood function I assign very low number to the likelihood whenever model_c<0 - don't know if it is a good practice but basically I tell the solver not to explore infeasible set of parameters. How to do this in stan?

Priors on B1, B2, B3 would be the way to go, or even hard constraints
making them positive if that makes sense for the model. Or equivalently,
putting them inside the exp() function as additive terms.

There's a square(x) function, by the way, which is a bit more
efficient than x * x. (And we really need to add normal distributions
parameterized with variance and precision, which would avoid that
sqrt operation, too --- it's just going to get squared again inside
the density!).

You could also constrain concentration to be positive, though you'll
still get odd behavior with model_c going below 0 (it's symmetric around
0 in the sigma term).

>
> Also, I generate predictive posterior in generated quantities block:
> for (t in 1 : N_pred) {
> model_c <- B1 * exp (-lambda1 * time_pred [t]) +
> B2 * exp (-lambda2 * time_pred [t]) +
> B3 * exp (-ka * time_pred [t]);
> sigma <- sqrt (sigma_m * sigma_m * model_c * model_c +
> sigma_a * sigma_a);
> concentration_pred [t] <- normal_rng (model_c, sigma);
> }
> For some reason traceplot gives for stanfit object:
> Error in mcmc.list(x) : Arguments must be mcmc objects

Could you write back with exactly what you're trying to do and the
full model? traceplot only works for quantities declared as parameters,
transformed parameters, or generated quantities.

- Bob

[Linas Mockus]

Thanks, Bob.

On 3/31/2014 5:25 AM, Bob Carpenter wrote:
> On Mar 29, 2014, at 9:45 PM, Linas Mockus <linasm...@gmail.com> wrote:
>
>> I found the solution for 2-compartmental model in WinBUGS:
>> for (t in 1 : N) {
>> model_c <- B1 * exp (-lambda1 * time [t]) +
>> B2 * exp (-lambda2 * time [t]) +
>> B3 * exp (-ka * time [t]);
>> sigma <- sqrt (sigma_m * sigma_m * model_c * model_c +
>> sigma_a * sigma_a);
>> concentration [t] ~ normal (model_c, sigma);
>> }
>>
>> However for some combination of parameters model_c is negative as well as concentration. When using explicit likelihood function I assign very low number to the likelihood whenever model_c<0 - don't know if it is a good practice but basically I tell the solver not to explore infeasible set of parameters. How to do this in stan?
> Priors on B1, B2, B3 would be the way to go, or even hard constraints
> making them positive if that makes sense for the model. Or equivalently,
> putting them inside the exp() function as additive terms.
>
> There's a square(x) function, by the way, which is a bit more
> efficient than x * x. (And we really need to add normal distributions
> parameterized with variance and precision, which would avoid that
> sqrt operation, too --- it's just going to get squared again inside
> the density!).
>
> You could also constrain concentration to be positive, though you'll
> still get odd behavior with model_c going below 0 (it's symmetric around
> 0 in the sigma term).

> B1, B2, B3 are functions of ka, lambda1, lambda2. Therefore the only
> way I see is to constrain the concentration. If I put hard constraint
> real <lower=0> B then I get:
>
>
>
> Invalid value of B2: Error in function validate transformed params
> N4stan5agrad3varE: B2 is -466692.48703089252:0, but must be greater than
> 0Rejecting proposed initial value with zero density.
>
>
>
> How to constrain concentration to be positive? For the time being I
> will assume that sigma_a=0 although I saw WinBUGS examples which have both multiplicative and additive noise.

>> Also, I generate predictive posterior in generated quantities block:
>> for (t in 1 : N_pred) {
>> model_c <- B1 * exp (-lambda1 * time_pred [t]) +
>> B2 * exp (-lambda2 * time_pred [t]) +
>> B3 * exp (-ka * time_pred [t]);
>> sigma <- sqrt (sigma_m * sigma_m * model_c * model_c +
>> sigma_a * sigma_a);
>> concentration_pred [t] <- normal_rng (model_c, sigma);
>> }
>> For some reason traceplot gives for stanfit object:
>> Error in mcmc.list(x) : Arguments must be mcmc objects
> Could you write back with exactly what you're trying to do and the
> full model? traceplot only works for quantities declared as parameters,
> transformed parameters, or generated quantities.

I am still working on traceplot. I found a way how to plot my own
traces using fit@sim$samples but traceplot is more convenient.
Traceplot stopped working when I upgraded to the latest stan version -
so perhaps generated quantities is not a culprit.

Linas
>