Skip to content

Commit

Permalink
Fix broken quantiles (#88)
Browse files Browse the repository at this point in the history 
Closes #88
  • Loading branch information
@jstockwin
jstockwin committed Mar 17, 2020
1 parent c41d6e8 commit 8a5c3c4
Showing 1 changed file with 33 additions and 33 deletions.
66 changes: 33 additions & 33 deletions R/estimate_r.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,8 +34,8 @@
#' pairs of infector/infected individuals to be used to estimate the serial
#' interval distribution (see details).
#'
#' @param config An object of class \code{estimate_R_config}, as returned by
#' function \code{make_config}.
#' @param config An object of class \code{estimate_R_config}, as returned by
#' function \code{make_config}.
#'
#' @return {
#' an object of class \code{estimate_R}, with components:
Expand Down Expand Up @@ -104,7 +104,7 @@
#' different time windows, hence avoiding to rerun the MCMC every time
#' estimate_R is called.}
#' }
#'
#'
#' @seealso \code{\link{discr_si}} \code{\link{make_config}}
#' @author Anne Cori \email{a.cori@imperial.ac.uk}
#' @references {
Expand All @@ -125,30 +125,30 @@
#'
#' ## estimate the reproduction number (method "non_parametric_si")
#' ## when not specifying t_start and t_end in config, they are set to estimate
#' ## the reproduction number on sliding weekly windows
#' res <- estimate_R(incid = Flu2009$incidence,
#' ## the reproduction number on sliding weekly windows
#' res <- estimate_R(incid = Flu2009$incidence,
#' method = "non_parametric_si",
#' config = make_config(list(si_distr = Flu2009$si_distr)))
#' plot(res)
#'
#' ## the second plot produced shows, at each each day,
#' ## the estimate of the reproduction number over the 7-day window
#' ## the estimate of the reproduction number over the 7-day window
#' ## finishing on that day.
#'
#'
#' ## to specify t_start and t_end in config, e.g. to have biweekly sliding
#' ## windows
#' t_start <- seq(2, nrow(Flu2009$incidence)-13)
#' t_end <- t_start + 13
#' res <- estimate_R(incid = Flu2009$incidence,
#' ## windows
#' t_start <- seq(2, nrow(Flu2009$incidence)-13)
#' t_end <- t_start + 13
#' res <- estimate_R(incid = Flu2009$incidence,
#' method = "non_parametric_si",
#' config = make_config(list(
#' si_distr = Flu2009$si_distr,
#' t_start = t_start,
#' si_distr = Flu2009$si_distr,
#' t_start = t_start,
#' t_end = t_end)))
#' plot(res)
#'
#' ## the second plot produced shows, at each each day,
#' ## the estimate of the reproduction number over the 14-day window
#' ## the estimate of the reproduction number over the 14-day window
#' ## finishing on that day.
#'
#' ## example with an incidence object
Expand Down Expand Up @@ -241,7 +241,7 @@
#' R_si_from_sample <- estimate_R(MockRotavirus$incidence,
#' method = "si_from_sample",
#' si_sample = si_sample,
#' config = make_config(list(n2 = 50,
#' config = make_config(list(n2 = 50,
#' seed = overall_seed)))
#' plot(R_si_from_sample)
#'
Expand All @@ -260,12 +260,12 @@ estimate_R <- function(incid,
si_data = NULL,
si_sample = NULL,
config = make_config(incid = incid, method = method)) {

method <- match.arg(method)
config <- make_config(incid = incid, method = method, config = config)
config <- process_config(config)
check_config(config, method)

if (method == "si_from_data") {
## Warning if the expected set of parameters is not adequate
si_data <- process_si_data(si_data)
Expand Down Expand Up @@ -354,7 +354,7 @@ estimate_R_func <- function(incid,

calc_incidence_per_time_step <- function(incid, t_start, t_end) {
nb_time_periods <- length(t_start)
incidence_per_time_step <- vnapply(seq_len(nb_time_periods), function(i)
incidence_per_time_step <- vnapply(seq_len(nb_time_periods), function(i)
sum(incid[seq(t_start[i], t_end[i]), c("local", "imported")]))
return(incidence_per_time_step)
}
Expand All @@ -374,7 +374,7 @@ estimate_R_func <- function(incid,
b_posterior <- vector()
a_posterior <- vnapply(seq_len(nb_time_periods), function(t) if (t_end[t] >
final_mean_si) {
a_prior + sum(incid[seq(t_start[t], t_end[t]), "local"])
a_prior + sum(incid[seq(t_start[t], t_end[t]), "local"])
## only counting local cases on the "numerator"
}
else {
Expand Down Expand Up @@ -411,7 +411,7 @@ estimate_R_func <- function(incid,
b_posterior <- vector()
a_posterior <- vnapply(seq_len(nb_time_periods), function(t) if (t_end[t] >
final_mean_si) {
a_prior + sum(incid[seq(t_start[t], t_end[t]), "local"])
a_prior + sum(incid[seq(t_start[t], t_end[t]), "local"])
## only counting local cases on the "numerator"
}
else {
Expand All @@ -424,7 +424,7 @@ estimate_R_func <- function(incid,
else {
NA
})
sample_r_posterior <- vapply(seq_len(nb_time_periods), function(t)
sample_r_posterior <- vapply(seq_len(nb_time_periods), function(t)
if (!is.na(a_posterior[t])) {
rgamma(sample_size,
shape = unlist(a_posterior[t]),
Expand Down Expand Up @@ -532,11 +532,11 @@ estimate_R_func <- function(incid,
na.rm = TRUE
)
median_posterior <- apply(r_sample, 2, median, na.rm = TRUE)
quantile_0.25_posterior <- apply(r_sample, 2, quantile,
quantile_0.75_posterior <- apply(r_sample, 2, quantile,
0.75,
na.rm = TRUE
)
quantile_0.25_posterior <- apply(r_sample, 2, quantile,
quantile_0.95_posterior <- apply(r_sample, 2, quantile,
0.95,
na.rm = TRUE
)
Expand All @@ -550,10 +550,10 @@ estimate_R_func <- function(incid,
mean_si_sample <- rep(-1, config$n1)
std_si_sample <- rep(-1, config$n1)
for (k in seq_len(config$n1)) {
mean_si_sample[k] <- sum((seq_len(dim(si_sample)[1]) - 1) *
mean_si_sample[k] <- sum((seq_len(dim(si_sample)[1]) - 1) *
si_sample[, k])
std_si_sample[k] <- sqrt(sum(si_sample[, k] *
((seq_len(dim(si_sample)[1]) - 1) -
std_si_sample[k] <- sqrt(sum(si_sample[, k] *
((seq_len(dim(si_sample)[1]) - 1) -
mean_si_sample[k])^2))
}
temp <- lapply(seq_len(config$n1), function(k) sample_from_posterior(config$n2,
Expand All @@ -562,7 +562,7 @@ estimate_R_func <- function(incid,
b_prior, config$t_start, config$t_end
))
config$si_distr <- cbind(
t(vapply(seq_len(config$n1), function(k) (temp[[k]])[[2]],
t(vapply(seq_len(config$n1), function(k) (temp[[k]])[[2]],
numeric(nrow(si_sample)))),
rep(0, config$n1)
)
Expand All @@ -587,11 +587,11 @@ estimate_R_func <- function(incid,
na.rm = TRUE
)
median_posterior <- apply(r_sample, 2, median, na.rm = TRUE)
quantile_0.25_posterior <- apply(r_sample, 2, quantile,
quantile_0.75_posterior <- apply(r_sample, 2, quantile,
0.75,
na.rm = TRUE
)
quantile_0.25_posterior <- apply(r_sample, 2, quantile,
quantile_0.95_posterior <- apply(r_sample, 2, quantile,
0.95,
na.rm = TRUE
)
Expand Down Expand Up @@ -636,11 +636,11 @@ estimate_R_func <- function(incid,
shape = a_posterior,
scale = b_posterior, lower.tail = TRUE, log.p = FALSE
)
quantile_0.25_posterior <- qgamma(0.75,
quantile_0.75_posterior <- qgamma(0.75,
shape = a_posterior,
scale = b_posterior, lower.tail = TRUE, log.p = FALSE
)
quantile_0.25_posterior <- qgamma(0.95,
quantile_0.95_posterior <- qgamma(0.95,
shape = a_posterior,
scale = b_posterior, lower.tail = TRUE, log.p = FALSE
)
Expand All @@ -653,8 +653,8 @@ estimate_R_func <- function(incid,
results <- list(R = as.data.frame(cbind(
config$t_start, config$t_end, mean_posterior,
std_posterior, quantile_0.025_posterior, quantile_0.05_posterior,
quantile_0.25_posterior, median_posterior, quantile_0.25_posterior,
quantile_0.25_posterior, quantile_0.975_posterior
quantile_0.25_posterior, median_posterior, quantile_0.75_posterior,
quantile_0.95_posterior, quantile_0.975_posterior
)))

names(results$R) <- c(
Expand Down

0 comments on commit 8a5c3c4

Please sign in to comment.