library(gbd)
library(deSolve)
library(tidyverse)
## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.1.0     ✔ purrr   0.3.2
## ✔ tibble  3.0.3     ✔ dplyr   1.0.2
## ✔ tidyr   1.1.2     ✔ stringr 1.4.0
## ✔ readr   1.3.1     ✔ forcats 0.3.0
## ── Conflicts ───────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
theme_set(theme_bw())

Examples

d1<-run_sim(gamma = 0.2,days1 = 100,days2= 100,R01 =  3,R02 = 0.3,R012 = 0.1,R11 = 2.5,R12 = 1,R112 = 1)
d1$R01=3
d2<-run_sim(gamma = 0.2,days1 = 100,days2= 100,R01 =  1.6,R02 = 0.3,R012 = 0.1,R11 = 2.5,R12 = 1,R112 = 1)
d2$R01=1.6
d3<-run_sim(gamma = 0.2,days1 = 100,days2= 100,R01 =  1.4,R02 = 0.3,R012 = 0.1,R11 = 2.5,R12 = 1,R112 = 1)
d3$R01=1.4
d4<-run_sim(gamma = 0.2,days1 = 100,days2= 100,R01 =  1.1,R02 = 0.3,R012 = 0.1,R11 = 2.5,R12 = 1,R112 = 1)
d4$R01=1.1
d<-rbind(d1,d2,d3,d4)
d %>% pivot_longer(c(2:7)) %>% separate(name,into=c("Class", "group"), sep=1,remove=FALSE) %>% filter(Class !="S") -> dd
ggplot(dd,aes(x=time,y=value,col=name)) + geom_line() +facet_wrap(~R01) + geom_vline(aes(xintercept=100),col="grey")

n=1000

gamma = 0.2
days1 = 100
days2= 1
R01 =  rep(seq(1,5,length=10),each=10)
R02 = 0.3
R012 = rep(seq(0,1,length=10),times=10)
R11 = 2.5
R12 = 1
R112 = 1
d<-data.frame(gamma=gamma,days1,days2,R01=R01,R02=R02,R012=R012, R11=R11,R12=R12,R112=R112)
d %>% rowwise %>%  mutate(R1=final_R(R01=R01,R02=R02,R012=R012)$R1, R2=final_R(R01=R01,R02=R02,R012=R012)$R2) -> d
ggplot(d, aes(R012,R01)) +geom_tile(aes(fill = R2))

ggplot(d, aes(x=R01,y=R2)) +  geom_point()

gamma = 0.2
days1 = 100
days2= 1
R01 =  runif(1000,0.3,8)
R02 =  0.3
R012 = runif(1000,0.3,1)
R11 = 2.5
R12 = 1
R112 = 1
d<-data.frame(gamma=gamma,days1,days2,R01=R01,R02=R02,R012=R012, R11=R11,R12=R12,R112=R112)
d %>% rowwise %>%  mutate(R1=final_R(R01=R01,R02=R02,R012=R012)$R1, R2=final_R(R01=R01,R02=R02,R012=R012)$R2) -> d
d %>% filter(R2<0.4) -> dd

ggplot(dd,aes(y=R01,x=R012)) +geom_point()