algo.hmm  R Documentation 
This function implements online HMM detection of outbreaks based on
the retrospective procedure described in Le Strat and Carret (1999).
Using the function msm
(from package msm) a specified HMM
is estimated, the decoding problem, i.e. the most probable state
configuration, is found by the Viterbi algorithm and the most
probable state of the last observation is recorded. Online
detection is performed by sequentially repeating this procedure.
Warning: This function can be very slow  a more efficient implementation would be nice!
algo.hmm(disProgObj, control = list(range=range, Mtilde=1,
noStates=2, trend=TRUE, noHarmonics=1,
covEffectEqual=FALSE, saveHMMs = FALSE, extraMSMargs=list()))
disProgObj 
object of class disProg (including the observed and the state chain) 
control 
control object:

For each time point t the reference values values are extracted. If the number of requested values is larger than the number of possible values the latter is used. Now the following happens on these reference values:
A noStates
State Hidden Markov Model (HMM) is used based on
the Poisson distribution with linear predictor on the loglink
scale. I.e.
Y_t  X_t = j \sim Po(\mu_t^j),
where
\log(\mu_t^j) = \alpha_j + \beta_j\cdot
t + \sum_{i=1}^{nH} \gamma_j^i \cos(2i\pi/freq\cdot (t1)) +
\delta_j^i \sin(2i\pi/freq\cdot (t1))
and nH=
noHarmonics
and freq=12,52
depending on the
sampling frequency of the surveillance data. In the above t1
is
used, because the first week is always saved as t=1
, i.e. we
want to ensure that the first observation corresponds to cos(0) and
sin(0).
If covEffectEqual
then all covariate effects parameters are
equal for the states, i.e. \beta_j=\beta, \gamma_j^i=\gamma^i,
\delta_j^i=\delta^i
for all j=1,...,\code{noStates}
.
In case more complicated HMM models are to be fitted it is possible to
modify the msm
code used in this function. Using
e.g. AIC
one can select between different models (see the
msm package for further details).
Using the Viterbi algorithms the most probable state configuration
is obtained for the reference values and if the most probable
configuration for the last reference value (i.e. time t) equals
control$noOfStates
then an alarm is given.
Note: The HMM is refitted from scratch every time, sequential updating schemes of the HMM would increase speed considerably! A major advantage of the approach is that outbreaks in the reference values are handled automatically.
algo.hmm
gives a list of class survRes
which includes the
vector of alarm values for every timepoint in range
. No
upperbound
can be specified and is put equal to zero.
The resulting object contains a list control$hmms
, which
contains the "msm"
objects with the fitted HMMs
(if saveHMMs=TRUE
).
M. Höhle
Y. Le Strat and F. Carrat, Monitoring Epidemiologic Surveillance Data using Hidden Markov Models (1999), Statistics in Medicine, 18, 3463–3478
I.L. MacDonald and W. Zucchini, Hidden Markov and Other Models for Discretevalued Time Series, (1997), Chapman & Hall, Monographs on Statistics and applied Probability 70
msm
#Simulate outbreak data from HMM
set.seed(123)
counts < sim.pointSource(p = 0.98, r = 0.8, length = 3*52,
A = 1, alpha = 1, beta = 0, phi = 0,
frequency = 1, state = NULL, K = 1.5)
## Not run:
#Do surveillance using a two state HMM without trend component and
#the effect of the harmonics being the same in both states. A sliding
#window of two years is used to fit the HMM
surv < algo.hmm(counts, control=list(range=(2*52):length(counts$observed),
Mtilde=2*52,noStates=2,trend=FALSE,
covEffectsEqual=TRUE,extraMSMargs=list()))
plot(surv,legend.opts=list(x="topright"))
## End(Not run)
if (require("msm")) {
#Retrospective use of the function, i.e. monitor only the last time point
#but use option saveHMMs to store the output of the HMM fitting
surv < algo.hmm(counts,control=list(range=length(counts$observed),Mtilde=1,noStates=2,
trend=FALSE,covEffectsEqual=TRUE, saveHMMs=TRUE))
#Compute most probable state using the viterbi algorithm  1 is "normal", 2 is "outbreak".
viterbi.msm(surv$control$hmms[[1]])$fitted
#How often correct?
tab < cbind(truth=counts$state + 1 ,
hmm=viterbi.msm(surv$control$hmm[[1]])$fitted)
table(tab[,1],tab[,2])
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.