R/shinyorv.r

In epimdr2: Functions and Data for "Epidemics: Models and Data in R (2nd Edition)"

#' Launch a shiny-app to study outbreak-response vaccination campaigns
#' @details
#' Launch app for details
#' @examples
#' if(interactive()){orv.app}
#' @export
orv.app=shinyApp(ui=navbarPage("ORV",
  tabPanel("Intervention day",
    sidebarLayout(
    sidebarPanel(
sliderInput("vaccine_target1", "target:", 0.7,
              min = 0, max = 1),
sliderInput("vaccine_efficacy1", "efficacy:", 0.9,
              min = 0, max = 1),
numericInput("intervention_length1", "duration:", 14,
              min = 1, max = 100),
numericInput("N1", "N:", 1E5,
              min = 1E2, max = 1E6),
numericInput("mtime1", "endtime:", 120,
              min = 10, max = 3*365),
    sliderInput("R1", "R", 
              min = 1, max = 20, value=4),
    sliderInput("IP1", "Infectious period (days)", 5,
              min = 1, max = 100),
    sliderInput("LP1", "Latent period (days):", 8,
              min = 1, max = 100)),
    mainPanel(plotOutput("plot1"))
  )),
   tabPanel("R sensitivty",
    sidebarLayout(
    sidebarPanel(
    sliderInput("R2", "R", 
              min = 1, max = 20, value=4),
numericInput("pm2", "+/-:", 0.5,
              min = 1, max = 10),
sliderInput("vaccine_target2", "target:", 0.7,
              min = 0, max = 1),
sliderInput("vaccine_efficacy2", "efficacy:", 0.9,
              min = 0, max = 1),
numericInput("intervention_length2", "duration:", 14,
              min = 1, max = 100),
numericInput("N2", "N:", 1E5,
              min = 1E2, max = 1E6),
numericInput("mtime2", "endtime:", 120,
              min = 10, max = 3*365),
    sliderInput("IP2", "Infectious period (days)", 5,
              min = 1, max = 100),
    sliderInput("LP2", "Latent period (days):", 8,
              min = 1, max = 100)),
    mainPanel(plotOutput("plot2"))
  )),
   tabPanel("Duration sensitivity",
    sidebarLayout(
    sidebarPanel(
numericInput("intervention_length3", "duration:", 14,
              min = 1, max = 100),
numericInput("pm3", "+/-:", 7,
              min = 1, max = 21),
sliderInput("vaccine_target3", "target:", 0.7,
              min = 0, max = 1),
sliderInput("vaccine_efficacy3", "efficacy:", 0.9,
              min = 0, max = 1),
    sliderInput("R3", "R", 
              min = 1, max = 20, value=4),
numericInput("N3", "N:", 1E5,
              min = 1E2, max = 1E6),
numericInput("mtime3", "endtime:", 120,
              min = 10, max = 3*365),
    sliderInput("IP3", "Infectious period (days)", 5,
              min = 1, max = 100),
    sliderInput("LP3", "Latent period (days):", 8,
              min = 1, max = 100)),
    mainPanel(plotOutput("plot3"))
  )),
   tabPanel("Cover sensitivty",
    sidebarLayout(
    sidebarPanel(
sliderInput("vaccine_target4", "target:", 0.7,
              min = 0, max = 1),
numericInput("pm4", "+/-:", 0.1,
              min = 0, max = .9),
sliderInput("vaccine_efficacy4", "efficacy:", 0.9,
              min = 0, max = 1),
numericInput("intervention_length4", "duration:", 14,
              min = 1, max = 100),
    sliderInput("R4", "R", 
              min = 1, max = 20, value=4),
numericInput("N4", "N:", 1E5,
              min = 1E2, max = 1E6),
numericInput("mtime4", "endtime:", 120,
              min = 10, max = 3*365),
    sliderInput("IP4", "Infectious period (days)", 5,
              min = 1, max = 100),
    sliderInput("LP4", "Latent period (days):", 8,
              min = 1, max = 100)),
    mainPanel(plotOutput("plot4"))
  )),
   tabPanel("Retrospective analysis",
    sidebarLayout(
    sidebarPanel(
numericInput("day5", "Start day:", 60,
              min = 10, max = 3*365),
sliderInput("vaccine_target5", "target:", 0.7,
              min = 0, max = 1),
sliderInput("vaccine_efficacy5", "efficacy:", 0.9,
              min = 0, max = 1),
numericInput("intervention_length5", "duration:", 14,
              min = 1, max = 100),
sliderInput("R5", "R", 
              min = 1, max = 20, value=4),
numericInput("N5", "N:", 1E5,
              min = 1E2, max = 1E6),
numericInput("mtime5", "endtime:", 120,
              min = 10, max = 3*365),
    sliderInput("IP5", "Infectious period (days)", 5,
              min = 1, max = 100),
    sliderInput("LP5", "Latent period (days):", 8,
              min = 1, max = 100)),
    mainPanel(plotOutput("plot5"))
  )),
  tabPanel("Summary")
),

server=function(input, output){
######################################################
#SEIR model 
######################################################
simod<-function(t,x,parms){
################
#Parameters
#B = transmission rate
#1/r = latent period
#1/g = infectious period
#q = vaccine efficacy
#
#P = target coverage
#Dt = length of vaccination campaign
#T = Day of campaign start
#USAGE:
# times<-1:100
# xstrt<-c(S=.999,E=0,I=.001,R=0,K=0)
# par<-c(B=.5, r=1/7, g = 1/7, q = .8, P = 0, Dt = 10, T = 80)
# out<-as.data.frame(lsoda(xstrt,times,simod,par))
# plot(out$time,out$I,type="l")
#
#
# par<-c(B=.5, r=1/7, g = 1/7, q = .8, P = .99, Dt = 10, T = 50)
# out<-as.data.frame(lsoda(xstrt,times,simod,par))
# lines(out$time,out$I,col="red")

    S<-x[1]
    E<-x[2]
    I<-x[3]
    R<-x[4]
    K<-x[5]
    #
    with(as.list(parms),{
      Q<- ifelse(t<T | t>T+Dt,0,(-log(1-P)/Dt))
      dS<- -B*S*I-q*Q*S
      dE<- B*S*I-r*E
      dI<- r*E - g*I
      dR<- g*I+q*Q*S
      dK<-r*E
      res<-c(dS,dE,dI,dR,dK)
      list(res)
    })
  }
######################################################
  
#####################################################
#% Intervention
#####################################################
p_red<-function(R,vaccine_efficacy,target_vaccination,intervention_length, mtime=120, LP=7, IP=7, N=10000, step=1){
  steps<-(0:mtime)[seq(1,mtime,by=step)]
  p_red<-rep(NA,length(steps))
  xstrt<-c(S=1-1/N,E=0,I=1/N,R=0,K=0)   #starting values
  beta<- R/IP         #transmission rate
  t<-1
  for(i in 1:length(steps)){
    par<-c(B=beta, r=1/LP, g = 1/IP, q = vaccine_efficacy,
                        P = target_vaccination, Dt = intervention_length, T = steps[i])
    out<-as.data.frame(lsoda(xstrt,steps,simod,par))
    p_red[t]<-out$K[dim(out)[1]]
    t<-t+1
  }
  par<-c(B=beta, r=1/LP, g = 1/IP, q = vaccine_efficacy,
                P = 0, Dt = 0, T = Inf)
  outv<-as.data.frame(lsoda(xstrt,steps,simod,par))
        #fs should really be the prediction with steps=Inf?
  fs<-max(out$K)
  res<-list(out=cbind(steps,p_red/max(p_red)),
                  R=R,
                  vaccine_efficacy=vaccine_efficacy,
                  target_vaccination=target_vaccination,
                  intervention_length=intervention_length,
                  mtime=mtime, LP=LP, IP=IP, N=N, step=step,
                  virgin=outv$I, vfs=fs)
  class(res)<-"p_red"
  return(res)
}
######################################################

#####################################################
#plotting pred objects
#####################################################
plot.p_red<-function(object){
      plot(object$out[,1],object$out[,2],type="l", xlab="First intervention day", ylab="% final epidemic", ylim=c(0,1))
      title(paste("target= ", round(100*object$target_vaccination,0), "% campaign = ", object$intervention_length,"d"))
par(new=TRUE)
      plot(object$out[,1], object$virgin, col='red', axes=FALSE, xlab="", ylab="", type="l")
      legend(x="topleft", legend=c("natural epidemic", "%final size"), col=c("red", "black"), lty=c(1,1))
}

#EX
#out<-p_red(R=4,vaccine_efficacy=.9,target_vaccination=.5,intervention_length=14, step=2)
#plot(out)

#####################################################
#% Sensitivity analysis on R
#####################################################
R_compare<-function(R=c(2,4,8),vaccine_efficacy=.9,target_vaccination, intervention_length,mtime=120, LP=7, IP=7, N=10000, step=7){
  out<-numeric(0)
  for(j in 1:length(R)){
    tmp<-p_red(R=R[j],vaccine_efficacy=vaccine_efficacy,
                target_vaccination,intervention_length, mtime=mtime, LP=LP, IP=IP, N=N, step=step)
    out<-cbind(out,tmp$out[,2])
  }
  res<-list(R=R, p_red=out, T=tmp$out[,1],
                 vaccine_efficacy=vaccine_efficacy,
                 target_vaccination=target_vaccination,
                 intervention_length=intervention_length,
                 mtime=mtime, LP=LP, IP=IP, N=N, step=step)
  class(res)<-"Rcomp"
  return(res)
}

#####################################################
#plotting pred objects
#####################################################
plot.Rcomp<-function(object){
        plot(NA,xlim=range(object$T),ylim=c(0,1), xlab="First intervention day", ylab="%, final epidemic")
        title(paste("% final size: target= ", round(100*object$target_vaccination,0), "% campaign = ", object$intervention_length,"d"))
  for(j in 1:length(object$R)){
    lines(object$T,object$p_red[,j],lty=j)
  }
      legend(x="right", legend=c(object$R), lty=1:length(object$R), title="R=")

}

#EX
#res<-R_compare(R=c(1.5, 2.5, 3.5), vaccine_efficacy=.9,target_vaccination=.5,intervention_length=7)
#plot(res)

#####################################################
#% Sensitivity analysis on length of intervention
#####################################################
Int_compare<-function(R,vaccine_efficacy,target_vaccination,intervention_length=c(7,10,14),mtime=120, LP=7, IP=7, N=10000, step=7){
  out<-numeric(0)
  for(j in 1:length(intervention_length)){
    tmp<-p_red(R=R,vaccine_efficacy=vaccine_efficacy,
                target_vaccination,intervention_length[j], mtime=mtime, LP=LP, IP=IP, N=N, step=step)
    out<-cbind(out,tmp$out[,2])
  }
  res<-list(R=R, p_red=out, T=tmp$out[,1],
                 vaccine_efficacy=vaccine_efficacy,
                 target_vaccination=target_vaccination,
                 intervention_length=intervention_length,
                 mtime=mtime, LP=LP, IP=IP, N=N, step=step)
  class(res)<-"Intcomp"
  return(res)
}

#####################################################
#plotting Intcomp objects
#####################################################
plot.Intcomp<-function(object){
        plot(NA,xlim=range(object$T),ylim=c(0,1), xlab="First intervention day", ylab="% final epidemic")
        title(paste("% final size: target= ", round(100*object$target_vaccination,0), "%, R = ", object$R))
  for(j in 1:length(object$intervention_length)){
    lines(object$T,object$p_red[,j],lty=j)
  }
      legend(x="right", legend=c(object$intervention_length), lty=1:length(object$intervention_length), title="Campaign D:")

}

#EX
#res<-Int_compare(intervention_length=c(7,10,14),R=3, vaccine_efficacy=.9, target_vaccination=.5)
#plot(res)

#####################################################
#% Sensitivity analysis on target Vaccination
#####################################################
Vacc_compare<-function(R,vaccine_efficacy,target_vaccination=c(.50,.70,.90),intervention_length=7,mtime=120, LP=7, IP=7, N=10000, step=7){
  out<-numeric(0)
  for(j in 1:length(target_vaccination)){
    tmp<-p_red(R=R,vaccine_efficacy=vaccine_efficacy,
                target_vaccination[j],intervention_length, mtime=mtime, LP=LP, IP=IP, N=N, step=step)
    out<-cbind(out,tmp$out[,2])
  }
  res<-list(R=R, p_red=out, T=tmp$out[,1],
                 vaccine_efficacy=vaccine_efficacy,
                 target_vaccination=target_vaccination,
                 intervention_length=intervention_length,
                 mtime=mtime, LP=LP, IP=IP, N=N, step=step)
  class(res)<-"Vacccomp"
  return(res)
}

#####################################################
#plotting Vaccomp objects
#####################################################
plot.Vacccomp<-function(object){
        plot(NA,xlim=range(object$T),ylim=c(0,1), xlab="First intervention day", ylab="% final epidemic")
        title(paste("% final size: R = ", object$R))
  for(j in 1:length(object$target_vaccination)){
    lines(object$T,object$p_red[,j],lty=j)
  }
      legend(x="right", legend=c(round(100*object$target_vaccination,0)), lty=1:length(object$target_vaccination), title="Target %")
}

#EX
#res<-Vacc_compare(target_vaccination=c(.50,.70,.90), R=4, vaccine_efficacy=.9,intervention_length=7)
#plot(res)

#####################################################
#Mortality analysis on R
#####################################################
M_red<-function(object, case_fatality){
oldpar <- par(no.readonly = TRUE) # code line i
on.exit(par(oldpar)) # code line i + 1
        par(mfrow=c(2,1))
        plot(object)
  out<-(1-object$p_red)*object$N*case_fatality
  plot(NA,xlim=c(1,max(object$T)),ylim=c(0,max(out)), xlab="First intervention day", ylab="Extra survivors")
  for(j in 1:dim(out)[2]){
    lines(object$T, out[,j], lty=j)
  }
  title(paste("Reduced burden of mortality; N = ", object$N))
oldpar <- par(no.readonly = TRUE) # code line i
on.exit(par(oldpar)) # code line i + 1
        par(mfrow=c(1,1))
}

#EX
#res<-Vacc_compare(target_vaccination=c(.50,.70,.90), R=4, vaccine_efficacy=.9,intervention_length=7)
#M_red(res, 0.5)


#####################################################
#Retrospective
#####################################################
retro<-function(R,day, vaccine_efficacy,target_vaccination,intervention_length, mtime=120, LP=7, IP=7, N=10000){
  steps<-1:mtime
        out<-matrix(NA,nrow=mtime, ncol=3)
  xstrt<-c(S=1-1/N,E=0,I=1/N,R=0,K=0)   #starting values
  beta<- R/IP         #transmission rate
  t<-1
  par<-c(B=beta, r=1/LP, g = 1/IP, q = vaccine_efficacy,
                P = 0, Dt = 0, T = Inf)
  outv<-as.data.frame(lsoda(xstrt,steps,simod,par))
        #fsv and fsi should really be with steps=Inf?
  fsv<-max(outv$K)

  par<-c(B=beta, r=1/LP, g = 1/IP, q = vaccine_efficacy,
             P = target_vaccination, Dt = intervention_length, T = day)

        outi<-as.data.frame(lsoda(xstrt,steps,simod,par))

  fsi<-max(outi$K)
        out[,1]<-steps
        out[,2]<-outv$I
        out[,3]<-outi$I
        res<-list(out=out,
                  redn=fsi/fsv,
                  R=R,
                  vaccine_efficacy=vaccine_efficacy,
                  target_vaccination=target_vaccination,
                  intervention_length=intervention_length,
                  mtime=mtime, LP=LP, IP=IP, N=N, day=day)
  class(res)<-"retro"
  return(res)
}
######################################################

#####################################################
#plotting pred objects
#####################################################
plot.retro<-function(object){
      plot(object$out[,1],object$out[,2],type="l", ylim=c(0,max(object$out[,2])), xlab='day', ylab='prevalence')
      polygon(c(object$day, object$day, object$day+object$intervention_length,
               object$day+object$intervention_length), c(-0.1,1,1,-.1), col="gray")
      lines(object$out[,1],object$out[,2])
      lines(object$out[,1], object$out[,3], col='red')
      title(paste("final size: ", round(100*(object$redn),1), "% (R=",
               object$R,", target=", 100*object$target_vaccination,"%)", sep=""))
      legend(x="topright", legend=c("natural epidemic", "w intervention"),
               col=c("black", "red"), lty=c(1,1))
      text(x=object$day+object$intervention_length, y=0, pos=4,
                labels=paste(object$intervention_length,
                "d intervention from", object$day))
}


output$plot1 <- renderPlot({
out<-p_red(R=input$R1,input$vaccine_efficacy1,input$vaccine_target1,input$intervention_length1, input$mtime1, input$LP1, input$IP1, input$N1, step=1)
plot(out)}
)

output$plot2 <- renderPlot({
R2=c(input$R2-input$pm2, input$R2, input$R2+input$pm2)
R2[R2<0]=0
out2=R_compare(R=R2, input$vaccine_efficacy2,input$vaccine_target2,input$intervention_length2, 
  input$mtime2, input$LP2, input$IP2, input$N2, step=1)
plot(out2)}
)

output$plot3 <- renderPlot({
il=c(input$intervention_length3-input$pm3,
input$intervention_length3, input$intervention_length3+input$pm3)
il[il<0]=0
out3=Int_compare(R=input$R3, input$vaccine_efficacy3,input$vaccine_target3,il,  
  input$mtime3, input$LP3, input$IP3, input$N3, step=1)
plot(out3)}
)

output$plot4 <- renderPlot({
vt=c(input$vaccine_target4-input$pm4,
input$vaccine_target4, input$vaccine_target4+input$pm4)
vt[vt<0]=0
vt[vt>1]=1
out4=Vacc_compare(R=input$R4, input$vaccine_efficacy4,vt,input$intervention_length4, 
 input$mtime4, input$LP4, input$IP4, input$N4, step=1)
plot(out4)}
)

output$plot5 <- renderPlot({
out5<-retro(R=input$R5, day=input$day5, input$vaccine_efficacy5,input$vaccine_target5,
  input$intervention_length5, input$mtime5, input$LP5, input$IP5, input$N5)
plot(out5)}
)

}
)