# R/excess_cumulative.R In excessmort: Excess Mortality

#### Documented in excess_cumulative

#' Compute cumulative excess deaths
#'
#' This function takes the output of the excess_model function, a start date, and
#' end date and calculates excess mortality and standard errors.
#'
#' @param fit The output of excess_model
#' @param start The start date
#' @param end The end date
#'
#' @return A data frame with the following columns
#' \describe{
#' \item{date}{The date}
#' \item{observed}{The observed excess mortality,which is the sum of observed minus expected until that date}
#' \item{sd}{The standard deviation for excess mortality for that date if year is typical}
#' \item{fitted}{The estimated of excess mortality based on the estimated smooth event effect curve}
#' \item{se}{The standard error for fitted}
#' }
#'
#' @examples
#'
#' data(new_jersey_counts)
#' exclude_dates <- as.Date("2012-10-29") + 0:180
#' control_dates <- seq(min(new_jersey_counts$date), min(exclude_dates) - 1, by = "day") #' f <- excess_model(new_jersey_counts, #' start = as.Date("2012-09-01"), #' end = as.Date("2013-09-01"), #' exclude = exclude_dates, #' model = "correlated", #' weekday.effect = TRUE, #' control.dates = control_dates) #' #' excess_cumulative(f, #' start = as.Date("2017-12-15"), #' end = as.Date("2017-12-21") ) #' #' @export excess_cumulative <- function(fit, start, end){ if(!"curve_fit" %in% attr(fit, "type")) stop("This is not the correct excess_model fit, needs curve fit.") ind <- which(fit$date %in% seq(start, end, by = "day"))
n               <- length(ind)
A               <- matrix(1, n, n)
A[upper.tri(A)] <- 0
A               <- sweep(A, 2, fit$expected[ind], FUN = "*") fhat <- matrix(fit$fitted[ind], ncol = 1)

fit_excess   <- A %*% fhat
obs_excess   <- cumsum(fit$observed[ind] - fit$expected[ind])
varcov       <- fit$x[ind,] %*% fit$betacov %*% t(fit$x[ind,]) diag(varcov) <- fit$se[ind]^2
fitted_se    <- sqrt(diag(A %*%  varcov %*% t(A)))
sd           <- sqrt(diag(A %*% (fit$cov[ind,ind] + diag(fit$log_expected_se[ind]^2)) %*% t(A)))
data.frame(date     = fit\$date[ind],
observed = obs_excess,
sd       = sd,
fitted   = fit_excess,
se       = fitted_se)
}


## Try the excessmort package in your browser

Any scripts or data that you put into this service are public.

excessmort documentation built on Oct. 11, 2021, 9:06 a.m.