Solution_predation_jhelam.Rmd

---
title: "Solution: predator-prey lotka-volterra"
author: "Jhelam N. Deshpande"
date: "02/04/2023"
output: html_document
---
# Import and visualise the data set
```{r}
rm(list=ls()) # clear workspace
# import required libraries
library(rstan)
library(coda)
library(deSolve)
```
```{r}
#read data set
data<-read.csv("Data/predator_prey.csv")
#view data set
head(data)
```
```{r}
#plot time series
plot(data$time,data$n,col="blue",xlim=c(0,100),ylim=c(0,50),ylab="N1 or N2",xlab="Time",pch=16)
points(data$time,data$p,col="red",pch=16)
```


# Formulate the model

We fit the following system of ODE to the data. The N(t) represent the dynamics of the prey and P(t) the predator. The following equation represents prey dynamics:
$$\frac{dN}{dt}=rN-aNP$$
And predator dynamics are given by:
$$\frac{dP}{dt}=eaNP-dP$$
$r$ is prey growth rate, $a$ is predator attack rate and $e$ is the assimilation efficiency. The solution to this ODE is denoted as $N(t)$ and $P(t)$ and the observed values as $N_{obs}(t)$ and $P_{obs}(t)$ We write likelihoods as follows $N_{obs}(t)\sim Normal(N(t),\sigma)$ and $P_{obs}(t)\sim Normal(P(t),\sigma)$. So we assume that the only source of error is observational error.

# Formatting data for model fitting

```{r}
# keep only one replicate
repl=2 # replicate number we are keeping
N=data$n[data$replicate==repl] # time series of for spp1
P=data$p[data$replicate==repl] #time series of spp2
t=data$time[data$replicate==repl] # times
n=length(data$time[data$replicate==repl]) #size of data set

#rstan reads data as a named list
data_rstan=list(n=n,t=t,N=N,P=P)
```

# Translate to rstan
```{r}
model_competition_str='
  //write function for ode
  functions{
    real[] odemodel(real t, real[] N, real[] p, real[] x_r, int[] x_i){
      real dNdt[2];
      //p[1]=r, p[2]=a,p[3]=e,p[4]=d
      dNdt[1]=p[1]*N[1]-p[2]*N[1]*N[2];
      dNdt[2]=p[2]*p[3]*N[1]*N[2]-p[4]*N[2];
      return dNdt;
    }
  }
  //data 
  data
  {
    //make sure the names are the same as the list in R
    int n;
    real t[n];
    real N[n];
    real P[n];
  }
  //parameters that have to be estimated go here
  parameters
  {
    real<lower=0> r;
    real<lower=0> a;
    real<lower=0> e;
    real<lower=0> d;
    real<lower=0> sigma;
    real<lower=0> Ninit;
    real<lower=0> Pinit;
  }
  //model
  model
  {
    real p[4]; //store parameters to pass to ode 
    real N_sim[n-1,2];  //store simulated values
    
    //write priors
    r~lognormal(-0.5,1);
    a~lognormal(-0.5,1);
    e~lognormal(-0.5,1);
    d~lognormal(-0.5,1);
    sigma~gamma(2,0.1);
    Ninit~normal(N[1],1);
    Pinit~normal(P[1],1);
    //parameters for integrator
    p[1]=r;
    p[2]=a;
    p[3]=e;
    p[4]=d;
    
    //integrate ode
    N_sim=integrate_ode_rk45(odemodel,{Ninit,Pinit},t[1],t[2:n],p,rep_array(0.0,0),rep_array(0,0));
    
    //likelihood for initial value
    N[1]~normal(Ninit,sigma);
    P[1]~normal(Pinit,sigma);
    for(i in 2:n)
    {
      N[i]~normal(N_sim[i-1,1],sigma);
      P[i]~normal(N_sim[i-1,2],sigma);
    }
  }
  generated quantities{
  }
  '
```
Compile the model
```{r}
model=stan_model(model_code=model_competition_str,auto_write = TRUE)
```

# Fit model to data using MCMC

```{r}
#stan options
chains=3
#rstan_options(auto_write=TRUE)
options(mc.cores=chains)
iter=6000
warmup=4000

#initial value for sampling
init=rep(list(list(r=0.1,a=0.1,e=0.1,d=0.1,sigma=2,Ninit=data_rstan$N[1],N2init=data_rstan$P[1])),chains)
fit=sampling(model,data=data_rstan,iter=iter,warmup=warmup,chains=chains,init=init)
```


# Model diagnostics
```{r}
print(fit,digits=3)
```
```{r}
params=c("r","a","e","d")
samples=As.mcmc.list(fit)
plot(samples[,params])
```
```{r}
pairs(fit,pars=params)
```

#Posterior predictions

We now solve the ODE for 1000 samples of parameter estimates.
```{r}
ode.model=function(t,N,p)
{
  r=p$r
  a=p$a
  e=p$e
  d=p$d
  dN=r*N[1]-a*N[1]*N[2]
  dP=a*e*N[1]*N[2]-d*N[2]
  return(list(c(dN,dP)))
}

posteriors=as.matrix(fit) #posterior predictions

n_post=1000 # number of samples drawn
times=seq(min(data_rstan$t),max(data_rstan$t),length.out=50)
predictions=data.frame()
for(k in 1:n_post)
{
  par=posteriors[sample(1:nrow(posteriors),1),]
  sim=ode(c(par['Ninit'],par['Pinit']),times,ode.model,list(r=par["r"],a=par["a"],e=par["e"],d=par["d"]))
  temp=data.frame(sample_no=k,time=sim[,1],N=sim[,2],P=sim[,3])
  predictions=rbind(predictions,temp)
}
```

```{r}
#plot raw data
plot(data_rstan$t,data_rstan$N,col="blue",pch=16,ylim=c(0,50))
points(data_rstan$t,data_rstan$P,col="red",pch=16,ylim=c(0,50))

#plot posterior predictions
for(k in 1:n_post)
{
  lines(predictions$time[predictions$sample_no==k],predictions$N[predictions$sample_no==k],col=rgb(0,0,1,0.1))
  lines(predictions$time[predictions$sample_no==k],predictions$P[predictions$sample_no==k],col=rgb(1,0,0,0.1))
}
```