SIR model in phydynR

sir_phydynR.Rmd

---
title: "SIR in phydynR"
author: "Simon Frost"
date: '`r Sys.Date()`'
output: github_document
---

```{r setup, include=FALSE}
library(simecol)
library(phydynR)
library(reshape2)
library(ggplot2)
library(ggtree)
library(phydynR)
```


## simecol ODE model

```{r}
sir <- new("odeModel",
  main = function(time, init, parms, ...){
    with(as.list(c(init,parms)),{
      dS <- mu*(S+I+R)-beta*S*I-mu*S
      dI <- beta*S*I-(mu+gamma)*I
      dR <- gamma*I-mu*R
      list(c(dS,dI,dR))
    })},
  parms = c(beta=0.001,gamma=0.1,mu=0.1),
  times = c(from=0,to=20,by=0.1),
  init = c(S=999,I=1,R=0),
  solver = "lsoda"
)
```

```{r}
sir <- sim(sir)
sir.out <- out(sir)
sir.out.long <- melt(sir.out,id="time")
```

```{r}
ggplot(sir.out.long,aes(x=time,y=value,group=variable))+
  geom_line(aes(color=variable))
```

## phydynR

All vectors and matrices have to be named.

```{r}
INFECTEDNAMES <- c("I")
```

```{r}
births <-  c('parms$beta*S*I')
names(births) <- INFECTEDNAMES
```

If we had more than one type of infected individual, `births` would be a matrix.

```{r,eval=FALSE}
births <-  matrix('0',nrow=k,ncol=k)
# fill in entries of birth matrix
rownames(births)=colnames(births)<-INFECTEDNAMES
```

In addition, if there is more than one type of infected individual, we would have to define a 'migration' matrix (G in Erik's nomenclature).

With deaths, we just refer to deaths of infected individuals.

```{r}
deaths <- c('(parms$mu+parms$gamma)*I')
names(deaths) <- INFECTEDNAMES
```

To fully specify the system, we have to give the dynamics of other compartments.

```{r}
nonDemeDynamics <- c('parms$mu*(S+I+R)-parms$beta*S*I-parms$mu*S',
                     'parms$gamma*I-parms$mu*R')
names(nonDemeDynamics) <- c("S","R")
```

We first build a demographic process.

```{r}
sir.dm <- build.demographic.process(births=births
  , deaths = deaths
  , nonDemeDynamics = nonDemeDynamics
  , parameterNames=c('beta', 'mu', 'gamma')
  , rcpp=FALSE
  , sde=FALSE)
```

This generates a function, which can be called as follows to generate trajectories that should match the vanilla ODE model. The fifth entry in the output is a data frame of variables.

```{r}
sir.dm.result <- sir.dm(
  theta = list(beta = 0.001, mu = 0.1, gamma = 0.1),
  x0 = c(S=999, I=1, R=0),
  t0=0, t1=20, res = 1000, integrationMethod = 'adams')
sir.dm.out.long <- melt(as.data.frame(sir.dm.result[[5]])[,c("time","S","I","R")],id="time")
```

```{r}
p <- ggplot(sir.dm.out.long,aes(x=time,y=value,group=variable))+
  geom_line(aes(color=variable))
p
```

What's special about this is that you can also simulate trees. You have to specify when you want sampling.

```{r}
set.seed(1234)
sir.tre <- sim.co.tree(
  list(beta = 0.001, mu = 0.1, gamma = 0.1),
  sir.dm,
  x0 = c(S=999, I=1, R=0),
  t0 = 0,
  sampleTimes = seq(10, 15, length.out=50),
  res = 1000
) 
```


```{r}
sir.tre2 <- sir.tre
class(sir.tre2) <- "phylo"
p2 <- ggtree(sir.tre2)+theme_tree2()
p2
```

```{r}
p3 <- plot_grid(p,p2,rel_widths=c(1.5,1))
ggsave("tree.png",p3,width=8,height=4,units="in")
p3
```