Type:	Package
Title:	Analyze Outbreak Models of Multi-Group Populations with Vaccination
Version:	0.1.1
Description:	Model infectious disease dynamics in populations with multiple subgroups having different vaccination rates, transmission characteristics, and contact patterns. Calculate final and intermediate outbreak sizes, form age-structured contact models with automatic fetching of U.S. census data, and explore vaccination scenarios with an interactive 'shiny' dashboard for a model with two subgroups, as described in Nguyen et al. (2024) <doi:10.1016/j.jval.2024.03.039> and Duong et al. (2026) <doi:10.1093/ofid/ofaf695.217>.
License:	MIT + file LICENSE
Encoding:	UTF-8
RoxygenNote:	7.3.3
Imports:	deSolve, graphics, shiny, stats, bslib (≥ 0.9.0), htmltools, socialmixr
Suggests:	knitr, rmarkdown, testthat (≥ 3.0.0),
Config/testthat/edition:	3
VignetteBuilder:	knitr
URL:	https://epiforesite.github.io/multigroup-vaccine/
Depends:	R (≥ 2.10)
NeedsCompilation:	no
Packaged:	2026-02-06 23:02:13 UTC; u0578572
Author:	Damon Toth [aut, cre, cph], Jake Wagoner [aut], Willy Ray [aut], George Vega Yon [ctb], Centers for Disease Control and Prevention's Center for Forecasting and Outbreak Analytics [fnd] (Cooperative agreement CDC-RFA-FT-23-0069)
Maintainer:	Damon Toth <damon.toth@hsc.utah.edu>
Repository:	CRAN
Date/Publication:	2026-02-09 20:10:05 UTC

Aggregate population counts into age groups

Description

Aggregate per-age population counts into coarser age groups defined by the (sorted) lower bounds in age_groups. When length(age_groups) == 1, all ages >= that value are aggregated into a single open-ended group ("Xplus"). When length(age_groups) > 1, groups are formed as age_groups[i] to age_groups[i+1] - 1 for i = 1:(n-1) and age_groups[n] and above for the final group. Human-readable labels are produced: "under1" for the 0–0 group, "ageX" for single-year groups, "XtoY" for ranges, and "Xplus" for the final open group. When verbose = TRUE, the function prints aggregation summaries to the console for each group using cat().

Usage

aggregateByAgeGroups(ages, pops, age_groups, verbose = FALSE)

Arguments

Details

• Group boundaries are inclusive at both ends for finite ranges (i.e. ages satisfying lower <= age <= upper). For the last group the upper bound is infinite.

• If no ages fall into a group the aggregated count for that group is 0 (because sum(numeric(0)) == 0).

• When verbose = TRUE, the function writes progress messages to the console with cat() for each aggregated group (useful for debugging / logging). By default (verbose = FALSE), the function is silent.

Value

A named list with components:

pops

Numeric vector of aggregated population counts, one element per group.

labels

Character vector of labels for each group (e.g. "under1", "age5", "0to4", "65plus").

age_ranges

List of numeric vectors of length 2 giving the inclusive lower and upper bounds for each group; the upper bound for the final group is Inf.

Examples

# Multiple groups example
ages <- 0:100
pops <- rep(100, length(ages))
aggregateByAgeGroups(ages, pops, c(0, 5, 18, 65))

# Single open-ended group (65plus)
aggregateByAgeGroups(ages, pops, 65)

Calculate a contact matrix for age groups and schools

Description

Calculate a contact matrix for age groups and schools

Usage

contactMatrixAgeSchool(
  agelims,
  agepops,
  schoolagegroups,
  schoolpops,
  schportion
)

Arguments

Value

a square matrix with the contact rate of each group (row) with members of each other group (column)

Examples

contactMatrixAgeSchool(agelims = c(0, 5, 18), agepops = c(500, 1300, 8200),
schoolagegroups = c(2, 2), schoolpops = c(600, 700), schportion = 0.7)

Calculate a contact matrix for age groups based on Polymod contact survey data

Description

Calculate a contact matrix for age groups based on Polymod contact survey data

Usage

contactMatrixPolymod(agelims, agepops = NULL)

Arguments

Details

The socialmixr contact_matrix function supports age limits up to 90. Any age limits above 90 will be adjusted to 90 with a warning, and the corresponding populations will be aggregated into a single "90+" group.

Value

A symmetric contact matrix with row and column names indicating the age groups.

Examples

#Default population distribution uses population data from POLYMOD survey locations:
contactMatrixPolymod(agelims = c(0, 5, 18))
#Specifying the age distribution will lead to an adjusted version:
contactMatrixPolymod(agelims = c(0, 5, 18), agepops = c(500, 1300, 8200))

Calculate group contact matrix with proportional mixing and preferential mixing within group

Description

Calculate group contact matrix with proportional mixing and preferential mixing within group

Usage

contactMatrixPropPref(popsize, contactrate, ingroup)

Arguments

Value

a square matrix with the contact rate of each group (row) with members of each other group (column)

Examples

contactMatrixPropPref(popsize = c(100, 150, 200), contactrate = c(1.1, 1, 0.9),
ingroup = c(0.2, 0.25, 0.22))

Disaggregate ACS 5-year age groups into single-year ages

Description

Takes population counts from ACS 5-year age groupings and uniformly distributes them into single-year ages. This allows for more flexible age group aggregations.

Usage

disaggregateCityAges(acs_age_pops)

Arguments

Value

A list containing:

Calculate final outbreak size or distribution of a multigroup transmission model for a given basic reproduction number, contact/transmission assumptions, and initial conditions

Description

Calculate final outbreak size or distribution of a multigroup transmission model for a given basic reproduction number, contact/transmission assumptions, and initial conditions

Usage

finalsize(
  popsize,
  R0,
  contactmatrix,
  relsusc,
  reltransm,
  initR,
  initI,
  initV,
  method = "ODE",
  nsims = 1
)

Arguments

Value

a vector (nsims = 1) or matrix (nsims > 1) with the final number infected from each group (column) in each simulation (row)

Examples

popsize <- c(800, 200)
R0 <- 2
contactmatrix <- contactMatrixPropPref(popsize = popsize, contactrate = c(1, 1),
ingroup = c(0.2, 0.2))
relsusc <- c(1, 1)
reltransm <- c(1, 1)
initR <- c(0, 0)
initI <- c(1, 0)
initV <- 0.2 * popsize
# Default method "ODE" numerical solves ordinary differential equations until infectious count
# is close to 0
finalsize(popsize, R0, contactmatrix, relsusc, reltransm, initR, initI, initV)
finalsize(popsize, R0, contactmatrix, relsusc, reltransm, initR, initI, initV,
method = "analytic")
finalsize(popsize, R0, contactmatrix, relsusc, reltransm, initR, initI, initV,
method = "stochastic", nsims = 10)
# All "escaped" outbreaks set to deterministic final size:
finalsize(popsize, R0, contactmatrix, relsusc, reltransm, initR, initI, initV,
method = "hybrid", nsims = 10)

Get Census Population Data by Age and County

Description

Downloads and processes U.S. Census Bureau population estimates for a specified state and county, organized by age groups. Supports single-year age data with optional sex disaggregation.

Usage

getCensusData(
  state_fips,
  county_name,
  year = 2024,
  age_groups = NULL,
  by_sex = FALSE,
  csv_path = NULL,
  cache_dir = NULL,
  verbose = FALSE
)

Arguments

Value

A list containing:

Examples

# Use the included example data (recommended for package examples)
slc_data <- getCensusData(
  state_fips = "49", 
  county_name = "Salt Lake County",
  year = 2024,
  csv_path = getCensusDataPath()
)

# Get age groups without sex disaggregation
slc_grouped <- getCensusData(
  state_fips = "49",
  county_name = "Salt Lake County", 
  year = 2024,
  age_groups = c(0, 5, 18, 65),
  csv_path = getCensusDataPath()
)

# Get age groups by sex
slc_by_sex <- getCensusData(
  state_fips = "49",
  county_name = "Salt Lake County",
  year = 2024, 
  age_groups = c(0, 5, 18, 65),
  by_sex = TRUE,
  csv_path = getCensusDataPath()
)


# Download from web (requires internet)
slc_web <- getCensusData(
  state_fips = "49",
  county_name = "Salt Lake County",
  year = 2024
)

# Use caching to avoid repeated downloads
slc_cached <- getCensusData(
  state_fips = "49",
  county_name = "Salt Lake County",
  year = 2024,
  cache_dir = "~/census_cache"
)

Get path to example census data file

Description

Returns the path to the example Utah census data CSV file included with the package. This is useful for examples, testing, and when internet access is not available.

Usage

getCensusDataPath()

Value

Character string with the path to the example census CSV file for Utah (FIPS 49)

Examples

# Get path to example Utah census file
utah_csv <- getCensusDataPath()

# Use it with getCensusData

slc_data <- getCensusData(
  state_fips = "49",
  county_name = "Salt Lake County",
  year = 2024,
  csv_path = getCensusDataPath()
)

Get City Population Data by Age

Description

Reads and processes population data for specific cities from ACS 5-year estimates, organized by age groups. The ACS data provides 5-year age groupings (0-4, 5-9, etc.) which can be disaggregated into single-year ages or aggregated into custom age groups.

Usage

getCityData(
  city_name,
  csv_path,
  age_groups = c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85),
  verbose = FALSE
)

Arguments

Value

A list containing:

Examples

# Load Hildale data with default 5-year age groups
hildale_data <- getCityData(
  city_name = "Hildale city, Utah",
  csv_path = system.file("extdata", "hildale_ut_2023.csv", package = "multigroup.vaccine")
)

# Load with single-year ages (disaggregated)
hildale_single <- getCityData(
  city_name = "Hildale city, Utah",
  csv_path = system.file("extdata", "hildale_ut_2023.csv", package = "multigroup.vaccine"),
  age_groups = NULL
)

# Load with custom age groups
hildale_custom <- getCityData(
  city_name = "Hildale city, Utah",
  csv_path = system.file("extdata", "hildale_ut_2023.csv", package = "multigroup.vaccine"),
  age_groups = c(0, 18, 65)
)

Calculate final size of outbreak: the total number of infections in each group, by solving the analytic final size equation

Description

Calculate final size of outbreak: the total number of infections in each group, by solving the analytic final size equation

Usage

getFinalSizeAnalytic(transmrates, recoveryrate, popsize, initR, initI, initV)

Arguments

Value

vector of final sizes (number of infected over whole outbreak) for each group

Examples

getFinalSizeAnalytic(transmrates = matrix(0.2, 2 ,2), recoveryrate = 0.3,
popsize = c(100, 150), initR = c(0, 0), initI = c(0, 1), initV = c(10, 10))

Estimate the distribution of final outbreak sizes by group using stochastic simulations of multi-group model

Description

Estimate the distribution of final outbreak sizes by group using stochastic simulations of multi-group model

Usage

getFinalSizeDist(n, transmrates, recoveryrate, popsize, initR, initI, initV)

Arguments

Value

a matrix with the final number infected from each group (column) in each simulation (row)

Examples

getFinalSizeDist(n = 10, transmrates = matrix(0.2, 2 ,2), recoveryrate = 0.3,
popsize = c(100, 150), initR = c(0, 0), initI = c(0, 1), initV = c(10, 10))

Estimate the distribution of final outbreak sizes by group using a hybrid model: stochastic simulations for smaller-sized outbreaks and deterministic ordinary differential equation model for "escaped" outbreaks

Description

Estimate the distribution of final outbreak sizes by group using a hybrid model: stochastic simulations for smaller-sized outbreaks and deterministic ordinary differential equation model for "escaped" outbreaks

Usage

getFinalSizeDistEscape(
  n,
  transmrates,
  recoveryrate,
  popsize,
  initR,
  initI,
  initV
)

Arguments

Value

a matrix with the final number infected from each group (column) in each simulation (row)

Examples

getFinalSizeDistEscape(n = 10, transmrates = matrix(0.2, 2 ,2), recoveryrate = 0.3,
popsize = c(100, 150), initR = c(0, 0), initI = c(0, 1), initV = c(10, 10))

Calculate outbreak final size, the total number of infections in each group, by numerically solving the multi-group ordinary differential equation

Description

Calculate outbreak final size, the total number of infections in each group, by numerically solving the multi-group ordinary differential equation

Usage

getFinalSizeODE(transmrates, recoveryrate, popsize, initR, initI, initV)

Arguments

Value

vector of final sizes (number of infected over whole outbreak) for each group

Examples

getFinalSizeODE(transmrates = matrix(0.2, 2 ,2), recoveryrate = 0.3,
popsize = c(100, 150), initR = c(0, 0), initI = c(0, 1), initV = c(10, 10))

Calculate outbreak size at a given time

Description

Calculate outbreak size at a given time

Usage

getSizeAtTime(time, transmrates, recoveryrate, popsize, initR, initI, initV)

Arguments

Value

a list with totalSize (total cumulative infections) and activeSize (total currently infected) in each group at the specified time

Examples

getSizeAtTime(time = 30, transmrates = matrix(0.2, 2 ,2), recoveryrate = 0.3,
popsize = c(100, 150), initR = c(0, 0), initI = c(0, 1), initV = c(10, 10))

Get state FIPS code by state name

Description

Get state FIPS code by state name

Usage

getStateFIPS(state_name)

Arguments

Value

Two-digit FIPS code as character

Examples

getStateFIPS("Utah")  # Returns "49"

List available counties for a state

Description

List available counties for a state

Usage

listCounties(
  state_fips,
  year = 2024,
  csv_path = NULL,
  cache_dir = NULL,
  verbose = FALSE
)

Arguments

Value

Character vector of county names

Examples

# Use the included example data
utah_counties <- listCounties(
  state_fips = "49", 
  year = 2024,
  csv_path = getCensusDataPath()
)


# Download from web (requires internet)
utah_counties_web <- listCounties(state_fips = "49", year = 2024)

# With caching
utah_counties_cached <- listCounties(state_fips = "49", cache_dir = "~/census_cache")

Ordinary differential equation function for multi-group susceptible-infectious-removed (SIR) model used as "func" argument passed to the ode() function from deSolve package

Description

Ordinary differential equation function for multi-group susceptible-infectious-removed (SIR) model used as "func" argument passed to the ode() function from deSolve package

Usage

odeSIR(time, state, par)

Arguments

Value

a list of three vectors of derivatives dS, dI, and dR for each group, evaluated at the given state values

Examples

# Intended only for use as the func argument to the ode() function from the deSolve package:
y0 <- c(S1 = 79999, S2 = 20000, I1 = 1, I2 = 0, R1 = 0, R2 = 0)
parms <- c(beta11 = 1.6e-6, beta21 = 1.5e-6, beta12 = 1.4e-6, beta22 = 8.7e-6, recoveryrate = 1/7)
times <- seq(0, 350, len = 10)
deSolve::ode(y0, times, odeSIR, parms)

Process census data with sex disaggregation

Description

Process census data with sex disaggregation

Usage

processCensusDataBySex(county_data, age_groups, verbose = FALSE)

Process census data without sex disaggregation

Description

Process census data without sex disaggregation

Usage

processCensusDataTotal(county_data, age_groups, verbose = FALSE)

Runs the shiny app

Description

Runs the shiny app

Usage

run_my_app(...)

Arguments

Details

The app featured in this package is the one presented in the shiny demo: https://shiny.posit.co/r/getstarted/shiny-basics/lesson1/index.html.

Value

Starts the execution of the app, printing the port on the console.

Examples

# To be executed interactively only
if (interactive()) {
  run_my_app()
}

Calculate transmission rate matrix for multi-group model with specified R0

Description

Calculate transmission rate matrix for multi-group model with specified R0

Usage

transmissionRates(R0, meaninf, reltransm)

Arguments

Value

a matrix of transmission rates to (row) and from (column) each group, in same time units as meaninf

Examples

transmissionRates(R0 = 15, meaninf = 7,
  reltransm = rbind(c(1, 0.5, 0.9), c(0.3, 1.9, 1), c(0.3, 0.6, 2.8)))

Calculate reproduction number for a multigroup model with a given state of vaccination and immunity

Description

Calculate reproduction number for a multigroup model with a given state of vaccination and immunity

Usage

vaxrepnum(meaninf, popsize, trmat, initR, initV, vaxeff)

Arguments

Value

the reproduction number

Examples

meaninf <- 7
popsize <- c(200, 800)
initR <- c(0, 0)
initV <- c(0, 0)
vaxeff <- 1
trmat <- matrix(c(0.63, 0.31, 0.19, 1.2), 2, 2)
vaxrepnum(meaninf, popsize, trmat, initR, initV, vaxeff)
vaxrepnum(meaninf, popsize, trmat, initR, initV = c(160, 750), vaxeff)