Hi all,

I’ve been trying to run an ODE model on Stan, but have been facing issues with the fit. I know that the y[3] output from the ode solver should have a peak within a certain time period. But because the data I have is obtained after this peak, I have been having trouble estimating beta (params[1]) and incubation period (inc_per). All my fits have either been ‘pushing’ inc_per to the upper limit accompanied by a low median value of beta, or inc_per has been way too low (close to or at 0 days) and median value of beta way too high and hitting its upper limit (causes peak to occur at time 0.1 days)

I was wondering if anyone had a suggestion on how to solve this issue. Or if there was a way to restrict the ode solver such that the peak for y[3] has to occur within a specific window of time. I think making inc_per an individual-specific parameter might solve this issue, but the dataset I’m using has over 30k individuals, so introducing an individual-specific level to the hierarchy wouldn’t be wise.

Thank you!

```
functions {
real[] OneSys(real t, real[] y, real[] params, real[] x_r, int[] x_i) {
real dydt[4];
dydt[1] = 0.14*( ((8*10^7)/exp(y[1])) - 1 ) - params[1] * exp(y[3]);
dydt[2] = (params[1] * exp(y[1]) * exp(y[3]))/exp(y[2]) - 1.07;
dydt[3] = (50 * exp(y[2]))/exp(y[3]) - params[2] - 0.001 * exp(y[4]);
dydt[4] = params[3] * exp(y[2]);
return dydt;
}
}
data {
int<lower = 0> L; //no. of rows
int<lower = 0> G; //no. of groups
int<lower = 0> AG; //total no. of agegr
int l_tsolve; //length of vector of times to solve
real logvirus[L]; //data
real t0; //initial value for t
int ts[L]; //all days (arranged by individual)
real t_solve[l_tsolve]; //vector of times to solve
int pos_vgofag[AG]; //vector of which vg each ag is in, indexing for tinc
int count_agegr[AG]; //vector of no. of rows in each agegr
int pos_agegr[AG]; //vector of indices of first row for each agegr
int count_vacgr[G]; //vector of no. of rows in each vacgr
int pos_vacgr[G]; //vector of indices of first row for each vacgr
real rel_tol;
real abs_tol;
int max_num_steps;
int which_fixed; //which parameter to fix as reference group
real fixed_value; //value to fix reference group at
}
transformed data {
real x_r[0];
int x_i[0];
}
parameters {
real<lower = 10^-7, upper = 5*10^-6> beta_par;
real<lower = 0.1, upper = 20> c_par;
real<lower = 10^-13, upper = 1> omega_par[G];
real<lower = 0, upper = 10> phi;
real<lower = 0, upper = 15> lambda;
real<lower = 1, upper = 10> inc_per;
real<lower = 0.999999, upper = 10.00000> imm_resp_mult[G-1];
}
transformed parameters {
real predVal[L];
real params[3];
real y0[4];
real<lower = 0.1, upper = 1> imm_resp[G];
{
real t_inc[l_tsolve];
for (l in 1:l_tsolve) {
t_inc[l] = t_solve[l] + inc_per;
}
y0[1] = log(8*10^7); //T0
y0[2] = log(10^-6); //I0
y0[3] = log(10); //V0
params[1] = beta_par;
params[2] = c_par;
imm_resp[which_fixed] = fixed_value;
for (vg in 2:G) {
imm_resp[vg] = imm_resp[which_fixed] * imm_resp_mult[vg-1];
}
// iterate over total age groups
for (ag in 1:AG) {
real yout[l_tsolve, 4];
y0[4] = log(imm_resp[pos_vgofag[ag]]);
params[3] = omega_par[pos_vgofag[ag]];
// params[3] = omega_par[ag];
yout = integrate_ode_bdf(OneSys, y0, t0, t_inc, params, x_r, x_i, rel_tol, abs_tol, max_num_steps);
predVal[pos_agegr[ag] : (pos_agegr[ag] + count_agegr[ag] - 1)] = yout[ts_ind[pos_agegr[ag] : (pos_agegr[ag] + count_agegr[ag] - 1)], 3];
}
}
}
model {
omega_par ~ beta(phi, lambda);
inc_per ~ normal(3,2);
// imm_resp ~ beta(1,1);
for (i in 1:L) {
logvirus[i] ~ normal(predVal[i], 1);
}
}
```