Hi,
Can anybody share experiences of fitting in stan time series (dissolution of a drug) data to Weibull function
concentration=m0 * (1 - exp(-((time - tlag) / scale) ^ shape)).
I have problems handling tlag parameter. However, when tlag=0 convergence is good.
Thanks a lot,
Linas
doasimcyp.stan (2.1 KB)
What sorta problems? Looks like there’s a lot more going on in this model than just y ~ normal(weibull(...), sigma).
Are you getting divergences? Or crazy looking parameter estimates?
Is tlag = 0 a reasonable estimate for tlag?
In a complicated model like this it’s always worth verifying that you can estimate the parameters of data you generate with your model. And maybe if possible try a simple y ~ normal(weibull(...), sigma) regression? I dunno if that makes sense in this case – if you can take a single time series and get super noisy estimates or whatever.
Well, it doesn’t converge. However, when I assume tlag=0 by removing
this parameter everything converges to sensible estimates.
real<lower=0,upper=1> tlag_raw;
If tlag == 0 produces gives reasonable estimates, shouldn’t the tlag parameter be allowed to go down to zero? I see there’s the multivariate normal stuff creeping in through d[4], but I don’t get this constraint still.
Does doing something like:
real tlag[NLOTS]; // Is NLOTS the right thing? is this how many curves you have?
And:
tlag ~ normal(0, 0.1) // (tight prior around a value you know works?)
Do anything? Or does that not make sense?
Yes, there is a random effect d. The random effect is on unit level -
each unit may have different tlag, and on lot level. Lots are comprised
of 12 units. If we start dissolving the unit it doesn’t start dissolving
until tlag. Thus tlag>=0. tlag_raw is a fixed effect. Thus real
constraint is tlag_raw+d[4]>=0.
Hmm, I see what you’re saying. It makes sense. But for lack of plots or data, what I’m trying to do is think of a small change in the model that might behave differently enough to give insight into the where the problem is :D.
I’m not really sure what you mean by not converging honestly, so hard to give specific advice. Do you have pairplots? Posterior distributions on the bad tlag? Is it that different chains aren’t converging? Or a multimodality? Or just lots of divergences?
If you can post a runnable model w/ data and stuff I wouldn’t mind running it.
Thanks for looking into it. I ran 4 chains and they don’t mix.
when tlag=0 I get a good mixing
doasimcyp.stan (2.67 KB)
doasimcyp.stan.data.r (38.4 KB)
At this point I’m probably gonna get around to trying this today, but I should be able to take a swing at it tomorrow ^^.
I worked on this awhile but didn’t get too far. I one issue might be with the if(time > tlag) conditional.
I don’t think the derivatives of the Weibull function are always zero at zero (going from pictures here: https://en.wikipedia.org/wiki/Weibull_distribution). Functions we put in HMC have to be continuous and have first derivatives, otherwise it can mess up the sampler, so it might not be totally okay to put a zero on the left side of the Weibull.
If we’re shifting the timescale, that’s a problem. It’s not a problem for the other parameters cause evaluations of that function either stay on one side or the other.
The most interesting thing I got was moving the tlag away from the shape/scale parameters:
y = m0 * (1 - exp(-(time / scale) ^ shape) - tlag_not_quite);
If it sits on the outside like this, the parameters seem to sample a lot better than if it sits inside with the time (you can verify this with a smaller model that only works on one lot if you want).
If you parameterize this way, you could do algebra to compute effectively what the t0 would be to reproduce the effective argument to the exp, but I dunno how far that gets you. The priors on tlag_not_quite wouldn’t be that easy to interpret, but if you’re just wanting to fit the curves, might be worth a try.
Or I guess add a constraint so that shape > 1.0, cause in those cases the derivative of the Weibull function wrt the parameter is actually zero at zero (so matches with the constant function zero). But maybe push it a bit away from 1.0 with the prior too. That transition looks really rough :/.
Well, I also felt that conditional if (time>tlag) creates problems. By the
way, derivative of Weibull function (not the distribution) is not defined
at time=tlag, and is positive for time>tlag, and 0 for time<tlag.
Just a random thought - if we go non-physical and make Weibull function
symmetric (at tlag) so it is negative at time<tlag, = 0 at time=tlag, >0
when time>tlag, and derivative is positive. Technically Weibull function is
always positive but maybe for modeling in stan we don’t have to adhere to
physics.
Hmm, I think the Weibull function here has a derivative at zero (of zero) if the shape parameter is greater than one (from the eqs. it looks like the CDF of the Weibull distribution, which is why I was looking at the plots on Wikipedia).
I don’t think symmetric helps here, cause the only way you don’t have discontinuous derivatives of a symmetric function at zero (or time = tlag) is if the derivative is zero (otherwise things kinda look like absolute value), which is what we want.
Well, for dissolution shape parameter is <1 so the derivative doesn’t
exist at t=tlag. I wonder if it makes sense (is it lawful) to extend
Weibull function to t<=tlag so that derivative always exist? Only caveat
that values t<tlag will be negative (non-physical). The similar problem
must occur for PK models with lag time. I will check how Torsten guys
did it in stan.
Aaah, I see. Yeah you could probably chop of the function at some point t > tlag and then just match the values and first derivative with a tail function that looks something like a * exp(-b * (t - tlag)). Doesn’t seem super elegant :P.
I will check how Torsten guys did it in stan.
Yeah, seems like a good plan if they’re fitting these things.