Running hierarchical model consisting of ODE

Modeling
1

Hi all,

This is my first time using Stan so forgive any dumb errors. I am trying to fit a hierarchical/random effects model consisting of a system of first order differential equation. One of the parameters (epsilon) in the ODE

comes from a gamma distribution and is unique for each individual. So (r,epsilon_i,K) for i = 1,…,J are ODE model parameters, alpha and beta are hyperparameters and sigma is the observance variance. The data is one observation for individual. Here’s my attempt in pystan:

obdata = {‘J’: ncats,‘ts’: daypool, ‘y1_data’ : countpool}
withinhostcode = “”"
functions {
real[] whmodel(real t, // time
real[] y, // state
real[] theta, // parameters
real[] x_r, // data (real)
int[] x_i) { // data (integer)
real dydt[2];
dydt[1] = theta[1] * y[1] * (1 - (y[1]/theta[2])) - theta[3] * y[1] * y[2];
dydt[2] = -theta[3] * y[1] * y[2];
return dydt;
}
}

data {
int<lower=1> J;
real ts[J,3];
real y1_data[J];
}

transformed data {
real y0[2];
real t0;
real x_r[0];
int x_i[0];

t0 = 0;
y0[1] = 10;
y0[2] = 100;
}

parameters {
real<lower=0> r;
real<lower=0> K;
real<lower=0> eps[J];
real<lower=0> alpha;
real<lower=0> beta;
real<lower=0> sigma;
}

transformed parameters{
real y[3,2];
real y_hat[J];
real theta[3];

for (i in 1:J)
{

theta[1] = r;
theta[2] = K;
theta[3] = eps[i];

y = integrate_ode_bdf(whmodel, y0, t0, ts[i], theta, x_r, x_i);
y_hat[i] = y[3,1]; //solution corresponding to observed data time

}
}

model{
eps ~ gamma(alpha, beta); //alpha, beta have uniform priors
r ~ gamma(100,0.00292);
K ~ normal(2000,2000);
sigma ~ normal(0,1000);
y1_data ~ normal(y_hat, sigma); //Gaussian error assumption
}
"""

However, when I run the model using
sm = StanModel(model_code = withinhostcode, verbose=True)
fit = sm.sampling(data = obdata, iter =100, chains = 1)

It successfully compiles but once it gets to fit, doesn’t really print anything or exit. It gets stuck even when I do a few iterations. Not quite sure if I’m doing anything wrong or it’s really that slow to sample from. Any suggestions would be highly appreciated!

2

Quickly going over your model I would say

  • make your priors weakly-informative. They are way too vague. You should simulate from the prior and check if the simulated ranges of observed data is at the same order of magnitude

  • make use of the advanced version of the integrate_ode_bdf function call which allows you to set relative, absolute tolerance and maximal number of integration steps. All of these are problem specific; in adequatly scaled problems a 10^-8 relative, 10^-5 absoulte and 1e3 maximal number of steps is often reasonable. These settings assume that your problem has mean values of at least O(1) roughly.

Best,
Sebastian

3

Thank you! I’ll try some of your suggestions. Does that mean you don’t see anything glaringly wrong in the code (eg syntax wise)?

1 Like
4

Does that mean you don’t see anything glaringly wrong in the code (eg syntax wise)?

Syntax is up to the compiler, which if it compiles, it likes :D.

What you describe sounds like Stan is exploring a bit of parameter space that is very hard on the ODE solver. With ODEs stuff, looking at crazy parameters can be a lot tougher than for models with closed form solutions cause running ODE solvers with crazy parameters at reasonable tolerances can be really difficult (take a long time to solve at the req’d tolerances, which is probably what you’re seeing).

That’s where @wds15’s first suggestion comes from (“You should simulate from the prior and check if the simulated ranges of observed data is at the same order of magnitude”)

Hopefully there’s a way to set up your priors to avoid these crazy bits of parameter spaces. If there’s nothing telling it not to try, the sampler will happily head off and explore these weird places.

r ~ gamma(100,0.00292);
K ~ normal(2000,2000);

And uniform [0, inf) on alpha and beta are super vague. The initial conditions on [0, inf) can be set very large I think, so especially if you think your parameter is less than millions (or some other number less than infinity), add a tighter prior :D.

I don’t know your problem or the scales of the numbers in them, but if I thought a parameter was somewhere in the range 1-5, I’d just slap a normal(0, 5) prior on it and see what happens from there. This is still a pretty weak prior, but if you’re scared this is biasing your results, you can always weaken the prior later and see what happens.

5

Thank you :) I’ll keep playing around with it.