Here’s a possibly odd question: does anyone know how to supply more than a single set of initial conditions – both for time t0 and the state of the system R0 – for the integrate_ode_rk45() function?
Context: I’m fitting ODEs to empirically observed experimental data, where we model each treatment as a separate ODE; each ODE is assumed to have the same parameter values *, which we want to estimate. We know initial conditions for both state and time, and we observed the system change through time. Easy enough so far.
Here’s where it gets difficult: we (intentionally) “re-initialized” our experimental conditions twice after first beginning the experiment, but we did this without shuffling treatments. In essence, we have a sequence of three temporally correlated experimental trials. The key for why we want to “re-initialize” the ODEs in Stan, and why I put the * asterisk above: we clearly have time-dependency in the parameters we’re wanting to estimate (i.e., the parameter values change across the three experimental trials), and this time-dependency is especially of interest for the biology we’re trying to model.
Here is simplified R code that uses an ODE to simulate data, add noise, and then set up a list to pass to Stan. This first example does NOT include multiple ODE initializations, and this first example works well:
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## ~~~~~~~~~~~ Simulate an ODE to fit in STAN ~~~ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## open ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
rm(list=ls())
library(rstan)
library(StanHeaders)
library(dplyr)
library(deSolve)
library(ggplot2)
library(reshape)
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Simulate ODE ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## params
params <- c(k=.025)
## create a list of initial conditions
R0 <- seq(50, 300, length.out=5)
## specify ODE function
resourceLoss <- function(t, R0, params) {
with(as.list(c(R0, params)),{
dR <- R0 * -k
list(c(dR))
})
}
## set vector of times
t.list <- seq(1, 50, length.out=3)
#t.list <- c(0, 1, 100)
## run ODE
out <- ode(y=c(R0), times=t.list, func=resourceLoss, parms=params)
## save as data frame
loss <- round(as.data.frame(out),3)
## save list without Times for later use
dRdt.list <- round(out[,-1],3)
## create and add prefix for col names
loss.dat <- as.data.frame(dRdt.list)
colnames(loss.dat) <- paste("T", colnames(loss.dat), sep="")
## plot all dRdt
all_dRdt <- melt(loss, id.vars="time")
tplot <- ggplot(all_dRdt, aes(time, value, col=variable,group=variable)) +
geom_point() + geom_line(color="black") +
xlab("Time") + ylab("Density") + ggtitle("Remaining density = R0 - dRdT")
print(tplot)
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## add noise ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## set dims for matrix and set sigma
rws=length(t.list) ## number of rows ie data points to simulate
cls=ncol(dRdt.list) ## number of columns ie separate FR runs
sd=5
set.seed(10)
## apply rgamma to deterministic data
mult_ODE <- function(n, dat, noise){
noise.out <- rnorm(n=length(dRdt.list), mean=dat, sd=sd)
return(matrix(noise.out, nrow=rws, ncol=cls, byrow=FALSE))
}
noisy.dat <- round(mult_ODE(n, dat=dRdt.list, noise=sig),digits=3)
## replace noisy row 1 with deterministic R0
noisy.dat[1,] <- dRdt.list[1,]
## create data frame for plotting, add time as a column
noise.dat.df <- as.data.frame(noisy.dat)
noise.dat.df <- cbind(loss$time, noise.dat.df)
names(noise.dat.df)[1]<-"time"
## plot
melt_noise <- melt(noise.dat.df, id.vars="time")
noise_plot <- ggplot(melt_noise, aes(time, value, col=variable, group=variable)) +
geom_point() + geom_line(data=all_dRdt, aes(time, value, group=variable), col="black") +
xlab("Time") + ylab("Density") + ggtitle("Density remaining = R0 - dRdT + noise")
print(noise_plot)
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## configure data list for STAN ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## extract deterministic list of initial conditions
R0.list <- as.numeric(loss.dat[1,])
## list for Stan
## standard time
resourceLoss.dat <- list(Tm = nrow(loss.dat)-1,
t0 = t.list[1],
th = t.list[-1],
nLev = ncol(noisy.dat),
R0_obs = R0.list,
R_obs = noisy.dat[-1,])
## run STAN
fit1 <- stan(
file = "ExpDecay_StanExample.stan",
data = resourceLoss.dat,
chains = 1,
warmup = 1000,
iter = 2000,
cores = 1,
refresh = 10)
## END standard time ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
And here’s the associated code in Stan, and again, this too works fine:
functions {
real[] resourceLoss(
real t, // time
real[] R, // state
real[] theta, // params
real[] x_r,
int[] x_i) {
real dRdt[x_i[1]];
real k;
k = theta[1];
for (i in 1:x_i[1]) {
dRdt[i] = R[i] * -k;
}
return dRdt;
}
}
data {
int <lower=1> Tm; // Tm = length(t.list - initial t point)
int <lower=1> nLev; // number of different R0, i.e. # of ODEs
real R0_obs[nLev]; // Intial conditions R0 list
real <lower=1> t0; // starting time t = 1
real <lower=t0> th[Tm]; // array of times for ODE, length Tm
real R_obs[Tm, nLev]; // simulated dat
}
transformed data{
real x_r[0];
int x_i[1];
x_i[1] = nLev;
}
parameters {
real <lower=0> k;
real <lower=0> sigma;
}
transformed parameters {
real theta[1];
theta[1] = k;
}
model {
real R[Tm, nLev];
R = integrate_ode_rk45(resourceLoss, R0_obs, t0, th, theta, x_r, x_i);
k ~ normal(0,5);
sigma ~ normal(0,5);
for (t in 1:Tm) {
for (i in 1:nLev) {
R_obs[t,i] ~ normal(R[t,i], sigma);
}
}
}
Here’s where I attempt to pass a vector for t0. Note that this code does not even include the “multiple re-initializations” of the ODE . . . I’m first simply trying to get Stan to accept a vector for t0, after which I’d then expand to include multiple re-initializations (this code replaces the end of the code above, starting at ## configure data list for STAN ~~~)
## treating intial time as a vector
t.list <- seq(1, 50, length.out=3)
t0.list <- rep(t.list[1], ncol(loss.dat))
t.obs_1 <- rep(t.list[2], ncol(loss.dat))
t.obs_2 <- rep(t.list[3], ncol(loss.dat))
## all time points -- partitioned and passed to Stan below
full.t.list <- rbind(t0.list, t.obs_1, t.obs_2)
## extract deterministic list of initial conditions
R0.list <- as.numeric(loss.dat[1,])
## intialize list for Stan
resourceLoss.dat <- list(Tm = nrow(loss.dat)-1,
t0 = full.t.list[1,],
th = full.t.list[-1,],
nLev = ncol(loss.dat),
R0_obs = R0.list,
R_obs = array(loss.dat[-1,]))
## run STAN
fit1 <- stan(
file = "ExpDecay_StanExample.stan",
data = resourceLoss.dat,
chains = 1,
warmup = 1000,
iter = 2000,
cores = 1,
refresh = 10)
Unsurprisingly, I ran into:
Wrong signature for function integrate_ode_rk45; first argument must be the name of a function with signature (real, real[ ], real[ ], data real[ ], data int[ ]): real[ ].
In case it’s at all useful, here’s the conceptual layout of what the initial state + observed values would look like, Where R0_X are the experimentally supplied initial conditions for the state of the system, and DX represent observed data points as the state of the system changes through time; we supply all these values; this represents four treatments and thus 4 ODEs “re-initialized” twice, i.e., three total temporal ODE “runs”:
R0_1, R0_1, R0_1, R0_1
D2, D2, D2, D2,
D3, D3, D3, D3,
R0_2, R0_2, R0_2, R0_2,
D4, D4, D4, D4,
D5, D5, D5, D5,
R0_3, R0_3, R0_3, R0_3
D6, D6, D6, D6,
D7, D7, D7, D7,
And where the first initializing time point is t0_1 = 1, the second initialization of the system is t0_2 = 4, and the third is t0_3 = 7. (what I haven’t included here is a vector of data with “cumulative time”, and it is this vector that will be the predictor for the parameter’s temporal dependency).
t0_1, t0_1, t0_1, t0_1,
2, 2, 2, 2,
3, 3, 3, 3,
t0_2, t0_2, t0_2, t0_2,
5, 5, 5, 5,
6, 6, 6, 6,
t0_3, t0_3, t0_3, t0_3,
8, 8, 8, 8,
9, 9, 9, 9
In essence, we want to string together three separate initializations of an ODE system. And it seems one way to do this is to pass a vector to t0 in Stan. Is this at all possible? If it’s not possible in the context of integrate_ode_rk45(), is there some other context that might work?
Thank you all in advance!!