surveillance (version 1.11.0)

bodaDelay: Bayesian aberration detection in presence of reporting delays.

Description

The function takes range values of the surveillance time series sts and for each time point uses a Bayesian model of the negative binomial family with log link inspired by the work of Noufaily et al. (2012) and of Manitz and H�hle{Hoehle} (2014). It allows delay-corrected aberration detection as explained in Salmon et al. (2015). A reportingTriangle has to be provided in the control slot.

Usage

bodaDelay(sts, control = list(range = NULL, b = 3, w = 3,
                                          mc.munu=100, mc.y=10,
                                          pastAberrations = TRUE, 
                                          verbose = FALSE,
                                          alpha = 0.01, trend = TRUE,
                                          limit54=c(5,4), 
                                          inferenceMethod=c("asym","INLA"),
                                          quantileMethod=c("MC","MM"),
                                          noPeriods = 1, 
                                          pastWeeksNotIncluded = 26,
                                          delay = TRUE))

Arguments

sts

sts-object to be analysed. Needs to have a reporting triangle.

control

list with control arguments

b

How many years back in time to include when forming the base counts.

w

Window's half-size, i.e. number of weeks to include before and after the current week in each year.

range

Specifies the index of all timepoints which should be tested. If range is NULL all possible timepoints are used.

pastAberrations

Boolean indicating whether to include an effect for past outbreaks in a second fit of the model. This option only makes sense if inferenceMethod is INLA, as it is not supported by the other inference method.

verbose

Boolean specifying whether to show extra debugging information.

alpha

An approximate (one-sided) $(1-\alpha)\cdot 100%$ prediction interval is calculated unlike the original method where it was a two-sided interval. The upper limit of this interval i.e. the $(1-\alpha)\cdot 100%$ quantile serves as an upperbound.

trend

Boolean indicating whether a trend should be included

noPeriods

Number of levels in the factor allowing to use more baseline. If equal to 1 no factor variable is created, the set of reference values is defined as in Farrington et al (1996).

inferenceMethod

Which inference method used, as defined in Salmon et al. (2015). If one chooses INLA then inference is performed with INLA. If one chooses asym then the asymptotic normal approximation of the posteriori is used.

pastWeeksNotIncluded

Number of past weeks to ignore in the calculation.

delay

Boolean indicating whether to take reporting delays into account.

mc.munu

Number of samples for the parameters of the negative binomial distribution for calculating a threshold

mc.y

Number of samples for observations when performing Monte Carlo to calculate a threshold

limit54

c(cases,period) is a vector allowing the user to change these numbers.

quantileMethod

Character, either MC or MM. Indicates how to compute the quantile based on the posterior distribution (no matter the inference method): either by sampling mc.munu values from the posterior distribution of the paramet

References

Farrington, C.P., Andrews, N.J, Beale A.D. and Catchpole, M.A. (1996): A statistical algorithm for the early detection of outbreaks of infectious disease. J. R. Statist. Soc. A, 159, 547-563. Noufaily, A., Enki, D.G., Farrington, C.P., Garthwaite, P., Andrews, N.J., Charlett, A. (2012): An improved algorithm for outbreak detection in multiple surveillance systems. Statistics in Medicine, 32 (7), 1206-1222. Salmon, M., Schumacher, D., Stark, K., H�hle{Hoehle}, M. (2015): Bayesian outbreak detection in the presence of reporting delays. Biometrical Journal, 57 (6), 1051-1067.

Examples

data(stsNewport)
salm.Normal <- list()
salmDelayAsym <- list()
for (week in 43:45){
  listWeeks <- as.Date(row.names(stsNewport@control$reportingTriangle$n), origin="1970-01-01")
  dateObs <- listWeeks[isoWeekYear(listWeeks)$ISOYear==2011&
                         isoWeekYear(listWeeks)$ISOWeek==(week)]
  stsC <- sts_observation(stsNewport,
                          dateObservation=dateObs,
                          cut=TRUE)
  inWeeks <- which(formatDate(epoch(stsC), "                     isoWeekYear(as.Date(epoch(stsC)))$ISOWeek>=40 & 
                     isoWeekYear(as.Date(epoch(stsC)))$ISOWeek<=48)
  
  rangeTest <- inWeeks
  alpha <- 0.07
  #  isoWeekYear(epoch(dataTrunc))$ISOWeek >39 & 
  # isoWeekYear(epoch(dataTrunc))$ISOWeek <=47)
  # Control slot for Noufaily method          
  controlNoufaily <- list(range=rangeTest,noPeriods=10,
                          b=4,w=3,weightsThreshold=2.58,pastWeeksNotIncluded=26,
                          pThresholdTrend=1,thresholdMethod="nbPlugin",alpha=alpha*2,
                          limit54=c(0,50))
  # Control slot for the Proposed algorithm with D=0 correction
  controlNormal <- list(range = rangeTest, b = 4, w = 3,
                        reweight = TRUE, mc.munu=10000, mc.y=100,
                        verbose = FALSE,
                        alpha = alpha, trend = TRUE,
                        limit54=c(0,50), 
                        noPeriods = 10, pastWeeksNotIncluded = 26,
                        delay=FALSE)
  # Control slot for the Proposed algorithm with D=10 correction
  
  controlDelayNorm <-  list(range = rangeTest, b = 4, w = 3,
                            reweight = FALSE, mc.munu=10000, mc.y=100,
                            verbose = FALSE,
                            alpha = alpha, trend = TRUE,
                            limit54=c(0,50), 
                            noPeriods = 10, pastWeeksNotIncluded = 26,
                            delay=TRUE,inferenceMethod="asym")
  set.seed(1)
  
  salm.Normal[[week]] <- farringtonFlexible(stsC, controlNoufaily)
  
  salmDelayAsym[[week]] <- bodaDelay(stsC, controlDelayNorm)
  
}

par(mfrow=c(2,3))
lapply(salmDelayAsym[c(43,44,45)],plot, legend=NULL, main="", ylim=c(0,35))
lapply(salm.Normal[c(43,44,45)],plot, legend=NULL, main="", ylim=c(0,35))