Hello all,
My goal is to fit an ODE model with 3 state variables to many days of hourly data. I have been working in R and using an L-BFGS-B algorithm with the “optim” function, but it occurred to me that a hierarchical MCMC approach would be better. The idea is that each day can have varying parameters (ideally smoothly varying) as shown with this nice example (in Hierarchical Models section at the end): Fitting differential equation models with Stan. Each day also has two time series of forcing inputs (light and temperature).
I found this paper that seems to be a good analogy of what I want to do (Supplement with Stan code here): https://ars.els-cdn.com/content/image/1-s2.0-S0378429022001204-mmc1.pdf
However, it appears to be set up quite differently than what is shown here in Section 4.3, for example: Differential Equations Based Models in Stan
I am entirely new to writing code in Stan. I’ve read the intro docs about the blocks and what goes in them, but my survey of examples indicates that there is more than one way to go about coding this. One of my main hangups is how to code it hierarchically. Do I use matrix format in the parameter block like matrix[T,N] or do I use for loops to loop through each day?
Would someone please advise me how to best code this? I’ve attached my initial effort below. Do I put priors in the data block or the model block…or both? I’ve seen both.
functions {
real[] myODE(real t, real[] y, real[] theta, real[] x_r, int[] x_i) {
real dydt[3]; // Number of state variables
// Unpack parameters and states
real VAR1 = y[1];
real VAR2 = y[2];
real VAR3 = y[3];
real A = theta[1];
real B = theta[2];
real C = theta[3];
real temp_forcing = x_i[1] == 1; // not actually sure how to implement forcing...
real light_forcing = x_i[2] == 1;
// Define other constants and variables
// ... (continue defining other constants and variables)
// Calculate derivatives
dydt[1] = ((A - B + fVAR1 + LVAR1) / z); // Derivative of VAR1
dydt[2] = ((-A + B + fVAR2 + LVAR2) / z); // Derivative of VAR2
dydt[3] = 2 * C + LVAR3 / z; // Derivative of VAR3
return dydt;
}
}
data {
int<lower=1> T; // Number of time points per day
int<lower=1> N; // Number of days (individuals/groups)
real<lower=0> timestep; // length of each timestep in days
int<lower=1> n; // number of observations per day
matrix[T, N] y_VAR1; // Observed VAR1 (T rows for hours, N columns for days)
matrix[T, N] y_VAR2; // Observed VAR2 (T rows for hours, N columns for days)
matrix[T, N] y_VAR3; // Observed VAR2 (T rows for hours, N columns for days)
vector[T] light; // Time-varying light data
vector[T] temp; // Time-varying temp data
real y_init[3]; // Initial values for states
int<lower=0, upper=1> temp_forcing[N]; // Binary indicator for temperature forcing
int<lower=0, upper=1> light_forcing[N]; // Binary indicator for light forcing
// ... (other data variables)
}
parameters {
real<lower=0> theta[1]; // parameter(s)
real<lower=0> theta[2]; // parameter(s)
real<lower=0> theta[3]; // parameter(s)
\\ or something like matrix[T, N] theta[1] ??
real<lower=0> sigma; // Measurement error
}
model {
real y_hat[N, 3]; // Predicted states
real y_init_local[3] = y_init; // Initialize states with provided initial values
// Loop through time points
for (n in 1:N) {
int x_i[2];
x_i[1] = temp_forcing[n];
x_i[2] = light_forcing[n];
y_hat[n] = integrate_ode_rk45(myODE, y_init_local, 0, {theta}, t[n], x_i, x_i);
// Likelihood
y[n, 1] ~ normal(y_hat[n, 1], sigma); // VAR1
y[n, 2] ~ normal(y_hat[n, 2], sigma); // VAR2
y[n, 3] ~ normal(y_hat[n, 3], sigma); // VAR3
// Update initial values for the next time step
y_init_local = y_hat[n];
}
}