In weinbergerlab/ExcessILI: A package and RShiny app to examine excess ILI from syndromic surveillance data

knitr::opts_chunk$set(
  collapse = TRUE,
    echo=F,
    warning=FALSE, 
    message=FALSE,
  comment = "#>"
)

library(ExcessILI)
library(cdcfluview)
library(reshape2)

Overview

The goal for this package is to facilitate the formatting of line list data from syndromic surveillance datasets into time series and then the analysis of these data to detect increases above the seasonal baseline. For US data, there is an option to automatically adjust the data for state-specific flu activity (using data from NREVSS and/or state-specific RSV activity (based on Google search volume). The user can either start with line list data or formatted time series data

In this example, we will analyze ILINet data with a simple seasonal baseline, adjusting for flu and year-to-year variations. The model is fit through end-of-February 2020 and then extrapolated forward based on the time of year and the amount of influenza activity. Influenza activity is captured by using the proportion of tests that were positive from the NREVSS system (log transformed in model)

## Download the data

ili.data <- ilinet(region = c("state"))

ili.data$state <- state.abb[match(ili.data$region, state.name)]
ili.data       <- ili.data[, c("state", "week_start", "ilitotal", "total_patients")]
ili.data       <- ili.data[!is.na(ili.data$total_patients),]
ili.data.spl   <- split(ili.data, ili.data$state)

min<-sapply(ili.data.spl, function(x)  min(x$total_patients))
min

state.select<-names(min)[which(min>0) ]
ili.data <- ili.data[ili.data$state %in% state.select,]

## Run the main analysis function, adjusting for flu using NREVSS data
ili.data$age <-1
excess_cases1 <-
  excessCases(ds = ili.data,
              datevar       = "week_start", 
              statevar      = "state",
              agevar        = 'age',
              denom.var     = "total_patients",
              adj.flu       = "auto",
              use.syndromes = c("ilitotal"),
              extrapolation.date = "2020-03-01",
              time.res='week')

excess_cases1 <-
  excessCases(ds = ili.data,
              datevar       = "week_start", 
              statevar      = "state",
              denom.var     = "total_patients",
              adj.flu       = "none",
              use.syndromes = c("ilitotal"),
              extrapolation.date = "2020-03-01",
              time.res='week')
## Run the main analysis function, adjusting for flu using NREVSS data
excess_cases.nb <-
  excessCases(ds = ili.data,
              datevar       = "week_start", 
              statevar      = "state",
              denom.var     = "total_patients",
              adj.flu       = "none",
              use.syndromes = c("ilitotal"),
              extrapolation.date = "2020-03-01",
              model.type = 'negbin',
              time.res='week')
pred1 <-
  excessExtract(ds = excess_cases1,
                syndrome = "ilitotal",
                extract.quantity = "pred")
pred2 <-
  excessExtract(ds = excess_cases.nb,
                syndrome = "ilitotal",
                extract.quantity = "pred")

pred.var1 <-
  excessExtract(ds = excess_cases1,
                syndrome = "ilitotal",
                extract.quantity = "upi")
pred.var2 <-
  excessExtract(ds = excess_cases.nb,
                syndrome = "ilitotal",
                extract.quantity = "upi")

plot(pred1, pred2)
abline(a=0, b=1)
plot(pred.var1, pred.var2)
abline(a=0, b=1)

## Plot the results in an interactive dashboard

# dashboardPlot(excess_cases1)

## Extract the quantities of interest

dates1 <-
  excess_cases1[[1]][[1]][[1]]$date

dates <-
  excess_cases1[[1]][[1]][[1]]$date

unexplained.cases <-
  excessExtract(ds = excess_cases1,
                syndrome = "ilitotal",
                extract.quantity = "unexplained.cases")

unexplained.log.rr <-
  excessExtract(ds = excess_cases1,
                syndrome = "ilitotal",
                extract.quantity = "resid1")

denom <-
  excessExtract(ds = excess_cases1,
                syndrome = "ilitotal",
                extract.quantity = "denom")


upi <-
  excessExtract(ds = excess_cases1,
                syndrome = "ilitotal",
                extract.quantity = "upi")

lpi <-
  excessExtract(ds = excess_cases1,
                syndrome = "ilitotal",
                extract.quantity = "lpi")

obs <-
  excessExtract(ds = excess_cases1,
                syndrome = "ilitotal",
                extract.quantity = "y")

pred <-
  excessExtract(ds = excess_cases1,
                syndrome = "ilitotal",
                extract.quantity = "pred")

result.object <-
  list('dates'=dates,
       'obs'=obs[,,1],
       'pred'=pred[,,1],
       'unexplained.cases'=unexplained.cases[,,1],
       'unexplained.log.rr'=unexplained.log.rr[,,1])
rr <-  excessExtract(ds = excess_cases1,
                syndrome = "ilitotal",
                extract.quantity = "resid1")
# saveRDS(result.object,'extracted.output.ilinet.rds')

#Extract ILI data
ili.log.rr <- 
  unexplained.log.rr[,,1]
ili.log.rr <- 
  as.data.frame(ili.log.rr)
ili.log.rr$date <- dates1
library(MMWRweek)
ili.mmwrweek <- MMWRweek(dates1)
ili.log.rr <- cbind(ili.log.rr, ili.mmwrweek)
write.csv(ili.log.rr,'C:/Users/dmw63/Desktop/ili.log.rr.csv')

n.days <- 52
select.indices <- (length(dates1)-n.days):length(dates1)
dates<-dates1[select.indices]
states <- dimnames(pred)[[2]]
ages <- dimnames(pred)[[3]]

Observed weekly proportion of visits due to ILI vs seasonal baseline (+/-95%PI)

  par(mfrow=c(4,4))

dates <- result.object$dates
states <- dimnames(pred)[[2]]
week.select.index<-select.indices

plot.state.rank <- cbind.data.frame(state.index=1:dim(rr)[2],state.rank= rank(-rr[dim(rr)[1],,1]))
plot.state.rank <- plot.state.rank[order(plot.state.rank$state.rank),]
plot.state.indices <- plot.state.rank$state.index

for(i in plot.state.indices){
  yrange1<-range(c(pred[week.select.index,i,1]/denom[week.select.index,i,1],obs[week.select.index,i,1]/denom[week.select.index,i,1],lpi[week.select.index,i,1]/denom[week.select.index,i,1],upi[week.select.index,i,1]/denom[week.select.index,i,1] ))

  plot(dates[week.select.index],
       pred[week.select.index,i,1]/denom[week.select.index,i,1],
       type='l',
       col='red',
       bty='l',
       ylim=yrange1,
       ylab='Proportion ILI',
       xlab='Date',
       main=states[i])

  points(dates[week.select.index],
         obs[week.select.index,i,1]/denom[week.select.index,i,1],
         type='l',
         col='black')

  polygon(c(dates[week.select.index],
            rev(dates[week.select.index])),
          c(lpi[week.select.index,i,1]/denom[week.select.index,i,1],
            rev(upi[week.select.index,i,1]/denom[week.select.index,i,1])),
          col = rgb(1, 0, 0, alpha = 0.1),
          border = NA)
}

date.mmwrdates <- mmwr_week(dates1)
mmwr.epiyr<- date.mmwrdates$mmwr_year
mmwr.epiyr[date.mmwrdates$mmwr_week<=26] <- mmwr.epiyr[date.mmwrdates$mmwr_week<=26] - 1

mmwr.epiwk <- date.mmwrdates$mmwr_week
mmwr.epiwk[date.mmwrdates$mmwr_week>=27]<-date.mmwrdates$mmwr_week[date.mmwrdates$mmwr_week>=27] - 52
mmwr.epiwk <- mmwr.epiwk +26
check<-cbind.data.frame(date.mmwrdates,mmwr.epiwk, mmwr.epiyr)

Observed ILI%/expected ILI% by state

  par(mfrow=c(4,4))
rr2<-rr[,,1]
for(i in plot.state.indices){
    y.range1<-c(0,5)
    ds2<-cbind.data.frame('epiwk'=mmwr.epiwk,'epiyr'=mmwr.epiyr, rr=rr2[,i])
    ds2.c<-dcast(ds2, epiwk~epiyr, value.var='rr', fun.aggregate = mean)
    cols1<-c(rep('grey',(ncol(ds2.c)-2) ),'red')
    matplot(ds2.c$epiwk         ,
       exp(ds2.c[,-1]),
       type='l',
       col=cols1,
       ylim=y.range1,
       bty='l',
       lty=1,
       ylab='Observed/Expected',
       main=paste(states[i]))
    abline(h=1, col='black')
    }

dashboardPlotOe(excess_cases1,
                            datevar='date',  
                            agevar='age',
                            statevar='state',  
                            outcome='ilitotal', 
                            yaxis='state',
                            facet='age')

weinbergerlab/ExcessILI documentation built on May 30, 2021, 10:57 a.m.