Remove all objects from workspace.

Load add-on packages - deSolve - contains lsoda function - differential equation solver.

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]        # exposure
  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)
         } 
        )
}

Define Parameters

transmission_rate = beta_value =  0.7   # transmission rate per day
infectious_period = 4.5                 # infectious period
latent_period = 1.9                     # latent period
gamma_value = 1 / infectious_period     # per day
delta_value = 1 / latent_period         # per day
vaccine_compliance = 0    # basecase 
vaccine_efficacy = 0.5    # vaccine efficacy

vac.total = vaccine_efficacy * vaccine_compliance

Disease dynamics parameters.

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

Initial values for sub-populations.

X = 9999        # susceptible hosts
Y = 1           # infectious hosts
Z = 0           # recovered hosts
W = 0           # exposure host

Compute total population.

N = X + Y + Z + W

Initial state values for the differential equations.

initial_values = c(S=(X*(1-vac.total))/N,E=W/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)

Basecase plot

plot(S~time,data=output,type="l",lwd=2,col=1,ylab="",ylim=c(0,1),main="Basecase")
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)

Base case attack rate and treatment costs

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

#Summary
print(paste("Basecase attack rate=",round(Attack.rate,3)))

## [1] "Basecase attack rate= 0.95"

print(paste("Basecase cost = $", round(basecase_cost,3)))

## [1] "Basecase cost = $ 4748990.569"

Description (Basecase)

Attack rate of base case is 0.95 and total health case cost of treatment is 4,748,990.569

Vaccination intervention

Difference 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=W/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("Basecase attack rate =",round(Attack.rate,3),"Basecase 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] "Basecase attack rate = 0.95 Basecase 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] "Basecase attack rate = 0.95 Basecase 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"

Description (vaccination compliance rate 40%)

If vaccine efficacy is 40%, the attack rate will be 0.764. The total health care cost (treatment+vaccine) is $3,853,018.431. Cost per case averted (ICER) is $ 482.76.
If vaccine efficacy is 60%, the attack rate will be 0.666. The total health care cost (treatment+vaccine) is $3,379,192.212. Cost per case averted (ICER) is $ 483.074.

Vaccination compliance 50% and 60%

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=W/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"

Description (Varies of vaccination compliance and vaccine efficacy)

If vaccine efficacy is 40% ;

vaccine compliance 50%

The attack rate is 0.716, The total health care cost (treatment+vaccine) is $3,618,964.699. cost per case averted (ICER) is $ 482.908.

Vaccine compliance 60%

The attack rate is 0.666, The total health care cost (treatment+vaccine) is $3,379,192.212. cost per case averted (ICER) is $ 483.074.

If vaccine efficacy is 60% ;

vaccine compliance 50%

The attack rate is 0.589, The total health care cost (treatment+vaccine) is $3,004,101.699. cost per case averted (ICER) is $ 483.38.

Vaccine compliance 60%

The attack rate is 0.505, The total health care cost (treatment+vaccine) is $2,594,860.887. cost per case averted (ICER) is $ 483.83.

Summary

ICER <- structure(list('vac efficacy 40%' = c(482.76, 482.908, 483.074),'vac efficacy 60%' =  c(483.074, 483.38, 483.83)), .Names = c("vaccine efficacy 40%", "vaccine efficary 60%"), row.names = c("compliance 40%", "compliance 50%", "compliance 60%"), class="data.frame")

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 efficary 60%"), row.names = c("compliance 40%", "compliance 50%", "compliance 60%"), class="data.frame")

Comparing cost($) per case averted (ICER)

##                vaccine efficacy 40% vaccine efficary 60%
## compliance 40%              482.760              483.074
## compliance 50%              482.908              483.380
## compliance 60%              483.074              483.830

Comparing attack rate

##                vaccine efficacy 40% vaccine efficary 60%
## compliance 40%                0.764                0.666
## compliance 50%                0.712                0.589
## compliance 60%                0.666                0.505

Description

In summary, when we implement influenza intervention with vaccination program, the attack rate decrease from 0.95 (base case) to 0.505-0.764 depening on vaccination efficacy and compliance. The higher vaccine efficacy and compliance result in the lower attack rate. The lowest attack rate is 0.505 which vaccination should be 60% efficacy and 60% compliance. Consideration cost, the lowest cost per case averted os $482.760 whem implement with vaccine efficacy 40% and compliance at 40%. When we implement intervention with vaccination efficacy 60% and compliance 60%, it seems to be highest cost per case averted ($483.830). However, cost per case averted does not show significant different ($482.760 to $483.830)

Lab9_1

Naree Ketusing

October 28, 2015