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[,2]*(q.max-q.min)+q.min,
p=lhs[,3]*(p.max-p.min)+p.min,
delta=lhs[,4]*(delta.max-delta.min)+delta.min,
b = lhs[,5]*(b.max-b.min)+b.min,
k = lhs[,6]*(k.max-k.min)+k.min,
beta = lhs[,7]*(beta.max-beta.min)+beta.min,
c = lhs[,8]*(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)
oseq<-seq(1,21, by=1)
levels<-length(oseq)
j <- 1
data <- data.frame(matrix(rep(NA,levels*h2*19),nrow=levels*h2))
for(i in 1:h2){
for (omega in oseq){
data[j,1:16] <- params <- as.list(c(params.set[i,], omega=omega, 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$omega, data$R0, pch=19, cex=0.3, col='blue', xlab='omega', ylab='R0')
plot(data$omega, data$outbreak.size/100000, pch=19, cex=0.3, col='blue', xlab='omega', ylab='Outbreak Size/100,000') #divide epidemic final size by 100000 to get proportion
boxplot(data$outbreak.size/100000~data$omega, border='blue', log='y', pch='.', xlab='omega', ylab='Outbreak Size/100,000')
boxplot(data$R0~data$omega, border='blue', log='y', pch='.', xlab='omega', 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.
## g -0.145088914 2.556024e-04 0.023143377 -0.21678571 -0.07800839
## q 0.820701727 -1.276399e-04 0.007919103 0.79710038 0.84628276
## p 0.035871491 -3.988143e-04 0.020311180 -0.02237742 0.10318970
## delta -0.483855071 -4.651496e-04 0.018130557 -0.54044473 -0.42685005
## b -0.897744666 -7.088542e-05 0.004740938 -0.91259608 -0.88468657
## k -0.565166138 -5.383087e-04 0.014169750 -0.61477632 -0.52350900
## beta 0.427770525 -2.773383e-04 0.021384636 0.36076987 0.49803480
## c 0.490376696 1.902655e-05 0.017380071 0.43705942 0.54293928
## aH -0.010342115 -9.710555e-06 0.022277591 -0.08111216 0.05226755
## aQ 0.152641599 -1.180474e-04 0.020964726 0.09278583 0.22745008
## h -0.126465858 -4.208648e-04 0.022608751 -0.19878782 -0.06377433
## gammah 0.048592676 -3.240498e-04 0.020695699 -0.01308232 0.11972484
## gammaQ 0.023225642 -3.278860e-04 0.022143607 -0.03869742 0.08508096
## lambda -0.008206427 -1.057907e-03 0.025944508 -0.09414757 0.07115779
## omega 0.717105779 9.018675e-04 0.011056895 0.68370967 0.74700458
## S0 -0.023998895 2.339156e-02 0.021273758 -0.14674483 0.02367160
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)
