Last chance! 50% off unlimited learning
Sale ends in
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"), xMax=1e4))algo.glrpois(disProgObj, control = list(range=range, c.ARL=5,
mu0=NULL, Mtilde=1, M=-1, change="intercept",
theta=NULL, dir=c("inc","dec"),
ret=c("cases","value"), xMax=1e4))
disProg to do surveillance
foralgo.glrpois simply calls algo.glrnb with
control$alpha set to 0.
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 time point
$<=n$. both="" lead="" to="" the="" same="" alarm="" timepoints,="" but="" 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.
In case is.null(theta) and alpha>0 as well as
ret="cases" then a brute-force search is conducted for each time
point in range in order to determine the number of cases necessary
before an alarm is sounded. In case no alarm was sounded so far by time
$t$, the function increases $x[t]$ until an alarm is sounded any
time before time point $t$. If no alarm is sounded by xMax, a return value
of 1e99 is given. Similarly, if an alarm was sounded by time $t$ the
function counts down instead. Note: This is slow experimental code!
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.
Salmon, M., Schumacher, D. and
##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)
## Poisson example
#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",ret="value",dir="inc")
glr.ts <- algo.glrpois(s.ts,control=cntrl)
plot(glr.ts,xaxis.years=FALSE)
lr.ts <- algo.glrpois(s.ts,control=c(cntrl,theta=0.4))
plot(lr.ts,xaxis.years=FALSE)
## NegBin example
#Generate data
set.seed(42)
alpha <- 0.2
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