algo.glrnb: Count Data Regression Charts

In jimhester/surveillance: Temporal and Spatio-Temporal Modeling and Monitoring of Epidemic Phenomena

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 range vector of indices in the observed vector to monitor (should be consecutive) mu0 A vector of in-control values of the mean of the negative binomial distribution with the same length as range. If NULL the observed values in 1:(min(range)-1) are used to estimate beta through a generalized linear model. To fine-tune the model one can instead specify mu0 as a list with two components: S number of harmonics to include trend include a term t in the GLM model alpha The (known) dispersion parameter of the negative binomial distribution. If alpha=0 then the negative binomial distribution boils down to the Poisson distribution and a call of algo.glrnb is equivalent to a call to algo.glrpois. If alpha=NULL the parameter is calculated as part of the in-control estimation. c.ARL threshold in the GLR test, i.e. c_gamma Mtilde number of observations needed before we have a full rank the typical setup for the "intercept" and "epi" charts is Mtilde=1 M number of time instances back in time in the window-limited approach, i.e. the last value considered is \max{1,n-M}. To always look back until the first observation use M=-1. change a string specifying the type of the alternative. Currently the two choices are intercept and epi. See the SFB Discussion Paper 500 for details. theta if NULL then the GLR scheme is used. If not NULL the prespecified value for κ or λ is used in a recursive LR scheme, which is faster. dir a string specifying the direction of testing in GLR scheme. With "inc" only increases in x are considered in the GLR-statistic, with "dec" decreases are regarded. ret a string specifying the type of upperbound-statistic that is returned. With "cases" the number of cases that would have been necessary to produce an alarm or with "value" the glr-statistic is computed (see below).

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(... >= c_gamma)

where instead of 1<= k <= n the GLR statistic is computed for all k \in {n-M, …, n-Mtilde+1}. To achieve the typical behaviour from 1<= k <= 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,… compared to time 1,2,… 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.

Value

survRes

algo.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.

Author(s)

M. Hoehle

Source

Count data regression charts for the monitoring of surveillance time series (2008), M. Höhle 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, M., SFB386 Discussion Paper 500.

See Also

algo.rkiLatestTimepoint

Examples

##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)

jimhester/surveillance documentation built on May 19, 2019, 10:33 a.m.