Hello! I’m having a lot of trouble fitting a model that (on the surface) should be pretty easy to fit…
To make a long story short, I began with Antia et al.'s (1994) model of within host population dynamics, which looks like this:
We can simulate this model deterministically, as follows:
Antia_Model <- function(t, y, p1){
r <- p1[1]; k <- p1[2]; p <- p1[3]; o <- p1[4]
P <- y[1]; I <- y[2]
# r = Replication rate of parasite (0.1-10)
# p = Max. growth rate of immune system (1.0)
# k = Rate of destruction of parasite by immune system (10^(-3))
# o = Parasite density at which growth rate of IS is half its max. (10^3)
dP = r*P - k*P*I
dI = p*I*(P/(P + o))
list(c(dP, dI))
}
r <- 0.2; k <- 0.001; p <- 1; o <- 1000
parms <- c(r, k, p, o)
P0 <- 1; I0 <- 1
N0 <- c(P0, I0)
TT <- seq(0, 50, 0.1)
results <- lsoda(N0, TT, Antia_Model, parms, verbose = TRUE)
P <- results[,2]; I <- results[,3];
Which gives us this:
And we can generate demographically stochastic data as follows:
x0 <- c(P = 1, I = 1)
a <- c("P*r",
"k*P*I",
"p*I*(P/(P + o))")
nu <- matrix(c(+1, -1, 0,
0, 0, +1), nrow = 2, byrow = TRUE)
r <- 0.2; k <- 0.001; p <- 1; o <- 1000
parms1 <- c(r = r, k = k, p = p, o = o)
tf = 50
method <- "OTL"
simName <- "Antia"
set.seed(6)
Antia_SLC <- suppressWarnings(ssa(x0, a, nu, parms1, tf, method, simName,
verbose = FALSE,
consoleInterval = 1,
censusInterval = 1,
maxWallTime = 30,
ignoreNegativeState = TRUE))
Antia_Sim <- Antia_SLC$data
Antia_Sim <- as.data.frame(Antia_Sim)
colnames(Antia_Sim) <- c("Time", "Parasites", "ImmuneCells")
Which gives us this:
The gist is, I’m trying to fit the model to the demographically stochastic data, and running into unforeseen problems.
Here’s the model I’m using to do the fitting (I have incredibly tight priors for troubleshooting reasons…):
functions {
real[] dz_dt(real t,
real[] z,
real[] theta,
real[] x_r,
int[] x_i) {
real P = z[1];
real I = z[2];
real r = theta[1];
real k = theta[2];
real p = theta[3];
real o = theta[4];
real dP_dt = r*P - k*P*I;
real dI_dt = p*I*(P/(P + o));
return { dP_dt, dI_dt };
}
}
data {
int<lower = 0> N;
real ts[N];
real y_init[2];
real<lower = 0> y[N, 2];
}
parameters {
real<lower = 0> r;
real<lower = 0, upper = 1> k;
real<lower = 0> p;
real<lower = 0> o;
real<lower = 0> z_init[2];
real<lower = 0> sigma[2];
}
transformed parameters {
real z[N, 2]
= integrate_ode_rk45(dz_dt, z_init, 0, ts, {r, k, p, o},
rep_array(0.0, 0), rep_array(0, 0),
1e-5, 1e-3, 5e2);
}
model {
r ~ uniform(0, 0.5); // r = 0.2
k ~ uniform(0, 0.1); // k = 0.001
p ~ uniform(0, 2); // p = 1
o ~ uniform(990, 1010); // o = 1000
sigma ~ lognormal(-1, 1);
z_init ~ normal(1, 1);
for (m in 1:2) {
y_init[m] ~ lognormal(log(z_init[m]), sigma[m]);
y[ , m] ~ lognormal(log(z[, m]), sigma[m]);
}
}
generated quantities {
real y_init_rep[2];
real y_rep[N, 2];
for (m in 1:2) {
y_init_rep[m] = lognormal_rng(log(z_init[m]), sigma[m]);
for (n in 1:N)
y_rep[n, m] = lognormal_rng(log(z[n, m]), sigma[m]);
}
}
And here’s what I’m doing to implement the model:
model <- stan_model(here("./03_Fit/Antia_Fit.stan")) # Which is just where I keep the above model
N <- length(Antia_Sim$Time) - 1
ts <- 1:N
y <- as.matrix(InfecData[2:(N + 1), 2:3])
y <- cbind(y[ , 1], y[ , 2]);
Data <- list(N = N, ts = ts, y_init = y_init, y = y)
fit <- sampling(model, data = Data, chains = 2, iter = 1000, cores = 2, seed = 1)
When I try to run the fitting line, both chains pretty immediately fail to initialise:
Chain 2: Rejecting initial value:
Chain 2: Log probability evaluates to log(0), i.e. negative infinity.
Chain 2: Stan can't start sampling from this initial value.
Chain 2:
Chain 2: Initialization between (-2, 2) failed after 100 attempts.
Chain 2: Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
Does anyone know what might be causing this error? I know it’s a pretty standard one, but this seems like a pretty straightforward question and I don’t know what I’m doing wrong.
Any help would be incredibly appreciated :)