algo.glrpois {surveillance} | R Documentation |
Poisson regression charts for the monitoring of surveillance time series.
algo.glrpois(disProgObj,control = list(range=range,c.ARL=5, mu0=NULL, Mtilde=1, M=-1, change="intercept",theta=NULL,dir="inc"))
disProgObj |
object of class disProg to do surveillance
for |
control |
A list controlling the behaviour of the algorithm
|
This function implements the seasonal Poisson charts based on generalized likelihood ratio (GLR) as described in the SFB Discussion Paper 500. 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.glrpois 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. |
M. Hoehle
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} 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 <- rpois(length(t),mu0*(exp(kappa)^(t>=tau))) 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 cntrl = list(range=t,c.ARL=5, Mtilde=1, mu0=mu0,change="intercept",dir="inc") glr.ts <- algo.glrpois(s.ts,control=cntrl) lr.ts <- algo.glrpois(s.ts,control=c(cntrl,theta=0.4)) plot(glr.ts,xaxis.years=FALSE) plot(lr.ts,xaxis.years=FALSE)