rm(list=ls())
require(deSolve)
## Loading required package: deSolve
require(lhs)
## Loading required package: lhs
## Warning: package 'lhs' was built under R version 4.0.2
SARSCOV2Model <- function (t, y, params) {
S.h<-y[1] #create local variable S, first element of y
E.h<-y[2]
I.h<-y[3]
S.l<-y[4]
E.l<-y[5]
I.l<-y[6]
Q<-y[7]
R<-y[8]
V<-y[9]
with(
as.list(params, y),
{
dS.h<--q*beta*(I.h+I.l)*S.h/(S.h+E.h+I.h+S.l+E.l+I.l+Q+R)-c*(1-exp((-1/k)*V))*S.h
dE.h<-q*beta*(I.h+I.l)*S.h/(S.h+E.h+I.h+S.l+E.l+I.l+Q+R)+c*(1-exp((-1/k)*V))*S.h-lambda*E.h
dI.h<-lambda*E.h-b*g*I.h-aH*h*(1-g)*I.h-gammah*(1-h)*(1-g)*I.h
dS.l<--(1-p)*q*beta*(I.h+I.l)*S.l/(S.h+E.h+I.h+S.l+E.l+I.l+Q+R)-(1-p)*c*(1-exp((-1/k)*V))*S.l
dE.l<-(1-p)*q*beta*(I.h+I.l)/(S.h+E.h+I.h+S.l+E.l+I.l+Q+R)+(1-p)*c*(1-exp((-1/k)*V))*S.l-lambda*E.l
dI.l<-lambda*E.l-b*I.l
dQ<-b*I.l+g*b*I.h-aQ*h*Q-gammaQ*(1-h)*Q
dR<-gammah*(1-h)*(1-g)*I.h+gammaQ*(1-h)*Q
dV<-omega*I.h+(1-p)*omega*I.l-delta*V
dy<-c(dS.h,dE.h,dI.h,dS.l,dE.l,dI.l,dQ,dR,dV) #combine results into one vector dy
list(dy)
}
)
}
plot.SARSCOV2Model<-function(out)
{
par(mfrow=c(3,2))
plot(out$time, out$S.h, xlab='Time', ylab='Susceptible(High Risk)',type='l',col='steelblue4', lwd=3)
plot(out$time, out$E.h, xlab='Time', ylab='Exposed(High Risk)', type='l', col='steelblue4', lwd=3)
plot(out$time, out$I.h, xlab='Time', ylab='Infectious(High Risk)', type='l', col='steelblue4', lwd=3)
plot(out$time, out$S.l, xlab='Time', ylab='Susceptible(Low Risk)', type='l', col='steelblue4', lwd=3)
plot(out$time, out$E.l, xlab='Time', ylab='Exposed(Low Risk)', type='l', col='steelblue4', lwd=3)
plot(out$time, out$I.l, xlab='Time', ylab='Infectious(Low Risk)', type='l', col='steelblue4', lwd=3)
plot(out$time, out$Q, xlab='Time', ylab='Quarantined', type='l', col='steelblue4', lwd=3)
plot(out$time, out$R, xlab='Time', ylab='Recovered', type='l', col='steelblue4', lwd=3)
plot(out$time, out$R, xlab='Time', ylab='V', type='l', col='steelblue4', lwd=3)
}
get.size <- function(out) {
100000-(tail(out$S.h,1)+tail(out$S.l,1)) #subtract number of susceptibles in high and low risk groups from total population size 100000
}
get.peak <- function (out)
{
max(out$I.h)
}
SIRR0<-function(b,k,g,p,beta,delta,c,S0=100000,omega,aH,gammaH,h,q){
with(
as.list(params),
{
R0<-(b*k*(g*(-1+ p)^2+p)*beta*q*delta+b*c*(-g*(-1+p)^3+p)*S0*omega+(-1+g)*(-1+p)^2*(aH*h+gammah-h*gammah)*(-k*beta*q*delta+c*(-1+p)*S0*omega))/(b*k*(b*g-(-1+g)*(aH*h+gammah-h*gammah))*delta)
}
)
}
times <- seq(0,1000,by=1)
params <- c(q=0.2,beta=13,p=.2,c=0,lambda=1/4.43,b=1/0.77,g=0,gammaQ=0.1,gammah=1/2.7,aQ=1/1.93,aH=1/2.7,h=0.00082,omega=0,delta=1,k=10000)
ystart <- c(S.h=(.2)*100000,E.h=0,I.h=1,S.l=(1-0.2)*100000,E.l=0,I.l=0,Q=0,R=0,V=0)
out <- as.data.frame(lsoda(ystart,times,SARSCOV2Model,params))
plot.SARSCOV2Model(out)
h1<-1000
set.seed(6242015) #this random seed can be changed
lhs<-maximinLHS(h1,14) #here 14 is number of parameters that we will generate
g.min<-0
g.max<-1
q.min<-0
q.max<-0.2
p.min<-0.1
p.max<-0.2
delta.min<-1/7
delta.max<-24
b.max <- 1/1.29 # from Peak et al. 2020
b.min <- 0
k.max <- 1000000
k.min <- 10000
beta.max <- 18
beta.min <- 10
c.max <- 5
c.min <- 0.01
#S0.max <- 100000 #no need for S0. S0 = 100000
#S0.min <- 1
omega.max <- 20
omega.min <- 1
aH.max <- 1/2
aH.min <- 1/10
aQ.max <- 1/1
aQ.min <- 1/9
h.max <- 200/100000 #default: 82/100000 are hospitalized
h.min <- 50/100000
gammah.max <- 1/2.5
gammah.min <- 1/10
gammaQ.max <- 1/9
gammaQ.min <- 1/15
lambda.max <- 1/4
lambda.min <- 1/6
params.set<-cbind(
#g=lhs[,1]*(g.max-g.min)+g.min,
q=lhs[,1]*(q.max-q.min)+q.min,
p=lhs[,2]*(p.max-p.min)+p.min,
delta=lhs[,3]*(delta.max-delta.min)+delta.min,
b = lhs[,4]*(b.max-b.min)+b.min,
k = lhs[,5]*(k.max-k.min)+k.min,
beta = lhs[,6]*(beta.max-beta.min)+beta.min,
c = lhs[,7]*(c.max-c.min)+c.min,
omega = lhs[,8]*(omega.max-omega.min)+omega.min,
#S0 = lhs[,9]*(S0.max-S0.min)+S0.min,
aH =lhs[,9]*(aH.max-aH.min)+aH.min,
aQ = lhs[,10]*(aQ.max-aQ.min)+aQ.min,
h = lhs[,11]*(h.max-h.min)+h.min,
gammah = lhs[,12]*(gammah.max-gammah.min)+gammah.min,
gammaQ = lhs[,13]*(gammaQ.max-gammaQ.min)+gammaQ.min,
lambda = lhs[,14]*(lambda.max-lambda.min)+lambda.min
)
h2<-100 #number of simulated parameter sets. We would like to increase this to at least 100 for our final plots
times <- seq(0,100,by=0.1)
gseq<-seq(0,1, by=.05)
levels<-length(gseq)
j <- 1
data <- data.frame(matrix(rep(NA,levels*h2*19),nrow=levels*h2))
for(i in 1:h2){
for (g in gseq){
data[j,1:16] <- params <- as.list(c(params.set[i,], g=g, S0=100000)) #16 parameters
out <- (as.data.frame(ode(ystart,times, SARSCOV2Model, params)))
data[j,17] <- get.size(out)
data[j,18] <- SIRR0(params)
data[j,19] <- get.peak(out)
j <- j+1
}
}
names(data) <- c(names(params),'outbreak.size', 'R0','get.peak')
#save(data, file='data.Rdata') remove the comment symbol when you want to save a large matrix (e.g. 100 simulations)
plot(data$g, data$R0, pch=19, cex=0.3, col='blue',xlab='g', ylab='R0')
plot(data$g, data$outbreak.size/100000, pch=19, cex=0.3, col='blue',xlab='g', ylab='Outbreak Size/100,000') #divide epidemic final size by 100000 to get proportion
boxplot(data$outbreak.size/100000~data$g, border='blue', log='y', pch='.', xlab='g', ylab='Outbreak Size')
boxplot(data$R0~data$g, border='blue', log='y', pch='.',xlab='g', ylab='R0')
library(sensitivity)
## Warning: package 'sensitivity' was built under R version 4.0.2
bonferroni.alpha <- 0.05/16
prcc <- pcc(data[,1:16], data[,18], nboot = 1000, rank=TRUE, conf=1-bonferroni.alpha) #how sensitive is R0 to each parameter?
save(prcc, file='prcc.Rdata')
library(sensitivity)
load('prcc.Rdata')
summary <- print(prcc)
##
## Call:
## pcc(X = data[, 1:16], y = data[, 18], rank = TRUE, nboot = 1000, conf = 1 - bonferroni.alpha)
##
## Partial Rank Correlation Coefficients (PRCC):
## original bias std. error min. c.i. max. c.i.
## q 0.914076252 2.106413e-04 0.002943550 0.90498581 0.92326115
## p -0.273916173 -6.597440e-05 0.019146192 -0.33516621 -0.22205080
## delta -0.586474441 -4.853018e-04 0.013823712 -0.62733894 -0.54493500
## b -0.905313229 -1.405890e-04 0.004020022 -0.91893796 -0.89392854
## k -0.645903963 -2.208778e-04 0.011567061 -0.68036608 -0.61333680
## beta 0.219664784 8.266802e-05 0.022918355 0.15370516 0.29120407
## c 0.569791284 2.502045e-04 0.014517861 0.53242990 0.61850087
## omega 0.470710241 -2.447438e-04 0.015995187 0.42614256 0.51594668
## aH 0.108919800 -1.504617e-03 0.022949911 0.03066017 0.17340690
## aQ 0.131598460 6.435486e-06 0.021077611 0.07096401 0.19797733
## h 0.054559609 -7.775627e-04 0.020526167 -0.01306810 0.13264090
## gammah -0.003205130 -1.557927e-03 0.021355621 -0.06547218 0.05823795
## gammaQ 0.098673801 -2.127998e-04 0.022926214 0.03812288 0.17081541
## lambda -0.045668786 -6.136419e-04 0.022491875 -0.11507821 0.02275751
## g -0.223115165 -4.922092e-04 0.021790051 -0.29138337 -0.16293282
## S0 -0.008848547 8.624981e-03 0.022065907 -0.09387836 0.04191612
par(mar=c(7,4,4,2)+0.1)
plot(summary$original, main='Partial rank correlation coefficients', ylim=c(-1,1),
xlab='', ylab='Coefficient',
axes=FALSE)
axis(2)
axis(1, at=seq(1:16), labels=row.names(summary), las=2)
mtext(text='Parameter', side=1, line=4.5)
box()
for(i in 1:16) lines(c(i,i),c(summary[i,4], summary[i,5]))
abline(h=0)
