## EXAMPLE 1 Foot and mouth disease, Cumbria 2001:
## Counts of incident FMD positive farms in Cumbria by date:
edate <- seq(from = as.Date("2001-02-28", format = "%Y-%m-%d"),
to = as.Date("2001-06-15", format = "%Y-%m-%d"), by = 1)
ncas <- c(1,2,0,0,1,1,0,0,4,2,3,3,5,2,8,2,5,0,5,7,15,13,6,7,11,8,7,11,
6,5,10,7,8,8,7,5,6,3,3,3,3,4,1,4,6,2,1,4,3,3,1,1,1,2,2,2,2,0,4,1,1,
0,0,2,1,2,1,0,1,2,2,4,0,1,0,1,0,0,0,1,0,3,1,1,3,0,0,1,1,2,0,0,1,1,1,
3,3,1,1,1,0,0,0,1,1,1,1,1)
dat.df01 <- data.frame(edate = edate, ncas = ncas)
## Seven day EDR:
edr <- epi.edr(dat.df01$ncas, n = 7, conf.level = 0.95, nsim = 99,
na.zero = FALSE)
dat.df01$edr.500 <- edr$est
dat.df01$edr.025 <- edr$lower
dat.df01$edr.975 <- edr$upper
## Smooth the EDRs:
dat.df01$sedr.500 <- lowess(dat.df01$edate, dat.df01$edr.500, f = 0.10)$y
dat.df01$sedr.025 <- lowess(dat.df01$edate, dat.df01$edr.025, f = 0.10)$y
dat.df01$sedr.975 <- lowess(dat.df01$edate, dat.df01$edr.975, f = 0.10)$y
## Not run:
if (FALSE) {
library(ggplot2); library(scales)
## Frequency histogram of the number of incident FMD positive farms as a
## function of date:
ggplot() +
theme_bw() +
geom_histogram(data = dat.df01, aes(x = edate, weight = ncas),
binwidth = 1, fill = "#008080", alpha = 0.45) +
scale_x_date(breaks = date_breaks("7 days"),
name = "Date of onset of signs") +
scale_y_continuous(limits = c(0,20), name = "Number of events") +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
## EDR line plot (and its 95
## date. Set primary axis breaks and labels:
x <- 2^(-5:5)
y <- seq(from = 0, to = 20, length = length(x))
pybreaks <- y; pylabels <- x
ggplot() +
theme_bw() +
geom_line(data = dat.df01, aes(x = edate,
y = 10 + (2 * log2(edr.500))), col = "black", linewidth = 1) +
geom_line(data = dat.df01, aes(x = edate,
y = 10 + (2 * log2(edr.025))), col = "black",
linetype = "dashed", linewidth = 0.5) +
geom_line(data = dat.df01, aes(x = edate,
y = 10 + (2 * log2(edr.975))), col = "black",
linetype = "dashed", linewidth = 0.5) +
scale_x_date(breaks = date_breaks("7 days"),
name = "Date of onset of signs") +
scale_y_continuous(breaks = pybreaks, labels = pylabels,
name = "Estimated dissemination ratio (EDR)") +
coord_cartesian(xlim = range(dat.df01$edate), ylim = c(0,20)) +
geom_hline(yintercept = 10, col = "red", linetype = "dashed") +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
## EDR line plot (and its 95
## date superimposed on the epidemic curve:
## Set the EDR values to appear on secondary vertical axis:
x <- 2^(-5:5)
## Set the number of FMD cases to appear on primary vertical axis:
y <- seq(from = 0, to = 20, length = length(x))
## Equation to use x (EDR estimate) to predict y (number of FMD cases):
dat.lm01 <- lm(y ~ log2(x))
summary(dat.lm01)
## Breaks and values for primary y axis:
pybreaks <- seq(from = 0, to = 20, by = 5)
pylabels <- pybreaks
## Breaks and values for secondary y axis:
sybreaks <- coefficients(dat.lm01)[1] + (coefficients(dat.lm01)[2] * log2(x))
sylabels <- x
ggplot() +
theme_bw() +
geom_histogram(data = dat.df01, aes(x = edate, weight = ncas),
binwidth = 1, fill = "#008080", alpha = 0.45) +
geom_line(data = dat.df01, aes(x = edate,
y = 10 + (2 * log2(edr.500))), col = "black", linewidth = 1) +
geom_line(data = dat.df01, aes(x = edate,
y = 10 + (2 * log2(edr.025))), col = "black", linewidth = 0.5,
linetype = "dashed") +
geom_line(data = dat.df01, aes(x = edate,
y = 10 + (2 * log2(edr.975))), col = "black", linewidth = 0.5,
linetype = "dashed") +
scale_x_date(breaks = date_breaks("1 month"),
labels = date_format("%b %y"), name = "Date of onset of signs") +
scale_y_continuous(breaks = pybreaks, labels = pylabels, limits = c(0,20),
name = "Number of FMD cases",
sec.axis = sec_axis(breaks = sybreaks, labels = sylabels,
transform = ~ ., name = "Estimated dissemination ratio (EDR)")) +
coord_cartesian(xlim = range(dat.df01$edate), ylim = c(0,20)) +
geom_hline(yintercept = 10, col = "red", linetype = "dashed")
}
Run the code above in your browser using DataLab