Forecast extension to Lotka Volterra Model

Modeling
, ,
1

I have spent some time looking at the Lotka-Volterra example, and I tried to build some of the Exercices and Extensions (btw it is a great example for beginners).

It was straightforward to extend the predictions into the future by simply changing the time parameter of the ODE. However, the result is not really satisfying because the confidence interval does not become wider over time. I feel that the further away we get from the data we should have less confidence in our forecats.

Is there any approach to incorporate this uncertainty into the model?

Here is my stan code, there is not much changed relative to the original model:

// Lotka Volterra Model

functions {
  real[] dz_dt(real t,       // time
               real[] z,     // system state {prey, predator}
               real[] theta, // parameters
               real[] x_r,   // unused data
               int[] x_i) {
    real u = z[1];
    real v = z[2];

    real alpha = theta[1];
    real beta = theta[2];
    real gamma = theta[3];
    real delta = theta[4];

    real du_dt = (alpha - beta * v) * u;
    real dv_dt = (-gamma + delta * u) * v;
    return { du_dt, dv_dt };
  }
}


data {
  int<lower = 0> N;           // number of measurement times
  int<lower = 0> Nf;          // number of periods to forecast
  real ts[N+Nf];              // measurement and forecast times > 0
  real y_init[2];             // initial measured populations
  real<lower = 0> y[N, 2];    // measured populations
}
parameters {
  real<lower = 0> theta[4];   // { alpha, beta, gamma, delta }
  real<lower = 0> z_init[2];  // initial population
  real<lower = 0> sigma[2];   // measurement errors
}
transformed parameters {
  real z[N+Nf, 2] = integrate_ode_rk45(dz_dt, z_init, 0, ts, theta,
                         rep_array(0.0, 0), rep_array(0, 0),
                         1e-5, 1e-3, 5e2);
}
model {
  theta[{1, 3}] ~ normal(1, 0.5);
  theta[{2, 4}] ~ normal(0.05, 0.05);
  sigma ~ lognormal(-1, 1);
  z_init ~ lognormal(log(10), 1);
  for (k in 1:2) {
    y_init[k] ~ lognormal(log(z_init[k]), sigma[k]);
    y[ , k] ~ lognormal(log(z[1:N, k]), sigma[k]);
  }
}
generated quantities {
  real y_init_rep[2];
  real y_rep[N+Nf, 2];
  for (k in 1:2) {
    y_init_rep[k] = lognormal_rng(log(z_init[k]), sigma[k]);
    for (n in 1:(N+Nf))
      y_rep[n, k] = lognormal_rng(log(z[n, k]), sigma[k]);
  }


}

And here is my R code:


N <- length(lynx_hare_df$Year) - 1

y_init <- c(lynx_hare_df$Hare[1], lynx_hare_df$Lynx[1])
y <- as.matrix(lynx_hare_df[2:(N + 1), 2:3])
y <- cbind(y[ , 2], y[ , 1]); # hare, lynx order

Nf <- 50  # years to forecast
ts2 <- 1:(N+Nf) # times

lynx_hare_data_fcast <- list(N = N, Nf = Nf, ts = ts2, y_init = y_init, y = y)

model_fcast <- stan_model("lotka-volterra-forecast.stan")

fit_fcast <- sampling(model_fcast, data = lynx_hare_data_fcast, seed = 123)

print(fit_fcast, pars=c("theta", "sigma", "z_init"),
      probs=c(0.1, 0.5, 0.9), digits = 3)

z_init_draws <- extract(fit_fcast)$z_init
z_draws <- extract(fit_fcast)$z
y_init_rep_draws <- extract(fit_fcast)$y_init_rep
y_rep_draws <- extract(fit_fcast)$y_rep
predicted_pelts <- matrix(NA, 1+N+Nf, 2)
min_pelts <- matrix(NA, 1+N+Nf, 2)
max_pelts <- matrix(NA, 1+N+Nf, 2)
for (k in 1:2) {
  predicted_pelts[1, k] <- mean(y_init_rep_draws[ , k])
  min_pelts[1, k] <- quantile(y_init_rep_draws[ , k], 0.25)
  max_pelts[1, k] <- quantile(y_init_rep_draws[ , k], 0.75)
  for (n in 2:(1+N+Nf)) {
    predicted_pelts[n, k] <- mean(y_rep_draws[ , n - 1, k])
    min_pelts[n, k] <- quantile(y_rep_draws[ , n - 1, k], 0.25)
    max_pelts[n, k] <- quantile(y_rep_draws[ , n - 1, k], 0.75)
  }
}

lynx_hare_melted_df <- melt(as.matrix(lynx_hare_df[, 2:3]))
colnames(lynx_hare_melted_df) <- c("year", "species", "pelts")
lynx_hare_melted_df$year <-
  lynx_hare_melted_df$year +
  rep(1899, length(lynx_hare_melted_df$year))

Nmelt <- dim(lynx_hare_melted_df)[1]
lynx_hare_observe_df <- lynx_hare_melted_df
lynx_hare_observe_df$source <- rep("measurement", Nmelt)

lynx_hare_predict_df <-
  data.frame(year = rep(1900:(1900+N+Nf), 2),
             species = c(rep("Lynx", 1+N+Nf), rep("Hare", 1+N+Nf)),
             pelts = c(predicted_pelts[, 2],
                       predicted_pelts[, 1]),
             min_pelts = c(min_pelts[, 2], min_pelts[, 1]),
             max_pelts = c(max_pelts[, 2], max_pelts[, 1]),
             source = rep("prediction", 2*(1+N+Nf)))

lynx_hare_observe_df$min_pelts = NA #lynx_hare_predict_df$min_pelts
lynx_hare_observe_df$max_pelts = NA #lynx_hare_predict_df$max_pelts

lynx_hare_observe_predict_df <-
  rbind(lynx_hare_observe_df, lynx_hare_predict_df)

population_plot2 <-
  ggplot(data = lynx_hare_observe_predict_df,
         aes(x = year, y = pelts, color = source)) +
  geom_vline(xintercept = 1900, color = "grey") +
  geom_hline(yintercept = 0, color = "grey") +
  facet_wrap( ~ species, ncol = 1) +
  geom_ribbon(data = lynx_hare_observe_predict_df %>% filter( source=='prediction' ),
              aes(ymin = min_pelts, ymax = max_pelts),
	      colour = NA, fill = "black", alpha = 0.1) +
  geom_line() +
  geom_point() +
  ylab("pelts (thousands)") +
  ggtheme_tufte() +
  theme(legend.position="bottom") + labs(x=NULL)

population_plot2
1 Like
3

Hi!

I think one assumption behind such a model is the stability of the attractor.

So if you know the parameters exactly and assume that everything will remain constant, there is no inherent reason to consider that longer term predictions should be less accurate than shorter ones.

For longer term predictions to be less accurate, you should integrate the possibility that the system will be disturbed, or put a long-term pressure on it.

I am not sure if I described my idea clearly enough, tell me :)

Lucas

2 Likes
4

ode’s are deterministic – given the model parameters, the system is perfectly predictable up to some constant due to measurement error. stochastic differential equations, or some form of discretized system that includes uncertainty, are often more realistic in that they acknowledge there is unknown structure in the system. There’s a lot of ways to do it, I don’t think any are too simple in stan code, but maybe someone else has good ideas. It’s a nice example, making it stochastic would be interesting, I might try post it as a ctsem ( https://cdriver.rbind.io/post/ctsem-quick-start/ ) example. This would also give some stan code, just very ugly and overcooked :)

1 Like
5

Thank you for your insight, I understand the intuition of why an ODE based model is not the right aproach for forecasting if we want to allow for perturbabtions over time.

To explain where I am coming from, I tried to fit an SIR model (which I guess is a fashionable thing to do these days), and I observed the same issue: the ODE approach does not allow for any uncertainty in the system, and the error interval does not grow over time

I would be curious to see if anyone has any suggestions of an alternative approach.

1 Like
6

Lotka Volterra case study, excercises - Predator-Prey Population Dynamics: the Lotka-Volterra model in Stan

Forecasting and backcasting . Extend predictions another 50 years into the future and plot as in the last plot. This can be done by extending the solution points in the transformed parameters, but is more efficiently done in the generated quantities block.

to roll the uncertainty in your estimates into the future, you use the y_rep values in the generated quantities block as follows, assuming additional data block variable Nf which is the number periods into the future to do the forecast:

generated quantities {
  real y_init_rep[2];
  real y_rep[N, 2];
  real z_forecast[Nf, 2];
  real y_forecast[Nf, 2];
  for (k in 1:2) {
    y_init_rep[k] = lognormal_rng(log(z_init[k]), sigma[k]);
    for (n in 1:N)
      y_rep[n, k] = lognormal_rng(log(z[n, k]), sigma[k]);
  }
  z_forecast = integrate_ode_rk45(dz_dt, y_rep[N], 0, ts_forecast, theta,
                         rep_array(0.0, 0), rep_array(0, 0),
                         1e-5, 1e-3, 5e2);
  for (k in 1:2) {
    for (n in 1:Nf)
      y_forecast[n, k] = lognormal_rng(log(z_forecast[n, k]), sigma[k]);
    
  }
}

the resulting predictions, while perhaps not good ones, do have greater uncertainty - grey marks central 50% of the credible interval:

Rplot01
Rplot011000×536 65.4 KB

1 Like
7

not sure if the above is what was intended - when computing z_forecast, should one use y_rep[N] for the initial state, as above, or z[N], i.e.:

  z_forecast = integrate_ode_rk45(dz_dt, z[N], 0, ts_forecast, theta, rep_array(0.0, 0), rep_array(0, 0), 1e-5, 1e-3, 5e2);

for the latter, forecasts for 50 and 100 years look like this:

lv_forecast_bb
lv_forecast_bb1000×536 68.2 KB

ls_forecast_bb2
ls_forecast_bb21000×301 50.9 KB

time to make a lifeline call: @Bob_Carpenter which is correct?

1 Like
8

the simple answer is yes, and there’s a new chapter in the Stan User’s Guide on this: Stan User’s Guide

specifically, section sampling from the posterior predictive distribution discusses the sources of uncertainty:

Randomly drawing \tilde{y} from the sampling distribution is critical because there are two forms of uncertainty in posterior predictive quantities, sampling uncertainty and estimation
uncertainty. Estimation uncertainty arises because \theta is being estimated based only on a sample of data y. Sampling uncertainty arises because even a known value of \theta leads to a sampling distribution p(\tilde{y} \mid \theta) with variation in \tilde{y}.
Both forms of uncertainty show up in the factored form of the posterior predictive distribution,

p(\tilde{y} \mid y) = \int \underbrace{p(\tilde{y} \mid \theta)}_{ \textrm{sampling uncertainty}} \cdot \underbrace{p(\theta \mid y)}_{ \textrm{estimation uncertainty}} \textrm{d}\theta

getting back to above comment, for the Lotka-Volterra case study, forecasts which take z[N] as the initial value of y when computing z_forecast don’t roll sampling uncertainty into the forecast, samples which take y_rep[N] as the initial value of y do - as seen in the above plots.

3 Likes
9

Thank you for following up, and indeed the code in my original question is equivalent to using z[N] as the initial value for the forecast (I just implemented it by simply calling the integrate_ode_rk45 for a longer period in the transformed parameters call).

FWIW, I think you are right, integrating the ODE again while rolling the uncertainty into the new previous final state gives a more realistic model.