Solution_competition_jhelam.Rmd

---
title: "Solution: Competitive 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/competition.csv")
#view data set
head(data)
```
```{r}
#plot time series
plot(data$time,data$n1,col="red",xlim=c(0,200),ylim=c(0,200),ylab="N1 or N2",xlab="Time",pch=16)
points(data$time,data$n2,col="green",pch=16)
```


# Formulate the model

We fit the following system of ODE to the data.
$$\frac{dN_1}{dt}=r_1N_1(1-\frac{N_1+\alpha_{21}N_2}{K_1})$$
$$\frac{dN_2}{dt}=r_2N_2(1-\frac{N_2+\alpha_{12}N_1}{K_1})$$
The solution to this ODE is denoted as $N_1(t)$ and $N_2(t)$ and the observed values as $N_{1,obs}(t)$ and $N_{2,obs}(t)$ We write likelihoods as follows $N_{1,obs}(t)\sim Normal(N_1(t),\sigma)$ and $N_{2,obs}(t)\sim Normal(N_2(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
N1=data$n1[data$replicate==repl] # time series of for spp1
N2=data$n2[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,N1=N1,N2=N2)
```

# 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]=r1, p[2]=K1,p[3]=alpha21,p[4]=r2, p[5]=K2,p[6]=alpha12
      dNdt[1]=p[1]*N[1]*(1-((N[1]+p[3]*N[2])/p[2]));
      dNdt[2]=p[4]*N[2]*(1-((N[2]+p[6]*N[1])/p[5]));
      return dNdt;
    }
  }
  //data 
  data
  {
    //make sure the names are the same as the list in R
    int n;
    real t[n];
    real N1[n];
    real N2[n];
  }
  //parameters that have to be estimated go here
  parameters
  {
    real<lower=0> r1;
    real<lower=0> K1;
    real<lower=0> alpha21;
    real<lower=0> r2;
    real<lower=0> K2;
    real<lower=0> alpha12;
    real<lower=0> sigma;
    real<lower=0> N1init;
    real<lower=0> N2init;
  }
  //model
  model
  {
    real p[6]; //store parameters to pass to ode 
    real N_sim[n-1,2];  //store simulated values
    
    //write priors
    r1~lognormal(-0.5,1);
    K1~lognormal(5,1);
    alpha21~lognormal(0.5,1);
    r2~lognormal(-0.5,1);
    K2~lognormal(5,1);
    alpha12~lognormal(0.5,1);
    sigma~gamma(2,0.1);
    N1init~normal(N1[1],1);
    N2init~normal(N2[1],1);
    //parameters for integrator
    p[1]=r1;
    p[2]=K1;
    p[3]=alpha21;
    p[4]=r2;
    p[5]=K2;
    p[6]=alpha12;
    
    //integrate ode
    N_sim=integrate_ode_rk45(odemodel,{N1init,N2init},t[1],t[2:n],p,rep_array(0.0,0),rep_array(0,0));
    
    //likelihood for initial value
    N1[1]~normal(N1init,sigma);
    N2[1]~normal(N2init,sigma);
    for(i in 2:n)
    {
      N1[i]~normal(N_sim[i-1,1],sigma);
      N2[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(r1=0.1,K1=300,alpha21=1,r2=0.2,K2=200,alpha12=1,sigma=1,N1init=data_rstan$N1[1],N2init=data_rstan$N2[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("r1","K1","alpha21","r2","K2","alpha12")
samples=As.mcmc.list(fit)
plot(samples[,params])
```
```{r}
pairs(fit,pars=params)
```

#Posterior predictions

We now solve the ODE for 1000 samples of parameters.
```{r}
ode.model=function(t,N,p)
{
  r1=p$r1
  K1=p$K1
  alpha21=p$alpha21
  r2=p$r2
  K2=p$K2
  alpha12=p$alpha12
  dN1=r1*N[1]*(1-(N[1]+alpha21*N[2])/K1)
  dN2=r2*N[2]*(1-(N[2]+alpha12*N[1])/K2)
  return(list(c(dN1,dN2)))
}

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=200)
predictions=data.frame()
for(k in 1:n_post)
{
  par=posteriors[sample(1:nrow(posteriors),1),]
  sim=ode(c(par['N1init'],par['N2init']),times,ode.model,list(r1=par["r1"],K1=par["K1"],alpha21=par["alpha21"],r2=par["r2"],K2=par["K2"],alpha12=par["alpha12"]))
  temp=data.frame(sample_no=k,time=sim[,1],N1=sim[,2],N2=sim[,3])
  predictions=rbind(predictions,temp)
}
```

```{r}
#plot raw data
plot(data_rstan$t,data_rstan$N1,col="red",pch=16,ylim=c(0,200))
points(data_rstan$t,data_rstan$N2,col="green",pch=16,ylim=c(0,200))

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