title: “R0 Sensitivity Analysis-c” output: html_notebook: default pdf_document: default —
rm(list=ls())
require(deSolve)## Loading required package: deSolverequire(lhs)## Loading required package: lhs## Warning: package 'lhs' was built under R version 4.0.2SARSCOV2Model <- 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 generateg.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/6params.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[,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)
cseq<-seq(.01,5.01, by=.25)
levels<-length(cseq)j <- 1
data <- data.frame(matrix(rep(NA,levels*h2*19),nrow=levels*h2))
for(i in 1:h2){
for (c in cseq){
data[j,1:16] <- params <- as.list(c(params.set[i,], c=c, 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$c, data$R0, pch=19, cex=0.3, col='blue', xlab='c', ylab='R0')
plot(data$c, data$outbreak.size/100000, pch=19, cex=0.3, col='blue', xlab='c', ylab='Outbreak Size/100,000') #divide epidemic final size by 100000 to get proportion
boxplot(data$outbreak.size/100000~data$c, border='blue', log='y', pch='.', xlab='c', ylab='Outbreak Size/100,000')
boxplot(data$R0~data$c, border='blue', log='y', pch='.', xlab='c', ylab='R0')
library(sensitivity)## Warning: package 'sensitivity' was built under R version 4.0.2bonferroni.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.145649085 2.822405e-04 0.023284401 -0.22121999 -0.07822053
## q 0.821431693 -1.329810e-04 0.007954702 0.79778226 0.84685369
## p 0.032347380 -4.281010e-04 0.020274926 -0.02634433 0.09866706
## delta -0.465041436 -4.989628e-04 0.019150508 -0.52355488 -0.40531261
## b -0.894894914 -6.515652e-05 0.004977381 -0.91043001 -0.88115496
## k -0.547475194 -5.275587e-04 0.015071983 -0.60265964 -0.50489270
## beta 0.424105730 -2.463223e-04 0.021540785 0.35753945 0.49397581
## omega 0.439177239 1.923740e-05 0.018253208 0.38194332 0.49501807
## aH -0.018339282 4.283471e-05 0.022247523 -0.08865049 0.04448960
## aQ 0.150122876 -2.180502e-04 0.020888493 0.09170843 0.22682984
## h -0.117798081 -4.528155e-04 0.023096407 -0.19125029 -0.05278914
## gammah 0.047804934 -3.137001e-04 0.020728295 -0.01635219 0.11819243
## gammaQ 0.028491411 -4.049251e-04 0.022082051 -0.03431225 0.09099086
## lambda -0.004873354 -1.019375e-03 0.025798081 -0.09044874 0.07426037
## c 0.735702243 8.382347e-04 0.010243044 0.70539026 0.76302176
## S0 -0.024111531 2.352208e-02 0.021338266 -0.14707972 0.02998535par(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)