USPopulationSampler

R-CMD-check

USPopulationSampler is an R package designed to generate geospatial locations within census block groups (BG; the smallest geographic unit for which population counts are available) within target counties, states, or across the entirety of the U.S randomly selected according to population counts using the Census Bureau reference data. The package optionally allows users to embed temporal assignments to sampled locations using a provided Covid-19 reference data which allows deterministic assignment of dates to sampled geospatial locations.

Installation

You can install the development version of USPopulationSampler from GitHub with:

# install.packages("pak")
pak::pak("Techavoan/USPopulationSampler")

Loading the package

library(USPopulationSampler)

Loading necessary libraries

library(tidyverse)
#> ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
#> ✔ dplyr     1.2.0     ✔ readr     2.1.6
#> ✔ forcats   1.0.1     ✔ stringr   1.6.0
#> ✔ ggplot2   4.0.2     ✔ tibble    3.3.1
#> ✔ lubridate 1.9.4     ✔ tidyr     1.3.2
#> ✔ purrr     1.2.1     
#> ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
#> ✖ dplyr::filter() masks stats::filter()
#> ✖ dplyr::lag()    masks stats::lag()
#> ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(arrow)
#> 
#> Attaching package: 'arrow'
#> 
#> The following object is masked from 'package:lubridate':
#> 
#>     duration
#> 
#> The following object is masked from 'package:utils':
#> 
#>     timestamp

# For visualization of sampled points
library(usmap)
library(ggplot2)
library(sf)
#> Linking to GEOS 3.13.1, GDAL 3.11.4, PROJ 9.7.0; sf_use_s2() is TRUE
library(dplyr)
library(tigris)
#> To enable caching of data, set `options(tigris_use_cache = TRUE)`
#> in your R script or .Rprofile.
library(patchwork)
library(lubridate)

Loading the data from the Zenodo Repository

Important to the package is the Zenodo Repository: https://zenodo.org/records/21442375 which contains the reference datasets of the package. This includes the 2020 Decennial Census - Census Block Group Maps dataset which goes down to the block group level which is the smallest geographical unit that the US Census Bureau publishes population data. The other dataset is the Processed Covid-19 dataset which is a cleaned historical record of the Covid-19 pandemic in the United States at the county level that spans from 2020 to 2023.

Users can download both reference datasets using the pre-built functions load_bg_data() and load_covid_data() which are part of the package. More information regarding the datasets can be found in the Zenodo link attached.

# Loading the block-group level data
df_bg <- load_bg_data(overwrite = TRUE)
#> Downloading df_bg from Zenodo...
#> Downloading...   0% (0.00 / 74.78 MB)
#> Download complete.
covid_data <- load_covid_data(overwrite = TRUE)
#> Starting download from Zenodo...
#> Downloading...   0% (0.00 / 14.54 MB)
#> Download complete

Preparation of U.S 2020 Decennial Block Group Input Reference Data

When the synthetic spatiotemporal Covid-19 data replicates were generated, 241,891 block groups from our extracted 2020 Decennial Block group data frame were up for sampling. To acquire the 2020 Decennial Block group data frame containing the polygons and population counts, users must first acquire a census application interface (API) key from the United States Census Bureau https://api.census.gov/data/key_signup.html.

Users provide their organization and name to apply for the API key which is issued by email upon registration. This API key needs to be activated and copied to be used.

# Importing Necessary Libraries
library(tidyverse)
library(dplyr)
library(purrr)
library(tidycensus)

# Census API Key
# census_api_key("YOUR KEY", install = TRUE)

# Manually add Puerto Rico and District of Columbia to territories
states <- c(state.abb, "DC", "PR")

# Compiling block groups from all the states and territories
bg_list <- map(states, function(st) {
  get_decennial(
    geography = "block group",
    variables = "P1_001N",
    year = 2020,
    state = st,
    geometry = TRUE
  )
})

Upon acquiring the key, users need to create a state vector that contains the canonical 50 states, the District of Columbia, and the territory of Puerto Rico. The state.abb built-in character vector automatically contains the official two-letter postal abbreviations for all the 50 United States. The map() function provided through the purrr package iterates over all the specified regions and maps them to the get_decennial() function provided by the tidycensus package. get_decennial() sends a query to the U.S Census Bureau API and retrieves the census data for specified geography targeting block groups for the year 2020. P1_001N is the designation for the 2020 Decennial Census Redistricting Variable (Public Law 94-171) that represents the total population count for each block group as of April 1st, 2020. Setting geometry = TRUE, instructs tidycensus to return the TIGER/line geometry, the map-based polygon shapefile data provided by the census, for each block group in addition to their population counts. As a result, our output is a simple features (sf) object that contains geographic identifiers, population counts, and polygon boundaries for all the block groups retrieved.

# Combining all states and territories
bg_raw <- bind_rows(bg_list)

# Creating the sampling frame
df_bg <- bg_raw %>%
  transmute(
    GEO_ID = paste0("1500000US", GEOID),
    NAMELSAD = sub(",. *$", "", NAME),
    STATEFP = substr(GEOID, 1, 2),
    COUNTYFP = substr(GEOID, 3, 5),
    TRACTCE = substr(GEOID, 6, 11),
    BLKGRPCE = substr(GEOID, 12, 12),
    geometry,
    pop = value,
    FIPS = paste0(
      substr(GEO_ID, 1, 2),
      substr(GEO_ID, 3, 5)
    )
  ) %>%
  arrange(GEO_ID)

# Saving for package use (unecessary as data is provided via Zenodo)
save(df_bg, file = "df_bg.rda")

After combining all state-level Census downloads, the next step is to convert these raw Census fields into a cleaner structure that can be used for population-weighted sampling. Table 1 showcases a typical observation from the extracted raw block group data downloaded from the 2020 United States Census Bureau. Note that geometry is a MULTIPOLYGON data type which contains the border region of the specified block group within that county.

It is necessary to engineer new variable names to organize the extracted information. The Census API returns a 12-digit geographic identifier (GEOID) like 0100100201001 and our code prepends the Census geography summary prefix “1500000US” extrapolating GEOID to be 1500000US10010201001 which creates a unique block-group identifier consistent with Census conventions and TIGER/line shapefiles. A variable column called NAMELSAD through the substr() function strips trailing city descriptive information from the original NAME variable and outputs only “Block Group 1” which creates a concise block group label for that county. The state identifier (STATEFP) is extracted from the first 2 digits of GEOID. The county identifier (COUNTYFP) is extracted from the 3rd to 5th digit of GEOID. The census tract code (TRACTCE) is extracted from the 6th to 11 digit of GEOID. The block group code is the last digit of GEOID. The extracted geometry column is unchanged. The census variable “P1_001N” which contains the population counts is renamed to pop for downstream use and clarity. The FIPS codes for counties are constructed by combining STATEFP and COUNTYFP which is required as subsequent Covid-19 reference dataset is indexed based on county FIPS code. The dataset is sorted alphabetically based on GEOID.

A code reconstruction of the 2020 Census Block group sampling dataframe should include a total of 242,335 rows with each row representing a unique block group. 441 out of those has population counts of 0 adding no meaningful value to the sampler.

head(df_bg)
#>                  GEO_ID      NAMELSAD STATEFP COUNTYFP TRACTCE BLKGRPCE
#> 1 1500000US010010201001 Block Group 1      01      001  020100        1
#> 2 1500000US010010201002 Block Group 2      01      001  020100        2
#> 3 1500000US010010202001 Block Group 1      01      001  020200        1
#> 4 1500000US010010202002 Block Group 2      01      001  020200        2
#> 5 1500000US010010203001 Block Group 1      01      001  020300        1
#> 6 1500000US010010203002 Block Group 2      01      001  020300        2
#>                         geometry  pop  FIPS
#> 1 MULTIPOLYGON (((-86.5091 32...  575 01001
#> 2 MULTIPOLYGON (((-86.50407 3... 1200 01001
#> 3 MULTIPOLYGON (((-86.48093 3...  974 01001
#> 4 MULTIPOLYGON (((-86.47611 3... 1081 01001
#> 5 MULTIPOLYGON (((-86.47087 3... 2377 01001
#> 6 MULTIPOLYGON (((-86.46746 3...  839 01001

Preparation of Covid-19 Event-Level Reference Data

Raw annual county level Covid-19 data were tracked live by The New York Times on https://github.com/nytimes/covid-19-data starting from January 21, 2020 to March 23, 2023, on a day-to-day basis. The New York Times stopped tracking after this and relied on federal government reported data instead. The annual county level Covid-19 data aggregated case counts data stored as CSV files per year were combined and standardized to ensure that date formats are date class objects. FIPS codes in those files were also standardized to their typical 5-digit representation with leading zeros to ensure consistency between the 2020 decennial data and our Covid-19 reference data. Duplicate county-level observations were collapsed by retaining only the observation with the maximum cumulative count. Records with missing FIPS codes, null values, or invalid FIPS formats were removed.

library(dplyr)
library(tidyr)
library(stringr)
library(tidyverse)

# Reading County CSV Files Into R
county_2020 <- read_csv("C:/Users/techa/Downloads/us-counties-2020.csv")
county_2021 <- read_csv("C:/Users/techa/Downloads/us-counties-2021.csv")
county_2022 <- read_csv("C:/Users/techa/Downloads/us-counties-2022.csv")
county_2023 <- read_csv("C:/Users/techa/Downloads/us-counties-2023.csv")

# Defining the FIPS standardization function
fix_fips <- function(df) {
  df %>%
    mutate(
      date = as.Date(date),
      fips = if_else(
        is.na(fips),
        NA_character_,
        str_pad(as.character(fips), 5, pad = "0")
      )
    )
}

# Applying FIPS Standardization to the County CSV Files
county_2020 <- fix_fips(county_2020)
county_2021 <- fix_fips(county_2021)
county_2022 <- fix_fips(county_2022)
county_2023 <- fix_fips(county_2023)

# Combining All Counties Dataset Using Row Binding
covid_all <- bind_rows(
  county_2020,
  county_2021,
  county_2022,
  county_2023
)

# Grouping by FIP codes and date
covid_clean <- covid_all %>%
  group_by(fips, date) %>%
  # Collapsing Duplicate Counties
  summarise(
    county = first(county),
    state = first(state),
    cases = max(cases, na.rm = TRUE),
    deaths = max(deaths, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  filter(!is.na(fips)) %>% # Removing Invalid FIPS
  arrange(fips, date) # Ordering based on FIPS and Date

covid_daily <- covid_clean %>%
  group_by(fips) %>%
  arrange(date, .by_group = TRUE) %>%
  # Calculating Daily New Cases and Truncating Negative Revisions to zero
  mutate(
    new_cases = pmax(cases - lag(cases, default = 0), 0),
    new_deaths = pmax(deaths - lag(deaths, default = 0), 0)
  ) %>%
  ungroup()

Daily incident cases were computed by first differencing cumulative counts within each county, using zero as the baseline cumulative count for the first observation. Negative daily changes that arise from retrospective reporting, a term describing when health agencies revising historical records, were truncated to zero. Daily counts were then expanded into event-level observations using the uncount() built-it function of R. Figure 1 illustrates the complete pipeline for how the raw data was cleaned, prepared, and filtered before being used as a reference data set for the USPopulationSampler package. The filtered and cleaned dataset contained 100,539,437 total recorded Covid-19 cases from January 21, 2020 to March 23, 2023.

FIPS codes that were missing in the 2020 Decennial Census such as the US Virgin Islands which includes 78010 (St. Croix Island), 78020 (St. John Island), and 78030 (St. Thomas Island) were removed. Moreover, FIPS codes such as 69010, 69020, and 69030, which are regional designations from Guam, the U.S territory in Micronesia, were also removed. In total, these 6 FIPS codes that were removed accounted for 38,403 Covid-19 cases. The total cases is 100,501,034 accounting for the 50 canonical states, District of Columbia, and Puerto Rico.

Discontinued Alaskan FIPS designation 02261 for the Valdez-Cordova Area, 02997 and 02998 placeholder FIPS for missing Alaskan counties, and 48999 placeholder FIPS for missing Texas counties were kept because the USPopulationSampler package allows for state fallback sampling if FIPS codes provided can be tracked to those states and the removal of these FIPS codes could lead to a loss of thousands of cases.

# Calculating the sum to debug
sum(covid_daily$new_cases) # This should be 100539437

# Uncounting Combined Cases Into Event-Level Data
covid_data <- covid_daily %>%
  filter(new_cases > 0) %>%
  uncount(new_cases)

# Removing FIPS that are Missing From 2020 US Decennial Block Maps Data
setdiff(covid_data_fips, df_bg$FIPS)

invalid_fips <- c(
  "78010", "78020", "78030", # U.S. Virgin Islands
  "69100", "69110", "69120" # Guam
)

covid_data <- covid_data %>%
  filter(!fips %in% invalid_fips)

# Calculating the new sum to debug
nrow(covid_data) # This should be 100501034

A glimpse of the Covid-19 event-level reference data is shown below:

head(covid_data)
#> # A tibble: 6 × 4
#>   fips  date       county  state  
#>   <chr> <date>     <chr>   <chr>  
#> 1 01001 2020-03-24 Autauga Alabama
#> 2 01001 2020-03-25 Autauga Alabama
#> 3 01001 2020-03-25 Autauga Alabama
#> 4 01001 2020-03-25 Autauga Alabama
#> 5 01001 2020-03-26 Autauga Alabama
#> 6 01001 2020-03-26 Autauga Alabama

Explanation of US_pop_sampler parameters

The US_pop_sampler function is the main function of the USPopulationSampler package. Users can sample points from across America or by mentioning specific states of interest. The following parameters are apart of the function:

Usage of the US_pop_sampler function

The following sample code showcases the different features that users can expect of the US_pop_sampler function.

Sampling points from across the United States

If only N parameter is specified, the US_pop_sampler function will sample points from across the United States of America.

N can be varied to an integer of interest. Larger N values will take more time to sample points. Parallel is recommended for larger N values to speed up run time if user operating system permits. Typically N of 1 million or higher should induce the use of parallel processing as long as the operating system permits.

pts1 <- US_pop_sampler(
  N = 10000,
  df = df_bg
)
#> Requested total: 10,000
#> Final sampled points: 10,000
#> Active polygons: 9761
#> Elapsed time: 394.48 sec
options(tigris_use_cache = TRUE)

# State boundaries
us_states <- states(
  cb = TRUE,
  year = 2020
) |>
  st_transform(4326)

# Split state polygons
conus_states <- us_states |>
  filter(!STUSPS %in% c("AK", "HI", "PR"))

ak_state <- us_states |>
  filter(STUSPS == "AK")

hi_state <- us_states |>
  filter(STUSPS == "HI")

pr_state <- us_states |>
  filter(STUSPS == "PR")

# Split sampled points
pts_conus <- pts1 |>
  filter(
    lon >= -125,
    lon <= -66,
    lat >= 24,
    lat <= 50
  )

pts_ak <- pts1 |>
  filter(
    lon >= -180,
    lon <= -130,
    lat >= 50
  )

pts_hi <- pts1 |>
  filter(
    lon >= -161,
    lon <= -154,
    lat >= 18,
    lat <= 23
  )

pts_pr <- pts1 |>
  filter(
    lon >= -68.5,
    lon <= -65,
    lat >= 17.5,
    lat <= 18.7
  )

# ------------------------------------------------------------------
# Conus
# ------------------------------------------------------------------

p_conus <- ggplot() +
  geom_sf(
    data = conus_states,
    fill = "white",
    color = "grey50",
    linewidth = 0.2
  ) +
  geom_point(
    data = pts_conus,
    aes(lon, lat),
    alpha = 0.25,
    size = 0.15
  ) +
  coord_sf(
    xlim = c(-125, -66),
    ylim = c(24, 50),
    expand = FALSE
  ) +
  theme_minimal() +
  labs(
    title = "US_pop_sampler Output",
    subtitle = "10,000 synthetic locations sampled proportional to block-group population",
    x = "Longitude",
    y = "Latitude"
  ) +
  theme(
    plot.title.position = "plot",
    plot.margin = margin(
      t = 20,
      r = 20,
      b = 5,
      l = 5
    )
  )

# ------------------------------------------------------------------
# Alaska
# ------------------------------------------------------------------

p_ak <- ggplot() +
  geom_sf(
    data = ak_state,
    fill = "white",
    color = "grey50",
    linewidth = 0.2
  ) +
  geom_point(
    data = pts_ak,
    aes(lon, lat),
    alpha = 0.25,
    size = 0.15
  ) +
  coord_sf(
    xlim = c(-170, -130),
    ylim = c(54, 72),
    expand = FALSE
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    plot.title = element_text(size = 11),
    axis.text = element_blank(),
    axis.ticks = element_blank()
  )

# ------------------------------------------------------------------
# Hawaii
# ------------------------------------------------------------------

p_hi <- ggplot() +
  geom_sf(
    data = hi_state,
    fill = "white",
    color = "grey50",
    linewidth = 0.2
  ) +
  geom_point(
    data = pts_hi,
    aes(lon, lat),
    alpha = 0.25,
    size = 0.15
  ) +
  coord_sf(
    xlim = c(-161, -154),
    ylim = c(18, 23),
    expand = FALSE
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    plot.title = element_text(size = 11),
    axis.text = element_blank(),
    axis.ticks = element_blank()
  )

# ------------------------------------------------------------------
# Puerto Rico
# ------------------------------------------------------------------

p_pr <- ggplot() +
  geom_sf(
    data = pr_state,
    fill = "white",
    color = "grey50",
    linewidth = 0.2
  ) +
  geom_point(
    data = pts_pr,
    aes(lon, lat),
    alpha = 0.25,
    size = 0.15
  ) +
  coord_sf(
    xlim = c(-67.5, -65),
    ylim = c(17.8, 18.6),
    expand = FALSE
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank()
  )

# ------------------------------------------------------------------
# Combining all three figures
# ------------------------------------------------------------------

insets <- p_ak /
  p_hi /
  p_pr +
  plot_layout(
    heights = c(2, 2, 1) # Hawaii gets twice the height
  )

p_conus +
  insets +
  plot_layout(
    widths = c(4, 1.2)
  )

State Filtering Using state Parameter

To filter based on specific states, users can specify the state parameter using either full state names, abbreviated state names, or using STATEFP numbers. The state parameter accepts a vector of states so users can specify one or more states of interest.

Singular State

pts_cali <- US_pop_sampler(
  N = 10000,
  df = df_bg,
  state = "California"
)
#> Requested total: 10,000
#> Final sampled points: 10,000
#> Active polygons: 8013
#> Elapsed time: 338.50 sec
ca_counties <- suppressMessages(
  counties(
    state = "CA",
    cb = TRUE,
    year = 2020
  )
) |>
  st_transform(4326)
ggplot() +
  geom_sf(
    data = ca_counties,
    fill = "white",
    color = "grey70",
    linewidth = 0.15
  ) +
  geom_point(
    data = pts_cali,
    aes(lon, lat),
    alpha = 0.4,
    size = 0.25
  ) +
  coord_sf(expand = FALSE) +
  theme_minimal() +
  labs(
    title = "California Synthetic Population Sample",
    subtitle = "County boundaries shown",
    x = "Longitude",
    y = "Latitude"
  )

# How many sampled location points are in each county in California if we sample 10,000 points from California
table(substr(pts_cali$fips, 1, 5)) |>
  sort(decreasing = TRUE) |>
  head(10)
#> 
#> 06037 06073 06059 06065 06071 06085 06001 06067 06013 06019 
#>  2516   830   812   579   566   498   417   392   312   239

For multiple states, users can provide a vector containing the states of interest as shown below in the examples provided.

Using full state names

As long as the state name is spelled correctly, the function will recognize that state. Misspelled state names will not work.

pts2 <- US_pop_sampler(
  N = 10000,
  df = df_bg,
  state = c("California", "Oregon", "District of Columbia", "Puerto Rico", "North Dakota", "Wisconsin")
)
#> Requested total: 10,000
#> Final sampled points: 10,000
#> Active polygons: 8518
#> Elapsed time: 180.77 sec

We can visualize each of the stated regions for a better understanding of how the sampled points are located.

options(tigris_use_cache = TRUE)

# Assign State Labels
pts2$state <- dplyr::case_when(
  substr(pts2$fips, 1, 2) == "06" ~ "California",
  substr(pts2$fips, 1, 2) == "41" ~ "Oregon",
  substr(pts2$fips, 1, 2) == "11" ~ "District of Columbia",
  substr(pts2$fips, 1, 2) == "72" ~ "Puerto Rico",
  substr(pts2$fips, 1, 2) == "38" ~ "North Dakota",
  substr(pts2$fips, 1, 2) == "55" ~ "Wisconsin",
  TRUE ~ "Other"
)

# Download State Boundaries
states_sf <- tigris::states(
  cb = TRUE,
  year = 2020
) |>
  st_transform(4326) |>
  filter(
    NAME %in% c(
      "California",
      "Oregon",
      "District of Columbia",
      "Puerto Rico",
      "North Dakota",
      "Wisconsin"
    )
  )

# Keep Matching Names
states_sf$state <- states_sf$NAME

# Convert polygons to coordinates
states_df <- states_sf |>
  sf::st_cast("MULTILINESTRING") |>
  sf::st_coordinates() |>
  as.data.frame()

# Attach state names back
states_lookup <- states_sf |>
  st_drop_geometry() |>
  mutate(id = row_number()) |>
  select(id, state)

states_df$id <- states_df$L2

states_df <- left_join(
  states_df,
  states_lookup,
  by = "id"
)

# Plot
ggplot() +
  geom_path(
    data = states_df,
    aes(
      X,
      Y,
      group = interaction(id, L1)
    ),
    linewidth = 0.5
  ) +
  geom_point(
    data = pts2,
    aes(
      lon,
      lat
    ),
    alpha = 0.4,
    size = 0.25
  ) +
  facet_wrap(
    ~state,
    scales = "free"
  ) +
  theme_minimal() +
  labs(
    title = "Synthetic Population Samples",
    subtitle = "US_pop_sampler() output by state",
    x = "Longitude",
    y = "Latitude"
  )

The spelling of the full state names can also be case insensitive. Users can capitalize the names entirely or capitalize certain parts of the name. However, ensure that dictionary spelling of the state names is correct. Otherwise, the function will not recognize the state that is specified.

pts3 <- US_pop_sampler(
  N = 10000,
  df = df_bg,
  state = c("TeXAS", "VIRGINIA", "dELAwARE")
)
#> Requested total: 10,000
#> Final sampled points: 10,000
#> Active polygons: 7900
#> Elapsed time: 219.13 sec
options(tigris_use_cache = TRUE)

# Assign state labels from FIPS
pts3$state <- dplyr::case_when(
  substr(pts3$fips, 1, 2) == "48" ~ "Texas",
  substr(pts3$fips, 1, 2) == "51" ~ "Virginia",
  substr(pts3$fips, 1, 2) == "10" ~ "Delaware",
  TRUE ~ "Other"
)

# Download state boundaries
states_sf <- suppressMessages(
  tigris::states(
    cb = TRUE,
    year = 2020
  )
) |>
  st_transform(4326) |>
  filter(
    NAME %in% c(
      "Texas",
      "Virginia",
      "Delaware"
    )
  )

# Keep matching names
states_sf$state <- states_sf$NAME

# Convert polygons to coordinates
states_df <- states_sf |>
  sf::st_cast("MULTILINESTRING") |>
  sf::st_coordinates() |>
  as.data.frame()

# Attach state names back
states_lookup <- states_sf |>
  st_drop_geometry() |>
  mutate(id = row_number()) |>
  select(id, state)

states_df$id <- states_df$L2

states_df <- left_join(
  states_df,
  states_lookup,
  by = "id"
)

# Plot
ggplot() +
  geom_path(
    data = states_df,
    aes(
      X,
      Y,
      group = interaction(id, L1)
    ),
    linewidth = 0.5
  ) +
  geom_point(
    data = pts3,
    aes(
      lon,
      lat
    ),
    alpha = 0.4,
    size = 0.25
  ) +
  facet_wrap(
    ~state,
    scales = "free"
  ) +
  theme_minimal() +
  labs(
    title = "Synthetic Population Samples",
    subtitle = "US_pop_sampler() output by state",
    x = "Longitude",
    y = "Latitude"
  )

Using abbreviated state names

The function will recognize the current two-letter U.S state abbreviations that were developed by the U.S Postal Service in October 1963. These standardized codes were developed to facilitate automated mail processing with new ZIP codes. In the example, MN is Minnesota, NY is New York, and MI is Michigan.

pts4 <- US_pop_sampler(
  N = 10000,
  df = df_bg,
  state = c("MN", "NY", "MI")
)
#> Requested total: 10,000
#> Final sampled points: 10,000
#> Active polygons: 8247
#> Elapsed time: 265.21 sec
options(tigris_use_cache = TRUE)

# Assign state labels from FIPS
pts4$state <- dplyr::case_when(
  substr(pts4$fips, 1, 2) == "36" ~ "New York",
  substr(pts4$fips, 1, 2) == "27" ~ "Minnesota",
  substr(pts4$fips, 1, 2) == "26" ~ "Michigan",
  TRUE ~ "Other"
)

# Download state boundaries
states_sf <- suppressMessages(
  tigris::states(
    cb = TRUE,
    year = 2020
  )
) |>
  st_transform(4326) |>
  filter(
    NAME %in% c(
      "New York",
      "Minnesota",
      "Michigan"
    )
  )

# Keep matching names
states_sf$state <- states_sf$NAME

# Convert polygons to coordinates
states_df <- states_sf |>
  sf::st_cast("MULTILINESTRING") |>
  sf::st_coordinates() |>
  as.data.frame()

# Attach state names back
states_lookup <- states_sf |>
  st_drop_geometry() |>
  mutate(id = row_number()) |>
  select(id, state)

states_df$id <- states_df$L2

states_df <- left_join(
  states_df,
  states_lookup,
  by = "id"
)

# Plot
ggplot() +
  geom_path(
    data = states_df,
    aes(
      X,
      Y,
      group = interaction(id, L1)
    ),
    linewidth = 0.5
  ) +
  geom_point(
    data = pts4,
    aes(
      lon,
      lat
    ),
    alpha = 0.4,
    size = 0.25
  ) +
  facet_wrap(
    ~state,
    scales = "free"
  ) +
  theme_minimal() +
  labs(
    title = "Synthetic Population Samples",
    subtitle = "US_pop_sampler() output by state",
    x = "Longitude",
    y = "Latitude"
  )

Using STATEFP numbers

The function accepts STATEFP codes, which are two-digit numeric identifiers assigned by the U.S. Census Bureau. Leading zeros must be included (e.g., “06” for California). Only the 50 states, District of Columbia (11), and Puerto Rico (72) are supported.

A useful STATEFP reference table is provided below for users that are interested in using STATEFP numbers:

STATE_FIPS <- c(
  "AL" = "01", "AK" = "02", "AZ" = "04", "AR" = "05", "CA" = "06", "CO" = "08",
  "CT" = "09", "DE" = "10", "DC" = "11", "FL" = "12", "GA" = "13", "HI" = "15",
  "ID" = "16", "IL" = "17", "IN" = "18", "IA" = "19", "KS" = "20", "KY" = "21",
  "LA" = "22", "ME" = "23", "MD" = "24", "MA" = "25", "MI" = "26", "MN" = "27",
  "MS" = "28", "MO" = "29", "MT" = "30", "NE" = "31", "NV" = "32", "NH" = "33",
  "NJ" = "34", "NM" = "35", "NY" = "36", "NC" = "37", "ND" = "38", "OH" = "39",
  "OK" = "40", "OR" = "41", "PA" = "42", "RI" = "44", "SC" = "45", "SD" = "46",
  "TN" = "47", "TX" = "48", "UT" = "49", "VT" = "50", "VA" = "51", "WA" = "53",
  "WV" = "54", "WI" = "55", "WY" = "56", "PR" = "72"
)
pts5 <- US_pop_sampler(
  N = 10000,
  df = df_bg,
  state = c("12", "48", "35")
)
#> Requested total: 10,000
#> Final sampled points: 10,000
#> Active polygons: 8386
#> Elapsed time: 267.43 sec
options(tigris_use_cache = TRUE)

# Assign state labels from FIPS
pts5$state <- dplyr::case_when(
  substr(pts5$fips, 1, 2) == "48" ~ "Texas",
  substr(pts5$fips, 1, 2) == "12" ~ "Florida",
  substr(pts5$fips, 1, 2) == "35" ~ "New Mexico",
  TRUE ~ "Other"
)

# Download state boundaries
states_sf <- suppressMessages(
  tigris::states(
    cb = TRUE,
    year = 2020
  )
) |>
  st_transform(4326) |>
  filter(
    NAME %in% c(
      "Texas",
      "Florida",
      "New Mexico"
    )
  )

# Keep matching names
states_sf$state <- states_sf$NAME

# Convert polygons to coordinates
states_df <- states_sf |>
  sf::st_cast("MULTILINESTRING") |>
  sf::st_coordinates() |>
  as.data.frame()

# Attach state names back
states_lookup <- states_sf |>
  st_drop_geometry() |>
  mutate(id = row_number()) |>
  select(id, state)

states_df$id <- states_df$L2

states_df <- left_join(
  states_df,
  states_lookup,
  by = "id"
)

# Plot
ggplot() +
  geom_path(
    data = states_df,
    aes(
      X,
      Y,
      group = interaction(id, L1)
    ),
    linewidth = 0.5
  ) +
  geom_point(
    data = pts5,
    aes(
      lon,
      lat
    ),
    alpha = 0.4,
    size = 0.25
  ) +
  facet_wrap(
    ~state,
    scales = "free"
  ) +
  theme_minimal() +
  labs(
    title = "Synthetic Population Samples",
    subtitle = "US_pop_sampler() output by state",
    x = "Longitude",
    y = "Latitude"
  )

County Filtering Using fip_codes Parameter

Filtering based upon FIPS codes for counties is also a feature of the package. Users can supply FIPS codes in vector or scalar format.

pts6 <- US_pop_sampler(
  N = 10000,
  df = df_bg,
  fip_codes = c("06037", "06075", "06073")
)
#> Matched FIPS: 3
#> Missing FIPS: 0
#> Requested total: 10,000
#> Final sampled points: 10,000
#> Active polygons: 5859
#> Elapsed time: 203.58 sec

Using vector N to sample points

Users can pass a vector of N values through the N parameter. The matching is positional if the vector supplied is not named. Recall that a named vector is a programming data structure that functions like an ordered list of values where each value is paired with an explicit, readable string label (key) rather than just an index number. In this case, Washington gets 100 points sampled, Wisconsin gets 200 points sampled, and Kansas gets 500 points sampled. The total points sampled is 800 points. This vector N method works assuming the number of elements supplied in the N vector matches with the state vector. Similarly, users can supply STATEFP or case-insensitive names in the vector as well.

pts7 <- US_pop_sampler(
  N = c(100, 200, 500),
  df = df_bg,
  state = c("Washington", "Wisconsin", "Kansas")
)
#> Requested total: 800
#> Final sampled points: 800
#> Active polygons: 766
#> Elapsed time: 41.02 sec
options(tigris_use_cache = TRUE)

# Assign state labels from FIPS
pts7$state <- dplyr::case_when(
  substr(pts7$fips, 1, 2) == "55" ~ "Wisconsin",
  substr(pts7$fips, 1, 2) == "53" ~ "Washington",
  substr(pts7$fips, 1, 2) == "20" ~ "Kansas",
  TRUE ~ "Other"
)

# Download state boundaries
states_sf <- suppressMessages(
  tigris::states(
    cb = TRUE,
    year = 2020
  )
) |>
  st_transform(4326) |>
  filter(
    NAME %in% c(
      "Wisconsin",
      "Washington",
      "Kansas"
    )
  )

# Keep matching names
states_sf$state <- states_sf$NAME

# Convert polygons to coordinates
states_df <- states_sf |>
  sf::st_cast("MULTILINESTRING") |>
  sf::st_coordinates() |>
  as.data.frame()

# Attach state names back
states_lookup <- states_sf |>
  st_drop_geometry() |>
  mutate(id = row_number()) |>
  select(id, state)

states_df$id <- states_df$L2

states_df <- left_join(
  states_df,
  states_lookup,
  by = "id"
)

# Plot
ggplot() +
  geom_path(
    data = states_df,
    aes(
      X,
      Y,
      group = interaction(id, L1)
    ),
    linewidth = 0.5
  ) +
  geom_point(
    data = pts7,
    aes(
      lon,
      lat
    ),
    alpha = 0.4,
    size = 0.25
  ) +
  facet_wrap(
    ~state,
    scales = "free"
  ) +
  theme_minimal() +
  labs(
    title = "Synthetic Population Samples",
    subtitle = "US_pop_sampler() output by state",
    x = "Longitude",
    y = "Latitude"
  )

Similarly, users can also supply a vector of FIPS codes instead of states. users can supplied either states or FIPS codes but not both. The vector N works assuming the number of elements supplied in N matches the number of elements in fip_codes. This is positional in nature and the order in which each element appears indicate their matching similar to the state matching example shown above. The FIPS codes chosen are Los Angeles (06037), San Francisco (06075), and San Diego (06073) which are all counties in California.

pts8 <- US_pop_sampler(
  N = c(400, 300, 600),
  df = df_bg,
  fip_codes = c("06037", "06075", "06073")
)
#> Matched FIPS: 3
#> Missing FIPS: 0
#> Requested total: 1,300
#> Final sampled points: 1,300
#> Active polygons: 1044
#> Elapsed time: 52.10 sec

Users can plot just counties directly if interested:

options(tigris_use_cache = TRUE)

# Extract sampled county FIPS codes
county_fips <- unique(pts8$fips)

# Download county boundaries
counties_sf <- tigris::counties(
  cb = TRUE,
  year = 2020
) |>
  st_transform(4326) |>
  filter(GEOID %in% county_fips) |>
  mutate(
    county = paste0(NAME, ", ", STUSPS)
  )

# Convert county polygons to coordinates
counties_df <- counties_sf |>
  mutate(id = row_number()) |>
  st_cast("MULTILINESTRING") |>
  st_coordinates() |>
  as.data.frame()

# Attach county names
lookup <- counties_sf |>
  st_drop_geometry() |>
  mutate(id = row_number()) |>
  select(id, county)

counties_df$id <- counties_df$L2

counties_df <- left_join(
  counties_df,
  lookup,
  by = "id"
)

# Attach county names to sampled points
pts_plot <- pts8 |>
  left_join(
    counties_sf |>
      st_drop_geometry() |>
      select(GEOID, county),
    by = c("fips" = "GEOID")
  )

# Plot
ggplot() +
  geom_path(
    data = counties_df,
    aes(
      X,
      Y,
      group = interaction(id, L1)
    ),
    linewidth = 0.4
  ) +
  geom_point(
    data = pts_plot,
    aes(lon, lat),
    alpha = 0.4,
    size = 0.25
  ) +
  facet_wrap(
    ~county,
    scales = "free"
  ) +
  theme_minimal() +
  labs(
    title = "Synthetic Population Samples",
    subtitle = "US_pop_sampler() output by county",
    x = "Longitude",
    y = "Latitude"
  )

Usage of the US_pop_sampler function with parallel processing and streaming

The US_pop_sampler function uses sequential processing by default when parallel is set to FALSE. However, for larger N values or if users would like to speed up sampled points generation, users have the option of using parallel processing and the user can assign how many cores to use. By default, if parallel is set to TRUE, the function will use all available cores subtracted by 1.

Note that setting up cores can take longer if every core is used. Small number of cores if the operating system permits can achieve quicker generation time than mobilizing every core.

pts10 <- US_pop_sampler(
  N = 1e7,
  df = df_bg,
  parallel = TRUE,
  n_cores = 4,
)

Since the stream threshold is set at 1 million points by default, any N greater than this value will activate streaming. pts10 is stored as a list of 3 variables (path, batches, and total_points). Users can read into memory using the open_synthetic_dataset function to access the path of the sampled points.

# Finding the temporary path of the sampled points
ds <- open_synthetic_dataset(pts10)

Working with subsets

Arrow Datasets support lazy filtering and aggregation, allowing you to work with large datasets without loading everything into memory.

It is advised to inspect the schema of the dataset first before reading into memory.

names(ds)

dplyr::glimpse(ds)

You can extract records for a single county:

la_county <- ds |>
  filter(fips == "06037") |>
  collect()

Filtering just by state:

california <- ds |>
  filter(state == "California") |>
  collect()

Retrieve a random sample for exploratory analysis:

sample_pts <- ds |>
  slice_head(n = 100000) |>
  collect()

Select only the variables needed for analysis:

coords_only <- ds |>
  select(lon, lat, date) |>
  collect()

Summarize cases by county:

county_counts <- ds |>
  count(fips, sort = TRUE) |>
  collect()

Aggregate without materializing the full dataset:

county_counts <- ds |>
  count(fips, sort = TRUE) |>
  collect()

However, if the sampled points are not too memory intensive such as N is around 100 million or less then users can read into memory to access the whole generated synthetic points. 100 million points generally occupies around 4 GB of RAM on RStudio.

# Reading into memory if sampled points is not too large
pts_mem <- ds |>
  dplyr::collect()

Generating Synthetic Spatiotemporal Covid-19 Datasets

One of the main appeals of this package is the ability to generate synthetic spatiotemporal Covid-19 datasets. Users can load the processed Covid-19 dataset as follows:

# Loading the Covid-19 case data
covid_data <- load_covid_data(overwrite = TRUE)
#> Starting download from Zenodo...
#> Downloading...   0% (0.00 / 14.54 MB)
#> Download complete

Users can prepare the Covid-19 data by finding out how many FIPS codes exist, how many cases per fipcode, and store them as a named vector. This is particularly useful for Covid-19 data because there are some fipcodes in the Covid-19 data such as 48999 (Texas missing placeholder fipcode due to reporting errors), and discontinued Alaskan counties (02261, 02997, and 02998) that have not been dropped.

Named vector ensures positional matching with each fipcode being paired with the exact number of Covid-19 case that appears in them. For the discontinued fipcodes, the function ensures reconciliation by sampling them based on the other existing blockgroups populations within those states of Texas and Alaska respectively.

# Count observed Covid-19 records by county FIPS
fips_counts <- table(covid_data$fips)

# Convert to named vector for FIPS-based allocation
N_vec <- as.integer(fips_counts)
names(N_vec) <- names(fips_counts)

# Counties observed in the COVID dataset
covid_fips <- names(fips_counts)

Once the number of cases (fips_counts) per fips, the list of fipcodes in Covid-19 (covid_fips), and the matching of the fipcodes to the number of cases (N_vec) have been generated, users can then generate their own synthetic Covid-19 dataset.

# Generate synthetic spatial locations
pts_large <- US_pop_sampler(
  N = N_vec,
  df = df_bg,
  fip_codes = covid_fips,
  parallel = TRUE,
  n_cores = 4 # users can increase cores depending on their operating system
)

To ensure that the synthetic Covid-19 dataset has the temporal components, users can use the assign_dates_to_pts function.

Assigning dates to pts based on Covid-19 data from 2020 to 2023

Assigning dates is achieved from a deterministic point of view using Covid-19 data through the assign_dates_to_pts function. The Covid-19 data contains 100,501,034 recorded cases of Covid-19 across the United States of America from the year 2020 to 2023. Each row represents a recorded case of Covid-19 and contains the date of the case, the FIPS code of the county, and the county name.

The assign_dates_to_pts function has the following parameters: -pts A synthetic_points object returned by US_pop_sampler() or a character path to a directory containing parquet files.

-covid_data A data frame containing the FIPS code, date, and county name of recorded Covid-19 cases from 2020 to 2023

-output_dir An optional parameter that allows the user to specify the directory where the dated parquet files will be written to in disk.

-seed An optional parameter that allows the user to specify a number to ensure reproducible results in date assignment.

Users can then assign dates to pts_large using the assign_dates_to_pts function. The output will be written to disk as parquet files. Users can specify the output directory for the parquet files or if not specified, the files will be written to a temporary directory in temp files. In the example below, the synthetic spatiotemporal Covid-19 dataset is written into a folder called Synthetic Datasets and stored in a sub-folder called sampled_dataset new.

dated_pts_large <- assign_dates_to_pts(
  pts = pts_large,
  covid_data = covid_data,
  "C:/Users/techa/Desktop/Synthetic Datasets/sampled_dataset_new"
)

The returned object contains metadata about the generated dataset including the location of the Parquet files on disk.

Use open_synthetic_dataset() to lazily access the dataset without loading it entirely into memory. This returns an Arrow Dataset object that supports dplyr and arrow operations.

ds <- open_synthetic_dataset(dated_pts_large)

You can also open a dataset directly from its folder path.

ds <- open_synthetic_dataset("C:/Users/techa/Desktop/Synthetic Datasets/sampled_dataset_new")

Similarly, you can also load the full dataset into memory. Do note that collecting very large datasets can exhaust available RAM which is the primary reason the package uses streaming in the first place. However, with Covid-19 dataset, it is still safe to read entirely into memory. Users can still follow the same operations as included in earlier subsection “Working with subsets” to explore the generated Covid-19 dataset.

sampled_dataset <- ds |> collect()

If you prefer a single Parquet file instead of partitioned files, first collect the dataset into memory and then write it using Arrow.

sampled_dataset <- ds |> collect()

write_parquet(sampled_dataset, "C:/Users/techa/Desktop/Synthetic Datasets/pts_large_with_dates.parquet")

A single Parquet file can be loaded directly into memory using read_parquet().

sampled_dataset <- read_parquet(
  "C:/Users/techa/Desktop/Synthetic Datasets/pts_large_with_dates.parquet"
)

Accessing Zenodo Repository of Pre-generated Synthetic Covid-19 data

28 synthetic pre-generated synthetic Covid-19 datasets are published on the Zenodo repository that users can access and use for required analyses. The function load_synthetic_covid enables users to do the following. It is recommended to load each data one at a time to avoid memory issues on local machines.

# Loading synthetic Covid-19 dataset
sampled_data <- load_synthetic_covid(18) # users can specify a number between 1 to 28 (add overwrite = TRUE for first time users)
#> Synthetic dataset 18 already exists. Skipping download.

Due to memory constraints in the RStudio IDE, it is reccomended to read only a subset of the entire dataset. Users can choose to read 10% of the generated dataset as shown and randomly permute the order of the points for clear visibility of the years.

plot_data <- sampled_data |>
  group_by(fips, date) |>
  slice_sample(prop = 0.1) |>
  ungroup() |>
  arrange(runif(n())) # randomly permute row order

Users can create a column to isolate the year component of the synthetic spatiotemporal dataset.

plot_data <- plot_data |>
  mutate(
    year = factor(year(date))
  )

Users can then plot the synthetic spatiotemporal Covid-19 provided by the package.

options(tigris_use_cache = TRUE)

# State boundaries

us_states <- states(
  cb = TRUE,
  year = 2020
) |>
  st_transform(4326)

conus_states <- us_states |>
  filter(!STUSPS %in% c("AK", "HI", "PR"))

ak_state <- us_states |>
  filter(STUSPS == "AK")

hi_state <- us_states |>
  filter(STUSPS == "HI")

pr_state <- us_states |>
  filter(STUSPS == "PR")

# Splitting sampled points

pts_conus <- plot_data |>
  filter(
    lon >= -125,
    lon <= -66,
    lat >= 24,
    lat <= 50
  )

pts_ak <- plot_data |>
  filter(
    lon >= -180,
    lon <= -130,
    lat >= 50
  )

pts_hi <- plot_data |>
  filter(
    lon >= -161,
    lon <= -154,
    lat >= 18,
    lat <= 23
  )

pts_pr <- plot_data |>
  filter(
    lon >= -68.5,
    lon <= -65,
    lat >= 17.5,
    lat <= 18.7
  )

# ------------------------------------------------------------------
# Stylistic Choices
# ------------------------------------------------------------------

year_colors <- c(
  "2020" = "Red",
  "2021" = "#009E73", # Green
  "2022" = "Dark Blue", # Vermillion
  "2023" = "Yellow" # Purple
)

inset_theme <- theme_minimal() +
  theme(
    axis.title = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    plot.title = element_text(
      size = 9,
      face = "bold",
      hjust = 0.5
    ),
    legend.position = "none"
  )

# CONUS Territories

p_conus <- ggplot() +
  geom_sf(
    data = conus_states,
    fill = "white",
    color = "grey35",
    linewidth = 0.3
  ) +
  geom_point(
    data = pts_conus,
    aes(lon, lat, color = year),
    alpha = 0.2,
    size = 0.05
  ) +
  coord_sf(
    xlim = c(-125, -66),
    ylim = c(24, 50),
    expand = FALSE
  ) +
  scale_color_manual(
    values = year_colors,
    name = "Year"
  ) +
  theme_minimal() +
  labs(
    x = "Longitude",
    y = "Latitude"
  ) +
  theme(
    # Pulls the map up by using a negative top margin
    plot.margin = margin(t = -15, r = 20, b = 5, l = 5, unit = "pt")
  )

# ------------------------------------------------------------------
# Alaska
# ------------------------------------------------------------------

p_ak <- ggplot() +
  geom_sf(
    data = ak_state,
    fill = "white",
    color = "grey35",
    linewidth = 0.3
  ) +
  geom_point(
    data = pts_ak,
    aes(
      lon,
      lat,
      color = year
    ),
    alpha = 0.2,
    size = 0.02
  ) +
  coord_sf(
    xlim = c(-170, -130),
    ylim = c(54, 72),
    expand = FALSE
  ) +
  scale_color_manual(
    values = year_colors,
    name = "Year"
  ) +
  labs(title = "Alaska") +
  inset_theme

# ------------------------------------------------------------------
# Hawaii
# ------------------------------------------------------------------

p_hi <- ggplot() +
  geom_sf(
    data = hi_state,
    fill = "white",
    color = "grey35",
    linewidth = 0.3
  ) +
  geom_point(
    data = pts_hi,
    aes(
      lon,
      lat,
      color = year
    ),
    alpha = 0.2,
    size = 0.02
  ) +
  coord_sf(
    xlim = c(-161, -154),
    ylim = c(18, 23),
    expand = FALSE
  ) +
  scale_color_manual(
    values = year_colors,
    name = "Year"
  ) +
  labs(title = "Hawaii") +
  inset_theme


# ------------------------------------------------------------------
# Puerto Rico
# ------------------------------------------------------------------

p_pr <- ggplot() +
  geom_sf(
    data = pr_state,
    fill = "white",
    color = "grey35",
    linewidth = 0.3
  ) +
  geom_point(
    data = pts_pr,
    aes(
      lon,
      lat,
      color = year
    ),
    alpha = 0.2,
    size = 0.02
  ) +
  coord_sf(
    xlim = c(-67.5, -65),
    ylim = c(17.8, 18.6),
    expand = FALSE
  ) +
  scale_color_manual(
    values = year_colors,
    name = "Year"
  ) +
  labs(title = "Puerto Rico") +
  inset_theme

# ------------------------------------------------------------------
# Insets
# ------------------------------------------------------------------

insets <- p_ak /
  p_hi /
  p_pr +
  plot_layout(
    heights = c(2, 2, 1)
  )

# Plotting Figure

final_plot <- p_conus +
  insets +
  plot_layout(
    # This allocates more horizontal canvas space to the main CONUS map.
    widths = c(4.5, 1),
    guides = "collect"
  ) &
  theme(
    legend.position = "right",
    # Shrined legend title font size
    legend.title = element_text(
      size = 10,
      face = "bold"
    ),
    # Shrined legend item text font size
    legend.text = element_text(
      size = 9
    ),
    # Reduced the physical size of the colored legend squares
    legend.key.height = unit(0.4, "cm"),
    legend.key.width = unit(0.4, "cm"),
    # Tightens up the spacing around the keys
    legend.spacing.y = unit(0.1, "cm")
  ) &
  guides(
    color = guide_legend(
      override.aes = list(
        size = 3.5, # Tunable to make legend keys smaller
        alpha = 1
      )
    )
  )

# Render the final plot with annotations
final_plot +
  plot_annotation(
    title = "Synthetic COVID-19 Observations Across the United States",
    theme = theme(
      plot.title = element_text(
        size = 14,
        face = "bold",
        hjust = 0,
        margin = margin(
          t = 15,
          r = 0,
          b = 5,
          l = 0,
          unit = "pt"
        )
      )
    )
  )

Acknowledgements

The development and deployment of this package would not have been possible without the contributions of Dr. Andrew Holbrook, Associate Professor of Biostatistics, at the Fielding School of Public Health, at University of California, Los Angeles (UCLA). He has been an advisor who I could lean on for interpersonal and professional development. His insights on user interface, statistics, and package development have been instrumental in finalizing and improving the package. His willingness to help Graduate students like myself gain research experience is a testament to how Professors help uplift others and impart knowledge to aspiring students.

Moreover, I am grateful to my mentor, Jasen Zhang, PhD Candidate in Biostatistics at UCLA. Jasen’s patience in helping me to learn the theory behind the Hawkes and Spatiotemporal Hawkes processes were inspiring and motivational. I truly thank both of them for helping to actualize this package.

Thank you to my family, my parents, siblings, aunts, uncles, cousins, friends, and grandparents for their continued support and belief in my potential.

Finally, thank you to “R Packages” authored by Hadley Wickham and Jennifer Bryan for providing the schema for how to create packages in R. Without these individuals, many of the tools used in R would simply cease to exist.