Nothing
## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.align = "center",
fig.width = 7,
fig.height = 5
)
## ----data---------------------------------------------------------------------
onset <- as.Date(c("2017-02-04", "2017-02-12", "2017-02-15",
"2017-02-23", "2017-03-01", "2017-03-01",
"2017-03-02", "2017-03-03", "2017-03-03"))
## ----incidence----------------------------------------------------------------
library(incidence)
i <- incidence(onset)
i
plot(i, border = "white")
## ----incidence2---------------------------------------------------------------
today <- as.Date("2017-03-21")
i <- incidence(onset, last_date = today)
i
plot(i, border = "white")
## ----si-----------------------------------------------------------------------
mu <- 15.3 # mean in days days
sigma <- 9.3 # standard deviation in days
## ----estimate-----------------------------------------------------------------
library(earlyR)
library(ggplot2)
res <- get_R(i, si_mean = mu, si_sd = sigma)
res
plot(res)
plot(res, "lambdas", scale = length(onset) + 1) +
geom_vline(xintercept = onset, col = "grey", lwd = 1.5) +
geom_vline(xintercept = today, col = "blue", lty = 2, lwd = 1.5)
## ----sample_R-----------------------------------------------------------------
R_val <- sample_R(res, 1000)
summary(R_val)
quantile(R_val)
quantile(R_val, c(0.025, 0.975))
hist(R_val, border = "grey", col = "navy",
xlab = "Values of R",
main = "Sample of likely R values")
## ----generate_si--------------------------------------------------------------
si <- res$si
si
## ----projections--------------------------------------------------------------
library(projections)
future_i <- project(i, R = R_val, n_sim = 1000, si = res$si, n_days = 30)
future_i
mean(future_i) # average incidence / day
plot(future_i)
## ----pred_size----------------------------------------------------------------
predicted_n <- colSums(future_i)
summary(predicted_n)
hist(predicted_n, col = "darkred", border = "white",
main = "Prediction: new cases in 30 days",
xlab = "Total number of new cases")
## ----alternative--------------------------------------------------------------
alt_i <- incidence(onset)
alt_res <- get_R(alt_i, si_mean = mu, si_sd = sigma)
alt_R_val <- sample_R(alt_res, 1000)
alt_future_i <- project(alt_i, R = alt_R_val, n_sim = 1000, si = res$si, n_days = 30)
alt_future_i
mean(alt_future_i)
plot(alt_future_i)
## alternative plot
col <- "#cc66991a"
matplot(alt_future_i, type = "l", col = col, lty = 1, lwd = 5,
xlab = "Day from today",
ylab = "Projected daily incidence")
alt_predicted_n <- colSums(alt_future_i)
summary(alt_predicted_n)
hist(alt_predicted_n, col = "darkred", border = "white",
main = "Prediction: new cases in 30 days",
xlab = "Total number of new cases")
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.