# R/excess_model.R In excessmort: Excess Mortality

#### Documented in excess_model

#' Fit excess count model
#'
#' This function estimates the increase in the rate for a count time series relative to
#' the rate for a typical year. Two options are available: 1 - model the rate increase as a
#' smooth function and estimate this function or 2 - estimate the total excess in intervals.
#' For 1m an event date can be provided and a discontinuity included in the model.
#' You can do either 1 or 2 or both.
#'
#' Three versions of the model are available: 1 - Assume counts are Poisson distributed,
#' 2 - assume counts are overdispersed Poisson, or 3 - assume a mixed model with
#' correlated errors. The second is the default and recommended for weekly count data. For daily counts we often find
#' evidence of correlation and recommend the third along with setting weekday.effect = TRUE.
#'
#' If the counts object includes a expected column produced by compute_expected these are used
#' as the expected counts. If not, then these are computed.
#'
#'
#' @param counts A data frame with date, count and population columns.
#' @param start First day of interval to which model will be fit
#' @param end Last day of interval to which model will be fit
#' @param knots.per.year Number of knots per year used for the fitted smooth function
#' @param event If modeling a discontinuity is desired, this is the day in which it happens
#' @param intervals Instead of start and end a list of time intervals can be provided and excess is computed in each one
#' @param discontinuity Logical that determines if discontinuity is allowed at event
#' @param model Which version of the model to fit
#' @param exclude Dates to exclude when computing expected counts
#' @param trend.knots.per.year Number of knots per year used by compute_expected to estimate the trend for the expected counts
#' @param harmonics Number of harmonics used by compute_expected to estimate seasonal trend
#' @param frequency Number of observations per year. If not provided an attempt is made to calculate it
#' @param weekday.effect Logical that determins if a day of the week effects is included in the model. Should be FALSE for weekly or monthly data.
#' @param control.dates When model is set to correlated, these dates are used to estimate the covariance matrix. The larger this is the slower the function runs.
#' @param max.control If the length of control.dates is larger than max.control the function stops.
#' @param order.max Larges order for the Autoregressive process used to model the covariance structure
#' @param aic A logical that determines if the AIC criterion is used to selected the order of the AR process
#' @param maxit Maxium number of iterations for the IRLS algorithm used when model is correlated
#' @param epsilon Difference in deviance requried to declare covergenace of IRLS
#' @param alpha Percentile used to define what is outside the normal range
#' @param min.rate The estimated expected rate is not permited to go below this value
#' @param verbose Logical that determines if messages are displayed
#'
#' @return If only intervals are provided a data frame with excess estimates described below for excess.
#' if start and end are provided the a list with the following components are included:
#' \describe{
#' \item{date}{The dates for which the estimate was computed}
#' \item{observed}{The observed counts}
#' \item{expected}{The expected counts}
#' \item{fitted}{The fitted curve for excess counts}
#' \item{se}{The point-wise standard error for the fitted curve}
#' \item{population}{The population size}
#' \item{sd}{The standard deviation for observed counts on a typical year}
#' \item{cov}{The estimated covariance matrix for the observed counts}
#' \item{x}{The design matrix used for the fit}
#' \item{betacov}{The covariance matrix for the estimated coefficients}
#' \item{dispersion}{The estimated overdispersion parameter}
#' \item{detected_intervals}{Time intervals for which the 1 - alpha confidence interval does not include 0}
#' \item{ar}{The estimated coefficients for the autoregressive process}
#' \item{excess}{A data frame with information for the time intervals provided in itervals. This includes start, end, observed death rate (per 1,000 per year), expected death rate, standard deviation for the death rate, observed counts, expected counts, excess counts, standard deviation}
#' }
#'
#' @examples
#' data(cdc_state_counts)
#' counts <- cdc_state_counts[cdc_state_counts$state == "Massachusetts", ] #' exclude_dates <- c(seq(as.Date("2017-12-16"), as.Date("2018-01-16"), by = "day"), #' seq(as.Date("2020-01-01"), max(cdc_state_counts$date), by = "day"))
#' f <- excess_model(counts,
#' exclude = exclude_dates,
#' start = min(counts$date), #' end = max(counts$date),
#' knots.per.year = 12)

#' 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) #' #' @export #' @importFrom stats ARMAacf glm poly qnorm fitted.values excess_model <- function(counts, start = NULL, end = NULL, knots.per.year = 12, event = NULL, intervals = NULL, discontinuity = TRUE, model = c("quasipoisson", "poisson", "correlated"), exclude = NULL, trend.knots.per.year = 1/7, harmonics = 2, frequency = NULL, weekday.effect = FALSE, control.dates = NULL, max.control = 5000, order.max = 14, aic = TRUE, maxit = 25, epsilon = 1e-8, alpha = 0.05, min.rate = 0.0001, verbose = TRUE){ if("compute_expected" %in% class(counts)){ if(attr(counts, "keep.components")){ counts <- counts$counts
}
} else{
if(verbose) message("Computing expected counts.")
counts <-  compute_expected(counts,
exclude = exclude,
trend.knots.per.year = trend.knots.per.year,
harmonics = harmonics,
frequency = frequency,
weekday.effect = weekday.effect,
keep.components = FALSE,
verbose = verbose)
}

correlated.errors <- match.arg(model) == "correlated"

## number of observations per year
frequency <- attr(counts, "frequency")
dispersion <- attr(counts, "dispersion")

## checks
if(any(counts$excluded)) exclude <- counts$date[counts$excluded] else exclude <- NULL if(correlated.errors & is.null(control.dates)){ if(!is.null(exclude)){ warning("No control region suplied, which is not recommended when correlated.errors = TRUE. Using data up to first excluded point.") control.dates <-seq(min(counts$date), min(exclude, na.rm = TRUE) - 1, by = "day")
} else{
warning("No control region suplied, which is not recommended when correlated.errors = TRUE. Using all the data")
control.dates <- counts$date } } if(length(control.dates) > max.control & correlated.errors) stop("Length of control longer than", max.control) if((is.null(start) & !is.null(end)) | (!is.null(start) & is.null(end))) stop("You must provide both start and end, not just one.") if(is.null(start) & is.null(end) & is.null(intervals)) stop("You must provide start and end or intervals.") ## check to see if counts per unit of time are high enough for model to work if(mean(counts$outcome, na.rm = TRUE) < 1 & correlated.errors)
warning("Low counts per unit of time, consider fitting model with no correlation.")

## compute_expected always uses quasipoisson
if(match.arg(model) == "poisson") dispersion <- 1

## Use control days to compute the autocorrelation function
if(correlated.errors){
arfit <- fit_ar(counts, control.dates, order.max = order.max, aic = aic)
s <- arfit$sigma if(verbose) message("Order selected for AR model is ",length(arfit$ar),
". Estimated residual standard error is ", signif(arfit$sigma, 2), ".") } ## Fit start and end provided, fit the curve if(!is.null(start) & !is.null(end)){ if(!is.null(event)){ if(event >= end | event <= start) stop("event must be between start and end.") } ## now fit the GLS to the relevant subset of data include_dates = seq(start, end, by = "day") if(frequency == 365) include_dates <- include_dates[!(lubridate::month(include_dates) == 2 & lubridate::day(include_dates) == 29)] ind <- which(counts$date %in% include_dates)
date <- counts$date[ind] n <- length(ind) x <- 0:(n-1) / frequency * 365 mu <- pmax(min.rate, counts$expected[ind])
if(any(mu == min.rate)) warning("Minimum expected rate reached and was set at ", min.rate, ".")
min.fhat <- min.rate / mu - 1
obs <- counts$outcome[ind] pop <- counts$population[ind]

## compute residuals to fit ar model
if(correlated.errors) y <- (obs - mu) / mu

## create the design matrix
nknots <- round(knots.per.year * as.numeric(max(date) - min(date)) / 365)
knots <- x[round(seq(1, n, length = nknots + 2))]
knots <- knots[-c(1, length(knots))]
if(!is.null(event)){
event_index <- x[which.min(abs(as.numeric(date - event)))]
i <- which.min(abs(knots - event_index))
##shift knots so that one of the internal knots falls on the event day
knots <- knots + (event_index -  knots[i])
X <- cbind(1, splines::ns(x, knots = knots))
## add parameters to account for discontinuity
if(discontinuity){
after_ind <- as.numeric(I(x >= event_index))
X <- cbind(X, after_ind, poly((x - event_index)*after_ind, degree = 2))
}
} else{
X <- cbind(1, splines::ns(x, knots = knots))
}

if(correlated.errors){
fhat <- 0
beta <- 0;
dev<- 2*sum(ifelse(obs == 0, 0, obs * log(obs / mu)) - (obs - mu))
## convinience function
mysolve <- function(x) chol2inv(chol(x))

## parameters for covariance matrix
if(length(arfit$ar) > 0 & arfit$sigma > 0){
rhos <- ARMAacf(ar = arfit$ar, ma = 0, lag.max = n) } ## start iterations count <- 0 flag <- TRUE log_mu_vari <- counts$log_expected_se[ind]^2
while(count < maxit & flag){
if(length(arfit$ar) > 0 & s > 0){ Sigma <- apply(abs(outer(1:n, 1:n, "-")) + 1, 1, function(i) rhos[i]) * outer(sqrt((1 + fhat)^2 * s^2 + (1 + fhat)/mu + (1 + fhat)^2 * log_mu_vari), sqrt((1 + fhat)^2 * s^2 + (1 + fhat)/mu + (1 + fhat)^2 * log_mu_vari)) Sigma_inv <- mysolve(Sigma) } else{ Sigma <- diag((1 + fhat)^2 * s^2 + (1 + fhat)/mu + (1 + fhat)^2 * log_mu_vari) Sigma_inv <- diag(1/((1 + fhat)^2 * s^2 + (1 + fhat)/mu + (1 + fhat)^2 * log_mu_vari)) } ## fit spline using weighted least squares xwxi <- mysolve(t(X) %*% Sigma_inv %*% X) beta <- xwxi %*% t(X) %*% Sigma_inv %*% y count <- count + 1 fhat <- pmax(as.vector(X %*% beta), min.fhat) devold <- dev dev <- 2*sum(ifelse(obs == 0, 0, obs * log(obs / (mu*(1 + fhat)))) - (obs - mu*(1 + fhat))) flag <- abs(dev - devold)/(0.1 + abs(dev)) >= epsilon } if(count > maxit) warning("No convergence after ", maxit, " imterations.") se <- sqrt(apply(X, 1, function(x) matrix(x, nrow = 1) %*% xwxi %*% matrix(x, ncol = 1))) betacov <- xwxi } else{ fit <- glm(obs ~ X-1, offset = log(mu), family = "poisson") tmp <- predict(fit, se = TRUE, type = "response") fhat <- pmax(tmp$fit/mu - 1, min.fhat)
lambda <- fitted.values(fit)
cova <- summary(fit)$cov.unscaled * summary(fit)$dispersion
lambda_vari <- lambda^2 * diag(X %*% cova %*% t(X))
mu_vari <- counts$expected[ind]^2 * counts$log_expected_se[ind]^2
se <- sqrt((lambda_vari / mu^2) + (lambda^2 * mu_vari / mu^4))
Sigma <- diag(n) *  summary(fit)$dispersion / mu betacov <- summary(fit)$cov.unscaled * summary(fit)$dispersion } ## Warning if minimum reached if(any(fhat == min.fhat)) warning("Minimum rate reached and was set so that estimated rate is ", min.rate, ".") ## Compute regions for which estimate was outside usual range excess_ind <- which(fhat - qnorm(1 - alpha/2) * se >= 0) if(length(excess_ind) > 0){ cluster <- cumsum(c(2, diff(excess_ind)) > 1) indexes <- split(excess_ind, cluster) excess <- lapply(indexes, function(i){ n <- length(i) excess_stats(min(counts$date[ind[i]]),
max(counts$date[ind[i]]), counts$outcome[ind[i]],
counts$expected[ind[i]], Sigma[i, i, drop = FALSE], counts$population[ind[i]],
frequency,
fhat[i],
X[i,,drop=FALSE],
betacov)
})
detected_intervals <- do.call(rbind, excess)
} else{
detected_intervals <- data.frame(start = NA, end = NA, total = 0, se = NA, natural= NA)
}
ret <- list(date = date,
observed = obs,
expected = mu,
log_expected_se = counts$log_expected_se[ind], fitted = fhat, se = se, population = pop, sd = sqrt(diag(Sigma)), cov = Sigma, x = X, betacov = betacov, dispersion = dispersion, detected_intervals = detected_intervals) attr(ret, "frequency") <- frequency attr(ret, "model") <- match.arg(model) attr(ret, "type") <- "curve_fit" if(correlated.errors){ ret$ar <-  arfit
}

}

## If intervals provided compute excess deaths
## The uncertainty is calculated under the null
if(!is.null(intervals)){
res <- lapply(intervals, function(dates){
ind         <- which(counts$date %in% dates) date <- counts$date[ind]
n           <- length(ind)
mu          <- counts$expected[ind] log_mu_vari <- counts$log_expected_se[ind]^2

if(correlated.errors){
if(length(arfit$ar) > 0 & arfit$sigma > 0){
rhos <- ARMAacf(ar = arfit$ar, ma = 0, lag.max = n) } if(length(arfit$ar) > 0 & arfit$sigma > 0){ Sigma <- apply(abs(outer(1:n, 1:n, "-")) + 1, 1, function(i) rhos[i]) * outer(sqrt(s^2 + 1/mu + log_mu_vari), sqrt(s^2 + 1/mu + log_mu_vari)) } else{ Sigma <- diag(s^2 + 1/mu + log_mu_vari) } } else{ Sigma <- diag(n) * dispersion / mu } excess_stats(min(dates), max(dates), counts$outcome[ind],
mu,
Sigma,
counts$population[ind], frequency) }) res <- do.call(rbind, res) if(!exists("ret")){ ret <- res attr(ret, "type") <- "excess" } else{ ret$excess <- res
attr(ret, "type") <- append(attr(ret, "type"), "excess")
}
}
return(ret)
}


## 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.