Clearing environment

rm(list = ls())

R0 function

SIRR0 <- function(b,K,g,p,q,beta,delta,c,S0,omega,aH,h,gammaH){
  R0 <- (b*K*(g*(p-1)^2 + p)*q*beta*delta + b*c*(-g*(p-1)^3 + p)*S0*omega + (g-1)*(p-1)^2*(aH*h + gammaH - h*gammaH)*(-K*q*beta*delta + c*(p-1)*S0*omega)) / (b*K*(b*g - (g-1)*aH*h + gammaH - h*gammaH)*delta)
}

Initial Values

R0 <- SIRR0(b=1/0.77, K=100000, g=0, p=0.2, q=0.2, beta=13, delta=1, c=0, S0=10000, omega=0, aH=1/2.7, h=0.00082, gammaH=1/2.7)

Plotting all parameters in R0

## b ##
vals.b <- seq(0,1,by=0.01)
R0vals.b <- SIRR0(b=vals.b,K=100000,g=0,p=0.2,q=0.2,beta=13,delta=1,c=0,S0=10000,omega=0,aH=1/2.7,h=0.00082,gammaH=1/2.7)
plot(R0vals.b~vals.b,pch=19,main="R0 of Changing b",xlab="b values",ylab="R0")

## K ##
vals.K <- seq(0,100000,by=1000)
R0vals.K <- SIRR0(b=1/0.77,K=vals.K,g=0,p=0.2,q=0.2,beta=13,delta=1,c=0,S0=10000,omega=0,aH=1/2.7,h=0.00082,gammaH=1/2.7)
plot(R0vals.K~vals.K,pch=19,main="R0 of Changing K",xlab="K values",ylab="R0")

## g ##
vals.g <- seq(0,1,by=0.01)
R0vals.g <- SIRR0(b=1/0.77,K=100000,g=vals.g,p=0.2,q=0.2,beta=13,delta=1,c=0,S0=10000,omega=0,aH=1/2.7,h=0.00082,gammaH=1/2.7)
plot(R0vals.g~vals.g,pch=19,main="R0 of Changing g",xlab="g values",ylab="R0")

## p ##
vals.p <- seq(0,1,by=0.01)
R0vals.p <- SIRR0(b=1/0.77,K=100000,g=0,p=vals.p,q=0.2,beta=13,delta=1,c=0,S0=10000,omega=0,aH=1/2.7,h=0.00082,gammaH=1/2.7)
plot(R0vals.p~vals.p,pch=19,main="R0 of Changing p",xlab="p values",ylab="R0")

## q ##
vals.q <- seq(0,1,by=0.01)
R0vals.q <- SIRR0(b=1/0.77,K=100000,g=0,p=0.2,q=vals.q,beta=13,delta=1,c=0,S0=10000,omega=0,aH=1/2.7,h=0.00082,gammaH=1/2.7)
plot(R0vals.q~vals.q,pch=19,main="R0 of Changing q",xlab="q values",ylab="R0")

## beta ##
vals.beta <- seq(0,20,by=0.2)
R0vals.beta <- SIRR0(b=1/0.77,K=100000,g=0,p=0.2,q=0.2,beta=vals.beta,delta=1,c=0,S0=10000,omega=0,aH=1/2.7,h=0.00082,gammaH=1/2.7)
plot(R0vals.beta~vals.beta,pch=19,main="R0 of Changing beta",xlab="beta values",ylab="R0")

## delta ##
vals.delta <- seq(0,1,by=0.01)
R0vals.delta <- SIRR0(b=1/0.77,K=100000,g=0,p=0.2,q=0.2,beta=13,delta=vals.delta,c=0,S0=10000,omega=0,aH=1/2.7,h=0.00082,gammaH=1/2.7)
plot(R0vals.delta~vals.delta,pch=19,main="R0 of Changing delta",xlab="delta values",ylab="R0")

## c ##
vals.c <- seq(0,1,by=0.01)
R0vals.c <- SIRR0(b=1/0.77,K=100000,g=0,p=0.2,q=0.2,beta=13,delta=1,c=vals.c,S0=10000,omega=0,aH=1/2.7,h=0.00082,gammaH=1/2.7)
plot(R0vals.c~vals.c,pch=19,main="R0 of Changing c",xlab="c values",ylab="R0")

## S0 ##
vals.S0 <- seq(0,10000,by=100)
R0vals.S0 <- SIRR0(b=1/0.77,K=100000,g=0,p=0.2,q=0.2,beta=13,delta=1,c=0,S0=vals.S0,omega=0,aH=1/2.7,h=0.00082,gammaH=1/2.7)
plot(R0vals.S0~vals.S0,pch=19,main="R0 of Changing S0",xlab="S0 values",ylab="R0")

## omega ##
vals.omega <- seq(0,1,by=0.01)
R0vals.omega <- SIRR0(b=1/0.77,K=100000,g=0,p=0.2,q=0.2,beta=13,delta=1,c=0,S0=10000,omega=vals.omega,aH=1/2.7,h=0.00082,gammaH=1/2.7)
plot(R0vals.omega~vals.omega,pch=19,main="R0 of Changing omega",xlab="omega values",ylab="R0")

## aH ##
vals.aH <- seq(0,1,by=0.01)
R0vals.aH <- SIRR0(b=1/0.77,K=100000,g=0,p=0.2,q=0.2,beta=13,delta=1,c=0,S0=10000,omega=0,aH=vals.aH,h=0.00082,gammaH=1/2.7)
plot(R0vals.aH~vals.aH,pch=19,main="R0 of Changing aH",xlab="aH values",ylab="R0")

## h ##
vals.h <- seq(0,0.001,by=0.00001)
R0vals.h <- SIRR0(b=1/0.77,K=100000,g=0,p=0.2,q=0.2,beta=13,delta=1,c=0,S0=10000,omega=0,aH=1/2.7,h=vals.h,gammaH=1/2.7)
plot(R0vals.h~vals.h,pch=19,main="R0 of Changing h",xlab="h values",ylab="R0")

## gammaH ##
vals.gammaH <- seq(0,1,by=0.01)
R0vals.gammaH <- SIRR0(b=1/0.77,K=100000,g=0,p=0.2,q=0.2,beta=13,delta=1,c=0,S0=10000,omega=0,aH=1/2.7,h=0.00082,gammaH=vals.gammaH)
plot(R0vals.gammaH~vals.gammaH,pch=19,main="R0 of Changing gammaH",xlab="gammaH values",ylab="R0")

## q & beta ##
R0vals.qbeta1 <- SIRR0(b=1/0.77,K=100000,g=0,p=0.2,q=0,beta=vals.beta,delta=1,c=0,S0=10000,omega=0,aH=1/2.7,h=0.00082,gammaH=1/2.7)
R0vals.qbeta2 <- SIRR0(b=1/0.77,K=100000,g=0,p=0.2,q=0.2,beta=vals.beta,delta=1,c=0,S0=10000,omega=0,aH=1/2.7,h=0.00082,gammaH=1/2.7)
R0vals.qbeta3 <- SIRR0(b=1/0.77,K=100000,g=0,p=0.2,q=0.4,beta=vals.beta,delta=1,c=0,S0=10000,omega=0,aH=1/2.7,h=0.00082,gammaH=1/2.7)
R0vals.qbeta4 <- SIRR0(b=1/0.77,K=100000,g=0,p=0.2,q=0.6,beta=vals.beta,delta=1,c=0,S0=10000,omega=0,aH=1/2.7,h=0.00082,gammaH=1/2.7)
R0vals.qbeta5 <- SIRR0(b=1/0.77,K=100000,g=0,p=0.2,q=0.8,beta=vals.beta,delta=1,c=0,S0=10000,omega=0,aH=1/2.7,h=0.00082,gammaH=1/2.7)
R0vals.qbeta6 <- SIRR0(b=1/0.77,K=100000,g=0,p=0.2,q=1,beta=vals.beta,delta=1,c=0,S0=10000,omega=0,aH=1/2.7,h=0.00082,gammaH=1/2.7)

op1 <- par(mar=c(6,6,2,2))
plot(R0vals.qbeta1~vals.q,type="b", col="blue", main="R0 of Changing q & beta", xlab="q & beta values", ylab = "R0", ylim=c(0,30))
points(R0vals.qbeta2~vals.q,type="b",col="red")
points(R0vals.qbeta3~vals.q,type="b",col="darkturquoise")
points(R0vals.qbeta4~vals.q,type="b",col="orange")
points(R0vals.qbeta5~vals.q,type="b",col="purple")
points(R0vals.qbeta6~vals.q,type="b",col="green")
legend(0, 30,legend = c("q=0","q=0.2","q=0.4","q=0.6","q=0.8","q=1"), col = c("blue", "red", "darkturquoise", "orange", "purple", "green"), lty=1, cex=0.8)

par(op1)

Results of plots

# columns are parameters; rows are highest, lowest, and ending

results.R0 <- data.frame(
  b = c(max(R0vals.b), min(R0vals.b), tail(R0vals.b,n=1)),
  K = c(max(R0vals.K), min(R0vals.K), tail(R0vals.K,n=1)),
  g = c(max(R0vals.g), min(R0vals.g), tail(R0vals.g,n=1)),
  p = c(max(R0vals.p), min(R0vals.g), tail(R0vals.p,n=1)),
  q = c(max(R0vals.q), min(R0vals.q), tail(R0vals.q,n=1)),
  beta = c(max(R0vals.beta), min(R0vals.beta), tail(R0vals.beta,n=1)),
  delta = c(max(R0vals.delta), min(R0vals.delta), tail(R0vals.delta,n=1)),
  c = c(max(R0vals.c), min(R0vals.c), tail(R0vals.c,n=1)),
  S0 = c(max(R0vals.S0), min(R0vals.S0), tail(R0vals.S0,n=1)),
  omega = c(max(R0vals.omega), min(R0vals.omega), tail(R0vals.omega,n=1)),
  aH = c(max(R0vals.aH), min(R0vals.aH), tail(R0vals.aH,n=1)),
  h = c(max(R0vals.h), min(R0vals.h), tail(R0vals.h,n=1)),
  gammaH = c(max(R0vals.gammaH), min(R0vals.gammaH), tail(R0vals.gammaH,n=1))
)

`.rowNamesDF<-`(results.R0,make.names=FALSE,c('Highest','Lowest','Ending'))