Remove all objects from workspace.

r remove (list = objects() ) Load add-on packages - deSolve - contains lsoda function - differential equation solver.

library (deSolve)
## 
## Attaching package: 'deSolve'
## 
## The following object is masked from 'package:graphics':
## 
##     matplot

Function to compute derivatives of the differential equations.

seir_model = function (current_timepoint, state_values, parameters)
{
  # create state variables (local variables)
  S = state_values [1]        # susceptibles
  E = state_values [2]        # exposed
  I = state_values [3]        # infectious
  R = state_values [4]        # recovered
  
  with ( 
    as.list (parameters),     # variable names within parameters can be used 
         {
           # compute derivatives
           dS = (-beta * S * I)
           dE = (beta * S * I) - (delta * E)
           dI = (delta * E) - (gamma * I)
           dR = (gamma * I)
           
           # combine results
           results = c (dS, dE, dI, dR)
           list (results)
         }
    )
}

Parameters

transmission_rate = beta_value = 0.7   # transmission rate
infectious_period = 4.5                 # infectious period
latent_period = 1.9                    # latent period
vaccine_compliance = 0                 # base case 
vaccine_efficacy = 0.5               # vaccine efficacy

vac.total = vaccine_efficacy * vaccine_compliance

Compute values of gamma (recovery rate) and delta (exposion rate)

beta_value = transmission_rate
gamma_value = 1 / infectious_period
delta_value = 1 / latent_period

Compute Ro - Reproductive number.

Ro = beta_value / gamma_value

Disease dynamics parameters.

parameter_list = c (beta = beta_value, gamma = gamma_value, delta = delta_value)

Initial values for sub-populations.

V = 0        # exposed hosts
X = 9999        # susceptible hosts
Y = 1           # infectious hosts
Z = 0           # recovered hosts

Compute total population.

N = V + X + Y + Z

Initial state values for the differential equations.

initial_values = c (S=(X*(1-vac.total))/N,E=V/N,I=Y/N,R=(Z+X*vac.total)/N)

Output timepoints.

timepoints = seq (0, 100, by=1)

Simulate the SEIR epidemic.

output = lsoda (initial_values, timepoints, seir_model, parameter_list)

Plot SEIR Epidemic Base Case

plot(S~time,data=output,type="l",lwd=2,col=1,ylab="",ylim=c(0,1),main="SEIR Epidemic: Base Case")
lines(E~time,data=output,type="l",lwd=2,col=2,ylab="")
lines(I~time,data=output,type="l",lwd=2,col=3,ylab="")
lines(R~time,data=output,type="l",lwd=2,col=4,ylab="")
legend("right",legend=c("S","E","I","R"),col=c(1,2,3,4),lwd=2,cex=0.7)

Attack rate and treatment costs of base case

treatment_cost = 500  # per infected person
Attack.rate=sum(output[,4])/(1/gamma_value)
basecase_cost = treatment_cost*Attack.rate*N

#Summary
print(paste("Base case attack rate=",round(Attack.rate,3)))
## [1] "Base case attack rate= 0.95"
print(paste("Basecase cost = $", round(basecase_cost,3)))
## [1] "Basecase cost = $ 4748990.569"

The attack rate of the base case is 0.95 and the total treatment cost is $4,748.990.57.

Vaccination Intervention: Vaccination Efficacy

vaccine_compliance = 0.4
vac.eff = seq(0.4,0.6,0.2) # Define different vaccine efficacy
cost_per_vaccine = 20 # avg. cost of vaccine per person

par(mfrow=c(2,2))
for (i in vac.eff) {
  vaccine_efficacy = i
  vac.total = vaccine_efficacy * vaccine_compliance
  initial_values = c(S=(X*(1-vac.total))/N,E=V/N,I=Y/N,R=(Z+X*vac.total)/N)
  output = lsoda (initial_values, timepoints, seir_model, parameter_list)
  
  # plot intervention
  plot(S~time,data=output,type="l",lwd=2,col=1,ylab="",ylim=c(0,1),main=paste("vaccine efficacy=",i))
  lines(E~time,data=output,type="l",lwd=2,col=2,ylab="")
  lines(I~time,data=output,type="l",lwd=2,col=3,ylab="")
  lines(R~time,data=output,type="l",lwd=2,col=4,ylab="")
  legend("right",legend=c("S","E","I","R"),col=c(1,2,3,4),lwd=2,cex=0.7)

   # Cost
   vaccine_cost = cost_per_vaccine * vac.total * X
   new_Attack.rate=sum(output[,4])/(1/gamma_value)
 
   new_treatment_cost = treatment_cost*new_Attack.rate*N 
   total_cost = vaccine_cost + new_treatment_cost
  
   # Incremental cost effectiveness ratio
   ICER = (total_cost - basecase_cost) / {(new_Attack.rate*N) - (Attack.rate*N)}
    
   #Summary
   print(paste("Base case attack rate =",round(Attack.rate,3),"Base case cost =",round(basecase_cost,3)))
   print(paste("Vaccine compliance =",vaccine_compliance, "Vaccine efficacy =",i,", Attack.rate=",round(new_Attack.rate,3), "Total Healthcare Cost= $", round(total_cost,3),"Incremental cost effectiveness ratio = ", round(ICER,3)))
}
## [1] "Base case attack rate = 0.95 Base case cost = 4748990.569"
## [1] "Vaccine compliance = 0.4 Vaccine efficacy = 0.4 , Attack.rate= 0.764 Total Healthcare Cost= $ 3853018.431 Incremental cost effectiveness ratio =  482.76"
## [1] "Base case attack rate = 0.95 Base case cost = 4748990.569"
## [1] "Vaccine compliance = 0.4 Vaccine efficacy = 0.6 , Attack.rate= 0.666 Total Healthcare Cost= $ 3379192.212 Incremental cost effectiveness ratio =  483.074"

When the vaccine efficacy is 40%, the attack rate is 0.764. Total health care cost is $3,853,018.43 and ICER is $482.76 per case averted. When the vaccine efficacy is 60%, the attack rate is 0.666. Total health care cost is $3,379,192.21 and ICER is $483.07 per case averted.

vac.com = seq(0.5,0.6,0.1) # Define different vaccine compliance
vac.eff = seq(0.4,0.6,0.2) # Define different vaccine efficacy
cost_per_vaccine = 20 # avg. cost of vaccine per person

par(mfrow=c(2,2))
for (i in vac.eff) 
for (ii in vac.com)  {
  vaccine_efficacy = i
  vaccine_compliance = ii 
  vac.total = vaccine_efficacy * vaccine_compliance
  initial_values = c(S=(X*(1-vac.total))/N,E=V/N,I=Y/N,R=(Z+X*vac.total)/N)
  output = lsoda (initial_values, timepoints, seir_model, parameter_list)
  
  # plot intervention
  plot(S~time,data=output,type="l",lwd=2,col=1,ylab="",ylim=c(0,1),main=paste("vaccine efficacy=",i, "compliance=",ii))
  lines(E~time,data=output,type="l",lwd=2,col=2,ylab="")
  lines(I~time,data=output,type="l",lwd=2,col=3,ylab="")
  lines(R~time,data=output,type="l",lwd=2,col=4,ylab="")
  legend("right",legend=c("S","E","I","R"),col=c(1,2,3,4),lwd=2,cex=0.7)

   # Cost
   vaccine_cost = cost_per_vaccine * vac.total * X
   new_Attack.rate=sum(output[,4])/(1/gamma_value)
 
   new_treatment_cost = treatment_cost*new_Attack.rate*N 
   total_cost = vaccine_cost + new_treatment_cost
  
   # Incremental cost effectiveness ratio
   ICER = (total_cost - basecase_cost) / {(new_Attack.rate*N) - (Attack.rate*N)}
   
   #Summary
   print(paste("Basecase attack rate =",round(Attack.rate,3),"Basecase cost =",round(basecase_cost,3)))
   print(paste("Vaccine compliance=",ii, "Vaccine efficacy=",i, "Attack.rate=",round(new_Attack.rate,3), "Total Healthcare Cost= $", round(total_cost,3),"Incremental cost effectiveness ratio = ", round(ICER,3)))
}
## [1] "Basecase attack rate = 0.95 Basecase cost = 4748990.569"
## [1] "Vaccine compliance= 0.5 Vaccine efficacy= 0.4 Attack.rate= 0.716 Total Healthcare Cost= $ 3618964.699 Incremental cost effectiveness ratio =  482.908"
## [1] "Basecase attack rate = 0.95 Basecase cost = 4748990.569"
## [1] "Vaccine compliance= 0.6 Vaccine efficacy= 0.4 Attack.rate= 0.666 Total Healthcare Cost= $ 3379192.212 Incremental cost effectiveness ratio =  483.074"
## [1] "Basecase attack rate = 0.95 Basecase cost = 4748990.569"
## [1] "Vaccine compliance= 0.5 Vaccine efficacy= 0.6 Attack.rate= 0.589 Total Healthcare Cost= $ 3004101.699 Incremental cost effectiveness ratio =  483.38"

## [1] "Basecase attack rate = 0.95 Basecase cost = 4748990.569"
## [1] "Vaccine compliance= 0.6 Vaccine efficacy= 0.6 Attack.rate= 0.505 Total Healthcare Cost= $ 2594860.887 Incremental cost effectiveness ratio =  483.83"

When the vaccine efficacy is 40% and compliance 50%, the attack rate is 0.716. Total health care cost is $3,618,964.70 and ICER is $482.91 per case averted. When the vaccine efficacy is 40% and compliance 60%, the attack rate is 0.666. Total health care cost is $3,379,192.21 and ICER is $483.07 per case averted.

When the vaccine efficacy is 60% and compliance 50%, the attack rate is 0.589. Total health care cost is $3,004,101.70 and ICER is $483.38 per case averted. When the vaccine efficacy is 60% and compliance 60%, the attack rate is 0.505. Total health care cost is $2,594,860.89 and ICER is $483.83 per case averted.

Summarize findings

#Summary
print(ICER <- structure (list('vac efficacy 40%' = c(482.76, 482.91, 483.07),'vac efficacy 60%' =  c(483.07, 483.38, 483.83)), .Names = c("vaccine efficacy 40%", "vaccine efficacy 60%"), row.names = c("compliance 40%", "compliance 50%", "compliance 60%"), class="data.frame"))
##                vaccine efficacy 40% vaccine efficacy 60%
## compliance 40%               482.76               483.07
## compliance 50%               482.91               483.38
## compliance 60%               483.07               483.83
print(Attack_rate <- structure (list('vac efficacy 40%' = c(0.764, 0.712, 0.666),'vac efficacy 60%' =  c(0.666, 0.589, 0.505)), .Names = c("vaccine efficacy 40%", "vaccine efficacy 60%"), row.names = c("compliance 40%", "compliance 50%", "compliance 60%"), class="data.frame"))
##                vaccine efficacy 40% vaccine efficacy 60%
## compliance 40%                0.764                0.666
## compliance 50%                0.712                0.589
## compliance 60%                0.666                0.505

Upon the implementation of the vaccination program, the attack rates decrease from 0.95 in the base case to a range of 0.505-0.764 with differing efficacies and compliances of the vaccinations. Higher efficacy and compliance resulthed in a lower attack rate. On the other hand, higher efficacy and compliance resulted in a higher cost per case averted (ICER). For the vaccine intervention with both 60% efficacy and compliance and an attack rate of 0.505, the ICER was $483.83. This intervention had the highest cost per case averted as opposed to the vaccine intervention with both 40% efficacy and compliance and an attack rate of 0.764, the ICER was $482.76, which was the lowest cost per case averted.