Interpolate demographic data

Description

Interpolate yearly population estimates so that a population estimate is provided for each day of the year. The function 'approx' is used with 'rule = 2' for extrapolation.

Usage

approx_demographics(
  demo,
  first_day,
  last_day,
  by,
  extrapolation.type = c("linear", "constant", "none"),
  ...
)

Arguments

Value

A data frame with dates and population estimates.

Weekly death counts for each USA state

Description

The Center for Disease Control (CDC) provides weekly estimated death counts for each state in the USA. This data frame includes these estimates for each state along with population sizes.

• state Name of USA state

• date Corresponding date of observation

• outcome Estimated number of deaths

• population Population estimate (from the Census)

Usage

data(cdc_state_counts)

Format

An object of class "data.frame"

Callapse age groups into broader ones

Description

Collapse a count or demographic data frame into a broader age strata.

Usage

collapse_age_dist(demo, breaks)

collapse_counts_by_age(counts, breaks)

Arguments

Value

A age groups represented in 'demo' or 'counts' are grouped using the new age breaks defined by breaks but containing the populations and counts, if applicable, for age groups defined by 'breaks'.

Examples

library(lubridate)
data(cook_records)
## define smaller subset for example
cook_demographics_subset <- cook_demographics[year(cook_demographics$date)==2021, ]
demo <- collapse_age_dist(cook_demographics_subset, 
                          breaks = c(0, 20, 40, 60, 80, Inf))

Compute counts

Description

Compute counts for groups from individual records

Usage

compute_counts(
  dat,
  by = NULL,
  demo = NULL,
  date = "date",
  age = "age",
  agegroup = "agegroup",
  breaks = NULL
)

Arguments

Details

This is helper function that helps convert individual records data, in which each death is a row, to a count data frame where each row is a date. It is particulalry helpful for defining agegroups. If you provide the 'breaks' argument it will automaticall divided the data and provide the counts for each age strata. You can also select variables to group by using the 'by' argument. One can provide a data frame with demogrpahic inform through the 'demo' argument. This tabe must have the population size for each group for each data.

Value

A data frame with counts for each group for each date with population sizes, if demo was provided.

Examples

library(lubridate)
data("cook_records")
the_breaks <- c(0, 20, 40, 60, 75, Inf)

## take subset for example
cook_records_subset <- cook_records[year(cook_records$date)==2021, ]
cook_demographics_subset <- cook_demographics[year(cook_demographics$date)==2021, ]

cook_counts <- compute_counts(cook_records_subset, 
                       demo = cook_demographics_subset, 
                       by = c("agegroup", "race", "sex"), 
                       breaks = the_breaks)

Compute expected counts for each day

Description

Compute the expected death count for each unit of time. We assume counts are over-dispersed Poisson distributed with a trend that accounts for slow, year-to-year, changes in death rate across time and a seasonal effect. The function takes a data frame with dates and counts and returns the data frame with the expected counts as a new column. It also returns a logical column that is 'TRUE' if that entry was used in the estimation procedure.

Usage

compute_expected(
  counts,
  exclude = NULL,
  include.trend = TRUE,
  trend.knots.per.year = 1/7,
  extrapolate = TRUE,
  harmonics = 2,
  frequency = NULL,
  weekday.effect = FALSE,
  keep.components = TRUE,
  verbose = TRUE
)

Arguments

Details

Periods for which excess deaths will be estimated should be excluded when estimating expected counts. These can be supplied via the 'exclude' argument. Note that If 'extrapolate' is 'TRUE', the default, the time trend will be extrapolated following the estimated trend. If 'extrapolate' is 'FALSE' the trend is assumed to be a constant equal to the estimate on the last day before extrapolation. If the period for which excess deaths are estimated is long, extrapolation should be used with caution. We highly recommend exploring the estimated expected counts with the 'expected_diagnostics' function.

Value

The 'counts' data.frame with two columns added: 'expected' and 'excluded'. The 'expected' column is the estimated expected value of the counts for that date. The 'excluded' column is a logical vector denoting if that date was excluded when estimating the expected value.

If the argument 'keep.components' is 'TRUE' a list is returned with 'counts' data.frame in the first component, the estimated trend in the second, the estimated seasonal effect in the third and the estimated weekday effects in the fourth.

Examples

data(new_jersey_counts)
exclude_dates <- as.Date("2012-10-29") + 0:180
counts <- compute_expected(new_jersey_counts, exclude = exclude_dates, weekday.effect = TRUE)
library(ggplot2)
expected_plot(counts)

Cook County Medical Examiner Records

Description

A data frame with each row represening a death. The data includes death date and demographic information. When you load this dataset you also load 'cook_demographics' which includes the population size for each demographice group by date.

Usage

data("cook_records")

Format

A data frame with these columns

sex

Sex of the deceased

age

Age of the deceased

race

Race of the deceased

residenceplace

Residence placed of the deceased

date

Date of the death

cause_1

Reported cause of death

type_of_death

Type of death

Compute cumulative excess deaths

Description

This function takes the output of the 'excess_model' function, a start date, and end date and calculates excess mortality and standard errors.

Usage

excess_cumulative(fit, start, end)

Arguments

Value

A data frame with the following columns

date

The date

observed

The observed excess mortality,which is the sum of observed minus expected until that date

sd

The standard deviation for excess mortality for that date if year is typical

fitted

The estimated of excess mortality based on the estimated smooth event effect curve

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") )

Fit excess count model

Description

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 1 an 'event' date can be provided and a discontinuity included in the model. You can do either 1 or 2 or both.

Usage

excess_model(
  counts,
  start = NULL,
  end = NULL,
  knots.per.year = 12,
  event = NULL,
  intervals = NULL,
  discontinuity = TRUE,
  model = c("quasipoisson", "poisson", "correlated"),
  exclude = NULL,
  include.trend = TRUE,
  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-08,
  alpha = 0.05,
  min.rate = 1e-04,
  keep.counts = FALSE,
  keep.components = TRUE,
  verbose = TRUE
)

Arguments

Details

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.

Value

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:

date

The dates for which the estimate was computed

observed

The observed counts

expected

The expected counts

fitted

The fitted curve for excess counts

se

The point-wise standard error for the fitted curve

population

The population size

sd

The standard deviation for observed counts on a typical year

cov

The estimated covariance matrix for the observed counts

x

The design matrix used for the fit

betacov

The covariance matrix for the estimated coefficients

dispersion

The estimated overdispersion parameter

detected_intervals

Time intervals for which the 1 - 'alpha' confidence interval does not include 0

ar

The estimated coefficients for the autoregressive process

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)

Plot results from fitted excess count model

Description

Plot results from fitted excess count model

Usage

excess_plot(fit, title = "", ylim = NULL, show.data = TRUE, alpha = 0.05)

Arguments

Value

An ggplot object containing the plot.

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)

library(ggplot2)
excess_plot(f)

Excess counts in an interval

Description

Helper function to compute excess deaths statistics for a

Usage

excess_stats(
  start,
  end,
  obs,
  mu,
  cov,
  pop,
  frequency,
  fhat = NULL,
  X = NULL,
  betacov = NULL
)

Arguments

Diagnostic Plots for Model Fit

Description

Check mean model fit via diagnostic figures of the model components

Usage

expected_diagnostic(
  expected,
  start = NULL,
  end = NULL,
  color = "#D22B2B",
  alpha = 0.5
)

Arguments

Value

A list with six ggplot objects: 'population' is a time series plot of the population. 'seasonal' is a plot showing the estimated seasonal effect. 'trend' is a time series plot showing the estimated trend. 'weekday' is a plot of the estimated weekday effects if they were estimated. 'expected'is a time series plot of the expected values. 'residual' is a time series plot of the residuals.

Examples

library(dplyr)
library(lubridate)
library(ggplot2)

flu_season <- seq(make_date(2017, 12, 16), make_date(2018, 1, 16), by = "day")

exclude_dates <- c(flu_season, seq(make_date(2020, 1, 1), today(), by = "day"))

res  <- cdc_state_counts %>%
 filter(state == "Massachusetts") %>%
 compute_expected(exclude = exclude_dates,
                  keep.components = TRUE)
                  
p <- expected_diagnostic(expected = res, alpha = 0.50)

p$population
p$seasonal
p$trend
p$expected
p$residual

Plot Expected Counts

Description

Check if expected counts fit data

Usage

expected_plot(
  expected,
  title = "",
  start = NULL,
  end = NULL,
  ylim = NULL,
  weekly = FALSE,
  color = "#3366FF",
  alpha = 0.5
)

Arguments

Value

A ggplot object containing a plot of the original counts and the estimated expected values.

Examples

data(new_jersey_counts)
exclude_dates <- as.Date("2012-10-29") + 0:180
e <- compute_expected(new_jersey_counts, exclude = exclude_dates, weekday.effect = TRUE)

library(ggplot2)
expected_plot(e, start = as.Date("2012-09-01"), end = as.Date("2013-09-01"))

Fit an ar model to residuals from expected counts

Description

Helper function to estimate autoregressive mode

Usage

fit_ar(counts, control.dates = NULL, order.max = 5, aic = FALSE, plot = FALSE)

Arguments

Get demographic data from Census

Description

Get demographic data from Census

Usage

get_demographics(
  geography = "state",
  state,
  county = NULL,
  years = 2018,
  vars = c("SEX", "AGEGROUP", "RACE", "HISP")
)

Arguments

Assign age to group

Description

Assign age to group

Usage

group_age(age, breaks)

Arguments

Value

A factor that groups the ages into the age groups defined by 'breaks'.

Louisiana daily mortality

Description

A data frame with Louisiana daily mortality counts from 2003-01-01 to 2016-12-31, which includes the day Katrina made landfall on 2015-08-29. The outcome column includes the number of deaths for that day.

Usage

data("louisiana_counts")

Format

An object of class tbl_df (inherits from tbl, data.frame) with 1461 rows and 3 columns.

New Jersey daily mortality

Description

A data frame with Louisiana daily mortality counts from 2007-01-01 to 2015-12-31 which includes the day Sandy made landfall on 2012-10-29. The outcome column includes the number of deaths for that day.

Usage

data("new_jersey_counts")

Format

An object of class tbl_df (inherits from tbl, data.frame) with 3287 rows and 3 columns.

Compute year of the day ingoring Feb 29

Description

Compute year of the day ingoring Feb 29

Usage

noleap_yday(x)

Arguments

Puerto Rico daily mortality

Description

A data frame with Puerto Rico daily mortality counts, stratified by agegroup, from 1985 to 2020. which includes the days hurricanes Hugo, Georges, and Maria made landfall on 1989-09-18, 1998-09-21, and 2017-09-20, respectively. The outcome column includes the number of deaths for that day. This data frame also includes population counts and estimates for the period. Population counts for 1985 to 2000 are interpolations of decennial census products since 1980 to 2000. For 2000 onward, we use interpolated population estimates (PEP).

Usage

data("puerto_rico_counts")

Format

An object of class data.table (inherits from data.frame) with 499644 rows and 5 columns.

Puerto Rico daily mortality by cause of death

Description

A data frame with Puerto Rico daily mortality counts, stratified by cause of death from 1999 to 2020. which includes the day hurricanes Maria made landfall on 2017-09-20. The outcome column includes the number of deaths for that day for that ICD-10 code. The object 'icd' is included to show the description of each ICD-10 code.

Usage

data("puerto_rico_icd")

Format

An object of class tbl_df (inherits from tbl, data.frame) with 262980 rows and 4 columns.

Weekly death counts for several countries

Description

A data frame with weekly death counts from multiple countries that includes the year 2020. Each country has a different range of data and most countries have at least 5 years of pre-2020 data. The original data was collated by the Financial Times .

Usage

data(world_counts)

Format

An object of class "data.frame"

Details

• country Name of the country

• date Corresponding date of observation

• outcome Estimated number of deaths

• population Population estimate (from the Census)