Nothing
## ----SETUP, include = FALSE--------------------------------------------------------
## create directories for plots and cache
dir.create("plots", showWarnings=FALSE)
dir.create("monitoringCounts-cache", showWarnings=FALSE)
## load packages
library('surveillance')
library('gamlss')
## ----echo=FALSE--------------------------------------------------------------------
data("salmNewport")
## ----echo=FALSE--------------------------------------------------------------------
stopifnot(
all.equal(observed(salmNewport),
observed(as(as(salmNewport, "ts"), "sts")))
)
## ----echo=FALSE--------------------------------------------------------------------
# This code is the one used for the Salmon et al. (2016) JSS article.
# Using this code all examples from the article can be reproduced.
# computeALL is FALSE to avoid the computationally intensive parts
# of the code (simulations to find a threshold value for categoricalCUSUM,
# INLA-driven BODA) but one can set it to TRUE to have it run.
computeALL <- FALSE
## ----NewportPlot-simple, fig.keep = 'none'-----------------------------------------
plot(salmNewport, type = observed ~ time,
xaxis.tickFreq = list("%m" = atChange, "%G" = atChange),
xaxis.labelFreq = list("%Y" = atMedian), xaxis.labelFormat = "%Y")
## ----unitPlot-simple, echo = FALSE, fig.keep = 'none'------------------------------
plot(salmNewport, units = 2:3)
## ----EARS, fig.keep='none'---------------------------------------------------------
in2011 <- which(isoWeekYear(epoch(salmNewport))$ISOYear == 2011)
salmNewportGermany <- aggregate(salmNewport, by = "unit")
control <- list(range = in2011, method = "C1", alpha = 0.05)
surv <- earsC(salmNewportGermany, control = control)
plot(surv)
## ----farHead-----------------------------------------------------------------------
con.farrington <- list(
range = in2011, noPeriods = 1,
b = 4, w = 3, weightsThreshold = 1,
pastWeeksNotIncluded = 3, pThresholdTrend = 0.05,
thresholdMethod = "delta"
)
con.noufaily <- list(
range = in2011, noPeriods = 10,
b = 4, w = 3, weightsThreshold = 2.58,
pastWeeksNotIncluded = 26, pThresholdTrend = 1,
thresholdMethod = "nbPlugin"
)
## ----echo=F------------------------------------------------------------------------
con.farrington$limit54 <- con.noufaily$limit54 <- c(0,50) # for the figure
## ----------------------------------------------------------------------------------
salm.farrington <- farringtonFlexible(salmNewportGermany, con.farrington)
salm.noufaily <- farringtonFlexible(salmNewportGermany, con.noufaily)
## ----farPlot-simple, echo = FALSE, fig.keep = 'none'-------------------------------
par(mfrow = c(1,2))
plot(salm.farrington)
plot(salm.noufaily)
## ----campyDE-simple, fig.keep='none'-----------------------------------------------
data("campyDE")
cam.sts <- sts(epoch = campyDE$date,
observed = campyDE$case, state = campyDE$state)
plot(cam.sts, col = "mediumblue")
lines(campyDE$hum * 50, col = "white", lwd = 2)
axis(4, at = seq(0, 2500, by = 500), labels = seq(0, 50, by = 10))
## ----boda-cache, echo = FALSE, results='hide'--------------------------------------
if (computeALL) {
##hoehle 2018-07-18: changed code to use NICELOOKINGboda, but that's iid. Reason:
##The option 'rw1' currently crashes INLA.
library("INLA")
rangeBoda <- which(epoch(cam.sts) >= as.Date("2007-01-01"))
control.boda <- list(range = rangeBoda, X = NULL, trend = TRUE,
season = TRUE, prior = "iid", alpha = 0.025,
mc.munu = 10000, mc.y = 1000,
samplingMethod = "marginals")
boda <- boda(cam.sts, control = control.boda)
save(list = c("boda", "control.boda", "rangeBoda"),
file = "monitoringCounts-cache/boda.RData")
} else {
load("monitoringCounts-cache/boda.RData")
}
## ----boda2-cache, echo = FALSE, results='hide'-------------------------------------
if (computeALL) {
covarNames <- c("l1.hum", "l2.hum", "l3.hum", "l4.hum",
"newyears", "christmas", "O104period")
control.boda2 <- modifyList(control.boda,
list(X = campyDE[, covarNames], season = FALSE))
boda.covars <- boda(cam.sts, control = control.boda2)
save(list = c("boda.covars", "covarNames", "control.boda2"),
file = "monitoringCounts-cache/boda.covars.RData")
} else {
load("monitoringCounts-cache/boda.covars.RData")
}
## ----bPlot-simple, echo = FALSE, fig.keep = 'none'---------------------------------
plot(boda.covars)
## ----boda3, echo = FALSE-----------------------------------------------------------
control.far <- list(range=rangeBoda,b=4,w=5,alpha=0.025*2)
far <- farrington(cam.sts,control=control.far)
#Both farringtonFlexible and algo.bayes uses a one-sided interval just as boda.
control.far2 <-modifyList(control.far,list(alpha=0.025))
farflex <- farringtonFlexible(cam.sts,control=control.far2)
bayes <- suppressWarnings(bayes(cam.sts,control=control.far2))
## ----boda4, echo = FALSE-----------------------------------------------------------
# Small helper function to combine several equally long univariate sts objects
combineSTS <- function(stsList) {
epoch <- as.numeric(epoch(stsList[[1]]))
observed <- NULL
alarm <- NULL
for (i in 1:length(stsList)) {
observed <- cbind(observed,observed(stsList[[i]]))
alarm <- cbind(alarm,alarms(stsList[[i]]))
}
colnames(observed) <- colnames(alarm) <- names(stsList)
res <- sts(epoch=as.numeric(epoch), epochAsDate=TRUE,
observed=observed, alarm=alarm)
return(res)
}
## ----alarmplot, fig.width=8, fig.height=4, out.width="\\linewidth", echo=FALSE-----
# Make an artificial object containing two columns - one with the boda output
# and one with the farrington output
cam.surv <- combineSTS(list(boda.covars=boda.covars,boda=boda,bayes=bayes,
farrington=far,farringtonFlexible=farflex))
par(mar=c(4,8,2.1,2),family="Times")
plot(cam.surv,type = alarm ~ time,lvl=rep(1,ncol(cam.surv)),
alarm.symbol=list(pch=17, col="red2", cex=1,lwd=3),
cex.axis=1,xlab="Time (weeks)",cex.lab=1,xaxis.tickFreq=list("%m"=atChange,"%G"=atChange),xaxis.labelFreq=list("%G"=at2ndChange),
xaxis.labelFormat="%G")
## ----glrnb, results='hide'---------------------------------------------------------
phase1 <- which(isoWeekYear(epoch(salmNewportGermany))$ISOYear < 2011)
phase2 <- in2011
control <- list(range = phase2, c.ARL = 4, theta = log(2), ret = "cases",
mu0 = list(S = 1, trend = TRUE, refit = FALSE))
salmGlrnb <- glrnb(salmNewportGermany, control = control)
## ----cat---------------------------------------------------------------------------
data("salmHospitalized")
isoWeekYearData <- isoWeekYear(epoch(salmHospitalized))
dataBefore2013 <- which(isoWeekYearData$ISOYear < 2013)
data2013 <- which(isoWeekYearData$ISOYear == 2013)
dataEarly2014 <- which(isoWeekYearData$ISOYear == 2014
& isoWeekYearData$ISOWeek <= 4)
phase1 <- dataBefore2013
phase2 <- c(data2013, dataEarly2014)
salmHospitalized.df <- cbind(as.data.frame(salmHospitalized),
weekNumber = isoWeekYearData$ISOWeek)
names(salmHospitalized.df) <- c("y", "t", "state", "alarm", "upperbound", "n",
"freq", "epochInPeriod", "weekNumber")
## ----catbis, results='hide'--------------------------------------------------------
vars <- c( "y", "n", "t", "epochInPeriod", "weekNumber")
m.bbin <- gamlss(cbind(y, n-y) ~ 1 + t
+ sin(2 * pi * epochInPeriod) + cos(2 * pi * epochInPeriod)
+ sin(4 * pi * epochInPeriod) + cos(4 * pi * epochInPeriod)
+ I(weekNumber == 1) + I(weekNumber == 2),
sigma.formula =~ 1,
family = BB(sigma.link = "log"),
data = salmHospitalized.df[phase1, vars])
## ----cat2--------------------------------------------------------------------------
R <- 2
h <- 2
pi0 <- predict(m.bbin, newdata = salmHospitalized.df[phase2, vars],
type = "response")
pi1 <- plogis(qlogis(pi0) + log(R))
pi0m <- rbind(pi0, 1 - pi0)
pi1m <- rbind(pi1, 1 - pi1)
## ----cat2bis-----------------------------------------------------------------------
populationHosp <- unname(cbind(
population(salmHospitalized),
population(salmHospitalized)))
observedHosp <- cbind(
"Yes" = as.vector(observed(salmHospitalized)),
"No" = as.vector(population(salmHospitalized) - observed(salmHospitalized)))
salmHospitalized.multi <- sts(
frequency = 52, start = c(2004, 1), epoch = epoch(salmHospitalized),
observed = observedHosp, population = populationHosp,
multinomialTS = TRUE)
## ----cat2ter-----------------------------------------------------------------------
dBB.cusum <- function(y, mu, sigma, size, log = FALSE) {
dBB(if (is.matrix(y)) y[1,] else y,
if (is.matrix(y)) mu[1,] else mu,
sigma = sigma, bd = size, log = log)
}
## ----cat3--------------------------------------------------------------------------
controlCat <- list(range = phase2, h = 2, pi0 = pi0m, pi1 = pi1m,
ret = "cases", dfun = dBB.cusum)
salmHospitalizedCat <- categoricalCUSUM(
salmHospitalized.multi, control = controlCat,
sigma = exp(m.bbin$sigma.coefficients))
## ----------------------------------------------------------------------------------
h.grid <- seq(1, 10, by = 0.5)
## ----cath-cache, echo = FALSE, results='hide'--------------------------------------
if (computeALL) {
simone <- function(sts, h) {
y <- rBB(length(phase2), mu = pi0m[1, , drop = FALSE],
bd = population(sts)[phase2, ],
sigma = exp(m.bbin$sigma.coefficients),
fast = TRUE)
observed(sts)[phase2, ] <- cbind(y, population(sts)[phase2, 1] - y)
one.surv <- categoricalCUSUM(
sts, control = modifyList(controlCat, list(h = h)),
sigma = exp(m.bbin$sigma.coefficients))
return(any(alarms(one.surv)[, 1]))
}
set.seed(123)
nSims <- 1000
pMC <- sapply(h.grid, function(h) {
mean(replicate(nSims, simone(salmHospitalized.multi, h)))
})
pMarkovChain <- sapply(h.grid, function(h) {
TA <- LRCUSUM.runlength(mu = pi0m[1,,drop = FALSE],
mu0 = pi0m[1,,drop = FALSE],
mu1 = pi1m[1,,drop = FALSE],
n = population(salmHospitalized.multi)[phase2, ],
h = h, dfun = dBB.cusum,
sigma = exp(m.bbin$sigma.coef))
return(tail(TA$cdf, n = 1))
})
save(pMC, file = "monitoringCounts-cache/pMC.RData")
save(pMarkovChain, file = "monitoringCounts-cache/pMarkovChain.RData")
} else {
load("monitoringCounts-cache/pMC.RData")
load("monitoringCounts-cache/pMarkovChain.RData")
}
## ----catF-simple, echo = FALSE, fig.keep = 'none'----------------------------------
plot(salmHospitalizedCat[,1])
## ----catARL-simple, echo = FALSE, fig.keep = 'none'--------------------------------
matplot(h.grid, cbind(pMC, pMarkovChain), type="l", lty=1:2, col=1)
abline(h=0.1, lty=5)
legend("center", c("Monte Carlo","Markov chain"), lty=1:2, bty="n")
## ----rotaPlot-simple, fig.keep='none'----------------------------------------------
data("rotaBB")
plot(rotaBB)
## ----------------------------------------------------------------------------------
rotaBB.df <- as.data.frame(rotaBB)
X <- with(rotaBB.df, cbind(intercept = 1, epoch,
sin1 = sin(2 * pi * epochInPeriod),
cos1 = cos(2 * pi * epochInPeriod)))
phase1 <- epoch(rotaBB) < as.Date("2009-01-01")
phase2 <- !phase1
library("MGLM")
## MGLMreg automatically takes the last class as ref so we reorder
order <- c(2:5, 1); reorder <- c(5, 1:4)
m0 <- MGLMreg(as.matrix(rotaBB.df[phase1, order]) ~ -1 + X[phase1, ],
dist = "MN")
## ----------------------------------------------------------------------------------
m1 <- m0
m1@coefficients[1, ] <- m0@coefficients[1, ] + log(2)
pi0 <- t(predict(m0, newdata = X[phase2, ])[, reorder])
pi1 <- t(predict(m1, newdata = X[phase2, ])[, reorder])
## ----CATCUSUM----------------------------------------------------------------------
dfun <- function(y, size, mu, log = FALSE) {
dmultinom(x = y, size = size, prob = mu, log = log)
}
h <- 2 # threshold for the CUSUM statistic
control <- list(range = seq(nrow(rotaBB))[phase2], h = h, pi0 = pi0,
pi1 = pi1, ret = "value", dfun = dfun)
surv <- categoricalCUSUM(rotaBB,control=control)
## ----include = FALSE---------------------------------------------------------------
alarmDates <- epoch(surv)[which(alarms(surv)[,1])]
format(alarmDates,"%b %Y")
## ----CATCUSUMMC,echo=FALSE,eval=FALSE----------------------------------------------
# #Number of MC samples
# nSamples <- 1e4
#
# #Do MC
# simone.stop <- function(sts, control) {
# phase2Times <- seq(nrow(sts))[phase2]
# #Generate new phase2 data from the fitted in control model
# y <- sapply(1:length(phase2Times), function(i) {
# rmultinom(n=1, prob=pi0[,i],size=population(sts)[phase2Times[i],1])
# })
# observed(sts)[phase2Times,] <- t(y)
# one.surv <- categoricalCUSUM(sts, control=control)
# #compute P(S<=length(phase2))
# return(any(alarms(one.surv)[,1]>0))
# }
#
# set.seed(1233)
# rlMN <- replicate(nSamples, simone.stop(rotaBB, control=control))
# mean(rlMN) # 0.5002
## ----------------------------------------------------------------------------------
m0.dm <- MGLMreg(as.matrix(rotaBB.df[phase1, 1:5]) ~ -1 + X[phase1, ],
dist = "DM")
c(m0@AIC, m0.dm@AIC)
## ----------------------------------------------------------------------------------
## Change intercept in the first class (for DM all 5 classes are modeled)
delta <- 2
m1.dm <- m0.dm
m1.dm@coefficients[1, ] <- m0.dm@coefficients[1, ] +
c(-delta, rep(delta/4, 4))
alpha0 <- exp(X[phase2,] %*% m0.dm@coefficients)
alpha1 <- exp(X[phase2,] %*% m1.dm@coefficients)
dfun <- function(y, size, mu, log = FALSE) {
dLog <- ddirmn(t(y), t(mu))
if (log) dLog else exp(dLog)
}
h <- 2
control <- list(range = seq(nrow(rotaBB))[phase2], h = h,
pi0 = t(alpha0), pi1 = t(alpha1),
ret = "value", dfun = dfun)
surv.dm <- categoricalCUSUM(rotaBB, control = control)
## ----echo=FALSE,eval=FALSE---------------------------------------------------------
# matplot(alpha0/rowSums(alpha0),type="l",lwd=3,lty=1,ylim=c(0,1))
# matlines(alpha1/rowSums(alpha1),type="l",lwd=1,lty=2)
## ----ctPlot-simple, echo = FALSE, fig.keep = 'none'--------------------------------
par(mfrow = c(1,2))
surv@multinomialTS <- surv.dm@multinomialTS <- FALSE # trick plot method ...
plot(surv[,1], col=c(NA,NA,4), ylab = expression(C[t]), ylim = c(0,33),
xaxis.tickFreq=list("%Y"=atChange, "%m"=atChange),
xaxis.labelFreq=list("%Y"=atMedian), xaxis.labelFormat="%Y")
abline(h=h, lwd=2, col="darkgrey")
plot(surv.dm[,1], col=c(NA,NA,4), ylab = expression(C[t]), ylim = c(0,33),
xaxis.tickFreq=list("%Y"=atChange, "%m"=atChange),
xaxis.labelFreq=list("%Y"=atMedian), xaxis.labelFormat="%Y")
abline(h=h, lwd=2, col="darkgrey")
## ----------------------------------------------------------------------------------
today <- which(epoch(salmNewport) == as.Date("2013-12-23"))
rangeAnalysis <- (today - 4):today
in2013 <- which(isoWeekYear(epoch(salmNewport))$ISOYear == 2013)
algoParameters <- list(range = rangeAnalysis, noPeriods = 10,
populationBool = FALSE,
b = 4, w = 3, weightsThreshold = 2.58,
pastWeeksNotIncluded = 26, pThresholdTrend = 1,
thresholdMethod = "nbPlugin", alpha = 0.05,
limit54 = c(0, 50))
results <- farringtonFlexible(salmNewport[, c("Baden.Wuerttemberg",
"North.Rhine.Westphalia")],
control = algoParameters)
## ----results='asis'----------------------------------------------------------------
start <- isoWeekYear(epoch(salmNewport)[min(rangeAnalysis)])
end <- isoWeekYear(epoch(salmNewport)[max(rangeAnalysis)])
caption <- paste0("Results of the analysis of reported S. Newport ",
"counts in two German federal states for the weeks ",
start$ISOYear, "-W", start$ISOWeek, " to ",
end$ISOYear, "-W", end$ISOWeek,
". Bold red counts indicate weeks with alarms.")
toLatex(results, caption = caption, label = "tableResults",
ubColumnLabel = "Threshold", include.rownames = FALSE,
sanitize.text.function = identity)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.