algo.glrnb {surveillance} | R Documentation |
Count data regression charts for the monitoring of surveillance time series.
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")))
disProgObj |
object of class disProg to do surveillance
for |
control |
A list controlling the behaviour of the algorithm
|
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.
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. |
M. Hoehle
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.
##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)