Learn R Programming

surveillance (version 1.5-4)

algo.glrnb: Cound data regression charts

Description

Count data regression charts for the monitoring of surveillance time series.

Usage

algo.glrnb(disProgObj,control = list(range=range,c.ARL=5, 
         mu0=NULL, alpha=0, Mtilde=1, M=-1, change="intercept",
         theta=NULL,dir=c("inc","dec"),ret=c("cases","value")))

Arguments

disProgObj
object of class disProg to do surveillance for
control
A list controlling the behaviour of the algorithm [object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object]

Value

  • survResalgo.glrnb returns a list of class survRes (surveillance result), which includes the alarm value for recognizing an outbreak (1 for alarm, 0 for no alarm), the threshold value for recognizing the alarm and the input object of class disProg. The upperbound slot of the object are filled with the current $GLR(n)$ value or with the number of cases that are necessary to produce an alarm at any timpoint $<=n$. both="" lead="" to="" the="" same="" alarm="" timepoints,="" but="" "cases" has an obvious interpretation.

encoding

latin1

source

Count data regression charts for the monitoring of surveillance time series (2008), M. H�hle{Hoehle} and M. Paul, Computational Statistics and Data Analysis, 52(9), pp. 4357--4368.

Poisson regression charts for the monitoring of surveillance time series (2006), H�hle{Hoehle}, M., SFB386 Discussion Paper 500.

Details

This function implements the seasonal cound data chart based on generalized likelihood ratio (GLR) as described in the Hoehle and Paul (2008) paper. A moving-window generalized likelihood ratio detector is used, i.e. the detector has the form$$N = \inf\left{ n : \max_{1\leq k \leq n} \left[ \sum_{t=k}^n \log \left{ \frac{f_{\theta_1}(x_t|z_t)}{f_{\theta_0}(x_t|z_t)} \right} \right] \geq c_\gamma \right}$$where instead of $1\leq k \leq n$ the GLR statistic is computed for all $k \in {n-M, \ldots, n-\tilde{M}+1}$. To achieve the typical behaviour from $1\leq k\leq n$ use Mtilde=1 and M=-1.

So $N$ is the time point where the GLR statistic is above the threshold the first time: An alarm is given and the surveillance is resetted starting from time $N+1$. Note that the same c.ARL as before is used, but if mu0 is different at $N+1,N+2,\ldots$ compared to time $1,2,\ldots$ the run length properties differ. Because c.ARL to obtain a specific ARL can only be obtained my Monte Carlo simulation there is no good way to update c.ARL automatically at the moment. Also, FIR GLR-detectors might be worth considering.

At the moment, window limited ``intercept'' charts have not been extensively tested and are at the moment not supported. As speed is not an issue here this doesn't bother too much. Therefore, a value of M=-1 is always used in the intercept charts.

See Also

algo.rkiLatestTimepoint

Examples

Run this code
##Simulate data and apply the algorithm
S <- 1 ; t <- 1:120 ; m <- length(t)
beta <- c(1.5,0.6,0.6)
omega <- 2*pi/52
#log mu_{0,t}
alpha <- 0.2
base <- beta[1] + beta[2] * cos(omega*t) + beta[3] * sin(omega*t) 
#Generate example data with changepoint and tau=tau
tau <- 100
kappa <- 0.4
mu0 <- exp(base)
mu1 <- exp(base  + kappa) 

#Generate data
set.seed(42)
x <- rnbinom(length(t),mu=mu0*(exp(kappa)^(t>=tau)),size=1/alpha)
s.ts <- create.disProg(week=1:length(t),observed=x,state=(t>=tau))

#Plot the data
plot(s.ts,legend=NULL,xaxis.years=FALSE)

#Run GLR based detection
cntrl = list(range=t,c.ARL=5, Mtilde=1, mu0=mu0, alpha=alpha, 
             change="intercept",ret="value",dir="inc")
glr.ts <- algo.glrnb(s.ts,control=c(cntrl))
plot(glr.ts,xaxis.years=FALSE)

#CUSUM LR detection with backcalculated number of cases
cntrl2 = list(range=t,c.ARL=5, Mtilde=1, mu0=mu0, alpha=alpha, 
              change="intercept",ret="cases",dir="inc",theta=1.2)
glr.ts2 <- algo.glrnb(s.ts,control=c(cntrl2))
plot(glr.ts2,xaxis.years=FALSE)

Run the code above in your browser using DataLab