A worked example with Alabama Pre-K sites

Steps 1 and 2: downloading and routing

The first call needs a Census API key, and both need an internet connection:

# Needs a Census API key and an internet connection.
acs <- cacs_acs_prefetch(state = "AL", year = 2023)

# Needs the osrm package and an internet connection; the public OSRM server
# needs no key.
iso <- cacs_isochrone(
  sites       = cacs_alabama_sites,
  drive_times = c(5, 10, 15),
  provider    = "osrm",
  osrm_mode   = "demo"
)

cacs_acs_prefetch() returns an sf table with one row for each tract and variable and the columns GEOID, NAME, variable, estimate, moe, and geometry. cacs_isochrone() returns an sf table with one row for each site and drive time. The two example files have the same columns, so the other three steps run on them as they would on downloaded data:

# The three example files described above
sites <- readRDS(system.file(
  "extdata", "legacy_2025_sites.rds", package = "catchmentACS"
))
acs <- readRDS(system.file(
  "extdata", "sample_alabama_subset.rds", package = "catchmentACS"
))
iso <- readRDS(system.file(
  "extdata", "legacy_2025_isochrones.rds", package = "catchmentACS"
))

# Keep the three sites and the drive times of 5, 10, and 15 minutes.
analysis_sites <- sites |>
  filter(site_id %in% site_ids)
iso <- iso |>
  filter(site_id %in% site_ids, drive_time_min %in% c(5L, 10L, 15L))

list(
  acs_dim         = dim(acs),
  acs_crs         = sf::st_crs(acs)$epsg,
  acs_variables   = length(unique(acs$variable)),
  acs_tracts      = length(unique(acs$GEOID)),
  iso_dim         = dim(iso),
  iso_crs         = sf::st_crs(iso)$epsg
)

#> $acs_dim
#> [1] 12754     6
#> 
#> $acs_crs
#> [1] 4269
#> 
#> $acs_variables
#> [1] 14
#> 
#> $acs_tracts
#> [1] 911
#> 
#> $iso_dim
#> [1]  9 16
#> 
#> $iso_crs
#> [1] 4326

The ACS data hold 14 variables for 911 squares in NAD83 (EPSG:4269), and the drive-time areas are in WGS 84 (EPSG:4326): the two coordinate reference systems that cacs_intersect_weight() requires.

Step 3: overlapping tracts and their weights

cacs_intersect_weight() finds the tracts, here the squares, that overlap each drive-time area and combines their estimates:

weighted <- cacs_intersect_weight(
  iso_sf        = iso,
  acs_sf        = acs,
  weight_method = "area",
  verbose       = FALSE
)

dim(weighted)

#> [1] 126  23

The result has a row for each site, drive time, and ACS variable. These are the rows for the 10-minute area of AL_SITE_11:

area_11 <- weighted |>
  filter(site_id == "AL_SITE_11", drive_time_min == 10) |>
  select(variable, estimate, weight_basis, weight_sum, n_tracts)

area_11

#> # A tibble: 14 × 5
#>    variable   estimate weight_basis weight_sum n_tracts
#>    <chr>         <dbl> <chr>             <dbl>    <int>
#>  1 B01003_001    5156. coverage           1.11        3
#>  2 B11001_001    1786. coverage           1.11        3
#>  3 B17001_001    3637. coverage           1.11        3
#>  4 B17001_002     579. coverage           1.11        3
#>  5 B19013_001   54801. area_mean          1.11        3
#>  6 B19056_001    1693. coverage           1.11        3
#>  7 B19056_002     143. coverage           1.11        3
#>  8 B19301_001   36264. area_mean          1.11        3
#>  9 B22003_001    1731. coverage           1.11        3
#> 10 B22003_002     403. coverage           1.11        3
#> 11 B23025_001    2515. coverage           1.11        3
#> 12 B23025_002    1338. coverage           1.11        3
#> 13 B23025_003    2188. coverage           1.11        3
#> 14 B23025_005     122. coverage           1.11        3

A count, such as the total population (B01003_001), is the sum of the tract estimates, each multiplied by the tract’s coverage weight: the share of the tract’s area that lies inside the drive-time area (weight_basis = "coverage"). This assumes that whatever a variable counts is spread evenly over each tract. Median household income (B19013_001) and per capita income (B19301_001) are instead averages of the tract values, weighted by each tract’s share of the overlapping area ("area_mean"). These weights ignore how many people live in each tract, so the average can be far from the median or per-person value of the drive-time area when the tracts differ in population density. The help page of cacs_intersect_weight() defines both weights, and vignette("theory-spatial-aggregation", package = "catchmentACS") explains their assumptions.

n_tracts counts the tracts in the area, and weight_sum adds up their coverage weights. The two are equal only when every tract lies wholly inside the area. Here n_tracts is 3 and weight_sum is 1.11: the square around the site lies inside the 10-minute area, and two neighboring squares lie mostly outside it.

Step 4: margins of error

cacs_intersect_weight() has already combined the margins of error of the tracts. cacs_propagate_moe() computes them again at the confidence level given by its argument level, 0.9 by default, and records the formula used for each row in moe_formula_effective:

with_moe <- cacs_propagate_moe(weighted, verbose = FALSE)

with_moe |>
  select(site_id, drive_time_min, variable, estimate, moe,
         moe_formula_effective) |>
  arrange(site_id, drive_time_min, variable) |>
  head(8)

#> # A tibble: 8 × 6
#>   site_id    drive_time_min variable   estimate     moe moe_formula_effective
#>   <chr>               <int> <chr>         <dbl>   <dbl> <chr>                
#> 1 AL_SITE_07              5 B01003_001    4543.   612.  weighted_sum         
#> 2 AL_SITE_07              5 B11001_001    1506.   309.  weighted_sum         
#> 3 AL_SITE_07              5 B17001_001    1867.   110.  weighted_sum         
#> 4 AL_SITE_07              5 B17001_002    1154.   181.  weighted_sum         
#> 5 AL_SITE_07              5 B19013_001   68369  16626   weighted_mean        
#> 6 AL_SITE_07              5 B19056_001     998.    76.0 weighted_sum         
#> 7 AL_SITE_07              5 B19056_002     167.    35.0 weighted_sum         
#> 8 AL_SITE_07              5 B19301_001   16355   2422   weighted_mean

At the default level, the estimates and margins of error are those of step 3. weighted_sum, used for counts, is the formula for the margin of error of a sum, applied to the tract margins of error multiplied by the coverage weights. weighted_mean, used for medians and per-person values, is the same formula with the area shares. Both treat the tract estimates as independent and the weights as fixed. vignette("theory-moe-propagation", package = "catchmentACS") gives the formulas and discusses these assumptions.

Step 5: rates

cacs_derive_rates() adds five rows for each site and drive time, one for each rate in cacs_acs_default_rates, after the rows of the ACS variables. Each rate divides one weighted count by another, and its margin of error is computed from those of the two counts:

final <- cacs_derive_rates(with_moe, verbose = FALSE)

dim(final)

#> [1] 171  27

The rate rows are those whose estimand_family, the kind of quantity in the row, is "derived_rate":

final |>
  filter(estimand_family == "derived_rate") |>
  select(site_id, drive_time_min, variable, estimate, moe,
         moe_formula_effective) |>
  arrange(site_id, drive_time_min, variable) |>
  head(10)

#> # A tibble: 10 × 6
#>    site_id    drive_time_min variable     estimate     moe moe_formula_effective
#>    <chr>               <int> <chr>           <dbl>   <dbl> <chr>                
#>  1 AL_SITE_07              5 labor_force…   0.851  0.115   general_ratio_conser…
#>  2 AL_SITE_07              5 poverty_rate   0.618  0.104   general_ratio_conser…
#>  3 AL_SITE_07              5 snap_rate      0.0399 0.00553 general_ratio_conser…
#>  4 AL_SITE_07              5 ssi_rate       0.167  0.0373  general_ratio_conser…
#>  5 AL_SITE_07              5 unemp_rate     0.0620 0.0151  general_ratio_conser…
#>  6 AL_SITE_07             10 labor_force…   0.846  0.112   general_ratio_conser…
#>  7 AL_SITE_07             10 poverty_rate   0.602  0.0984  general_ratio_conser…
#>  8 AL_SITE_07             10 snap_rate      0.0426 0.00557 general_ratio_conser…
#>  9 AL_SITE_07             10 ssi_rate       0.165  0.0360  general_ratio_conser…
#> 10 AL_SITE_07             10 unemp_rate     0.0630 0.0149  general_ratio_conser…

By default, the margins of error of all five rates come from the ratio formula, recorded as general_ratio_conservative. The Census Bureau’s handbook gives this formula for a ratio whose numerator is not part of its denominator (U.S. Census Bureau 2020, chap. 8). The numerator of each rate is part of its denominator, so each rate is a proportion, for which the handbook gives the proportion formula; that formula never gives a wider margin of error than the ratio formula. The argument formula_dispatch of cacs_derive_rates() and cacs_run() chooses between the two: with "auto" or "proportion_subset", poverty_rate and labor_force_participation use the proportion formula, and snap_rate, ssi_rate, and unemp_rate keep the ratio formula. A row for which the value under the square root of the proportion formula is negative gets the ratio formula instead, with moe_fallback = TRUE. vignette("theory-derived-rates", package = "catchmentACS") describes the five rates, gives both formulas, and explains why unemp_rate keeps the ratio formula.

A plain table with the main columns

tibble::as_tibble() turns result into a plain tibble, which keeps its rows in their order. The tables use seven of its columns:

result_tbl <- tibble::as_tibble(result)

report_cols <- c(
  "site_id", "drive_time_min", "variable",
  "estimate", "moe", "n_tracts", "failure_origin"
)

result_tbl |>
  select(all_of(report_cols)) |>
  head(12)

#> # A tibble: 12 × 7
#>    site_id    drive_time_min variable   estimate     moe n_tracts failure_origin
#>    <chr>               <int> <chr>         <dbl>   <dbl>    <int> <chr>         
#>  1 AL_SITE_07              5 B01003_001   4543.   6.12e2        1 none          
#>  2 AL_SITE_07              5 B11001_001   1506.   3.09e2        1 none          
#>  3 AL_SITE_07              5 B17001_001   1867.   1.10e2        1 none          
#>  4 AL_SITE_07              5 B17001_002   1154.   1.81e2        1 none          
#>  5 AL_SITE_07              5 B19013_001  68369    1.66e4        1 none          
#>  6 AL_SITE_07              5 B19056_001    998.   7.60e1        1 none          
#>  7 AL_SITE_07              5 B19056_002    167.   3.50e1        1 none          
#>  8 AL_SITE_07              5 B19301_001  16355    2.42e3        1 none          
#>  9 AL_SITE_07              5 B22003_001   1704.   1.58e2        1 none          
#> 10 AL_SITE_07              5 B22003_002     68.0  7.00e0        1 none          
#> 11 AL_SITE_07              5 B23025_001   2020.   1.14e2        1 none          
#> 12 AL_SITE_07              5 B23025_002   1720.   2.12e2        1 none

In this result failure_origin, which records the step at which a row failed, is "none" on every row; ?cacs_run lists its other values and explains how missing tract estimates show up in the result.

Rates by site and drive time

tidyr::pivot_wider() turns the rate rows into a table with one row for each site and drive time and one column for each rate:

rate_matrix <- result_tbl |>
  filter(variable %in% names(cacs_acs_default_rates)) |>
  select(site_id, drive_time_min, variable, estimate) |>
  pivot_wider(names_from = variable, values_from = estimate) |>
  arrange(site_id, drive_time_min)

print(rate_matrix, width = Inf)

#> # A tibble: 9 × 7
#>   site_id    drive_time_min poverty_rate snap_rate ssi_rate unemp_rate
#>   <chr>               <int>        <dbl>     <dbl>    <dbl>      <dbl>
#> 1 AL_SITE_07              5        0.618    0.0399   0.167      0.0620
#> 2 AL_SITE_07             10        0.602    0.0426   0.165      0.0630
#> 3 AL_SITE_07             15        0.360    0.111    0.124      0.0922
#> 4 AL_SITE_08              5        0.457    0.334    0.0626     0.0331
#> 5 AL_SITE_08             10        0.428    0.315    0.0636     0.0378
#> 6 AL_SITE_08             15        0.182    0.182    0.0772     0.117 
#> 7 AL_SITE_11              5        0.161    0.244    0.0882     0.0499
#> 8 AL_SITE_11             10        0.159    0.233    0.0846     0.0559
#> 9 AL_SITE_11             15        0.149    0.155    0.0649     0.0984
#>   labor_force_participation
#>                       <dbl>
#> 1                     0.851
#> 2                     0.846
#> 3                     0.726
#> 4                     0.614
#> 5                     0.614
#> 6                     0.614
#> 7                     0.527
#> 8                     0.532
#> 9                     0.561

summary(result) holds the same rates in its element rates_per_site. Its element rates_per_site_moe writes each rate with its margin of error, rounded to three decimals:

print(summary(result)$rates_per_site_moe, width = Inf)

#> # A tibble: 9 × 7
#>   site_id    drive_time_min poverty_rate  snap_rate     ssi_rate     
#>   <chr>               <int> <chr>         <chr>         <chr>        
#> 1 AL_SITE_07              5 0.618 ± 0.104 0.040 ± 0.006 0.167 ± 0.037
#> 2 AL_SITE_07             10 0.602 ± 0.098 0.043 ± 0.006 0.165 ± 0.036
#> 3 AL_SITE_07             15 0.360 ± 0.047 0.111 ± 0.012 0.124 ± 0.017
#> 4 AL_SITE_08              5 0.457 ± 0.142 0.334 ± 0.077 0.063 ± 0.017
#> 5 AL_SITE_08             10 0.428 ± 0.125 0.315 ± 0.069 0.064 ± 0.016
#> 6 AL_SITE_08             15 0.182 ± 0.030 0.182 ± 0.026 0.077 ± 0.011
#> 7 AL_SITE_11              5 0.161 ± 0.045 0.244 ± 0.048 0.088 ± 0.021
#> 8 AL_SITE_11             10 0.159 ± 0.040 0.233 ± 0.044 0.085 ± 0.018
#> 9 AL_SITE_11             15 0.149 ± 0.020 0.155 ± 0.021 0.065 ± 0.010
#>   unemp_rate    labor_force_participation
#>   <chr>         <chr>                    
#> 1 0.062 ± 0.015 0.851 ± 0.115            
#> 2 0.063 ± 0.015 0.846 ± 0.112            
#> 3 0.092 ± 0.012 0.726 ± 0.096            
#> 4 0.033 ± 0.005 0.614 ± 0.054            
#> 5 0.038 ± 0.006 0.614 ± 0.051            
#> 6 0.117 ± 0.012 0.614 ± 0.068            
#> 7 0.050 ± 0.012 0.527 ± 0.090            
#> 8 0.056 ± 0.012 0.532 ± 0.082            
#> 9 0.098 ± 0.015 0.561 ± 0.076

A rate that could not be computed for a site and drive time is NA in these tables.

Sites ranked by a rate

The code below orders the three sites by the poverty rate of their 15-minute areas, highest first. This rate is the share of the people for whom poverty status is determined who are below the poverty level. The filter on failure_origin leaves out a site whose poverty rate is NA because routing failed or because a count or margin of error that the rate needs is missing:

ranked_15 <- result_tbl |>
  filter(
    drive_time_min == 15,
    variable == "poverty_rate",
    failure_origin == "none"
  ) |>
  arrange(desc(estimate)) |>
  transmute(
    rank         = row_number(),
    site_id,
    poverty_rate = estimate,
    moe_90       = moe
  )

ranked_15

#> # A tibble: 3 × 4
#>    rank site_id    poverty_rate moe_90
#>   <int> <chr>             <dbl>  <dbl>
#> 1     1 AL_SITE_07        0.360 0.0472
#> 2     2 AL_SITE_08        0.182 0.0296
#> 3     3 AL_SITE_11        0.149 0.0201

Two estimates $\hat{p}_j$ and $\hat{p}_k$ differ at the 90 percent level when $|\hat{p}_j - \hat{p}_k| > 1.645 \sqrt{\mathrm{SE}_j^2 + \mathrm{SE}_k^2}$, where the standard error $\mathrm{SE}$ of each estimate is its margin of error at the 90 percent level divided by 1.645. This is the Census Bureau’s test for comparing two estimates (U.S. Census Bureau 2020, chap. 7). The code below applies it to each site and the next one in the ranking:

# Each site against the next one in the ranking, at the 90 percent level
ranking_test <- ranked_15 |>
  mutate(
    se         = moe_90 / 1.645,
    next_site  = lead(site_id),
    difference = poverty_rate - lead(poverty_rate),
    threshold  = 1.645 * sqrt(se^2 + lead(se)^2),
    differ     = difference > threshold
  ) |>
  filter(!is.na(next_site)) |>
  select(site_id, next_site, difference, threshold, differ)

ranking_test

#> # A tibble: 2 × 5
#>   site_id    next_site  difference threshold differ
#>   <chr>      <chr>           <dbl>     <dbl> <lgl> 
#> 1 AL_SITE_07 AL_SITE_08     0.178     0.0557 TRUE  
#> 2 AL_SITE_08 AL_SITE_11     0.0328    0.0358 FALSE

The poverty rate of AL_SITE_07 exceeds that of AL_SITE_08 by 17.8 percentage points, more than the threshold of 5.6 points. The difference between AL_SITE_08 and AL_SITE_11, 3.3 points, is below its threshold of 3.6 points, so the test does not show that their poverty rates differ.

The margins of error in this test come from the ratio formula, which cacs_run() uses for all five rates by default. With formula_dispatch = "auto", poverty_rate uses the proportion formula where it can, and that formula never gives a wider margin of error. With the default margins of error, the test therefore finds no more differences than with those from "auto", and it can find fewer. The code below computes the rates again with "auto" and repeats the test:

# suppressWarnings() hides a warning that counts the rows for which the ratio
# formula replaced the proportion formula.
rates_auto <- suppressWarnings(
  cacs_derive_rates(with_moe, formula_dispatch = "auto", verbose = FALSE)
)

ranking_auto <- rates_auto |>
  filter(
    drive_time_min == 15,
    variable == "poverty_rate",
    failure_origin == "none"
  ) |>
  arrange(desc(estimate)) |>
  mutate(
    se         = moe / 1.645,
    next_site  = lead(site_id),
    difference = estimate - lead(estimate),
    threshold  = 1.645 * sqrt(se^2 + lead(se)^2),
    differ     = difference > threshold
  ) |>
  select(site_id, moe, moe_fallback, next_site, difference, threshold, differ)

ranking_auto

#> # A tibble: 3 × 7
#>   site_id       moe moe_fallback next_site  difference threshold differ
#>   <chr>       <dbl> <lgl>        <chr>           <dbl>     <dbl> <lgl> 
#> 1 AL_SITE_07 0.0472 TRUE         AL_SITE_08     0.178     0.0483 TRUE  
#> 2 AL_SITE_08 0.0102 FALSE        AL_SITE_11     0.0328    0.0225 TRUE  
#> 3 AL_SITE_11 0.0201 TRUE         <NA>          NA        NA      NA

With "auto", the margin of error of the 15-minute poverty rate of AL_SITE_08 is 0.0102 instead of 0.0296. For AL_SITE_07 and AL_SITE_11, the value under the square root of the proportion formula is negative, so their margins of error still come from the ratio formula (moe_fallback = TRUE). The difference between AL_SITE_08 and AL_SITE_11, 3.3 points, is now above its threshold of 2.3 points. Whether the test shows that these two poverty rates differ therefore depends on the formula for their margins of error.

The test treats the two estimates as independent, which is reasonable when the two areas share no tract, as here, where the sites are far apart. Areas that overlap share tracts, and so do the areas of one site for different drive times; the Limitations section of vignette("methodology", package = "catchmentACS") explains what this means for the test. With many sites, the test is run on many pairs, so some of the differences it shows may be due to chance.

One site at three drive times

The five rates of AL_SITE_08 for its three drive-time areas, with their margins of error and 90 percent confidence intervals:

focal_site <- "AL_SITE_08"

result_tbl |>
  filter(
    site_id == focal_site,
    variable %in% names(cacs_acs_default_rates)
  ) |>
  arrange(variable, drive_time_min) |>
  transmute(
    variable,
    drive_time_min,
    estimate,
    moe_90  = moe,
    ci_low  = estimate - moe,
    ci_high = estimate + moe
  )

#> # A tibble: 15 × 6
#>    variable                  drive_time_min estimate  moe_90 ci_low ci_high
#>    <chr>                              <int>    <dbl>   <dbl>  <dbl>   <dbl>
#>  1 labor_force_participation              5   0.614  0.0536  0.560   0.667 
#>  2 labor_force_participation             10   0.614  0.0514  0.562   0.665 
#>  3 labor_force_participation             15   0.614  0.0682  0.546   0.682 
#>  4 poverty_rate                           5   0.457  0.142   0.315   0.600 
#>  5 poverty_rate                          10   0.428  0.125   0.302   0.553 
#>  6 poverty_rate                          15   0.182  0.0296  0.153   0.212 
#>  7 snap_rate                              5   0.334  0.0769  0.257   0.411 
#>  8 snap_rate                             10   0.315  0.0685  0.247   0.384 
#>  9 snap_rate                             15   0.182  0.0263  0.156   0.208 
#> 10 ssi_rate                               5   0.0626 0.0166  0.0460  0.0792
#> 11 ssi_rate                              10   0.0636 0.0161  0.0475  0.0796
#> 12 ssi_rate                              15   0.0772 0.0113  0.0659  0.0885
#> 13 unemp_rate                             5   0.0331 0.00525 0.0279  0.0384
#> 14 unemp_rate                            10   0.0378 0.00563 0.0322  0.0434
#> 15 unemp_rate                            15   0.117  0.0115  0.105   0.128

ci_low and ci_high are the ends of the 90 percent confidence interval, the estimate minus and plus its margin of error. Because each of the three areas lies inside the next larger one, the test above does not apply to the differences between them as it stands.

The labor force participation rate of AL_SITE_08 happens to be almost the same at all three drive times. The made-up values of each variable were drawn separately, so in one square of the 15-minute area the labor force (B23025_002) exceeds the population 16 years and over (B23025_001), of which it is part.

A table to save

The last table keeps, for each rate, the estimate and its margin of error as proportions and as percentages. Its columns failure_origin, moe_formula_effective, and moe_fallback record whether and how the rate and its margin of error were computed:

export_tbl <- result_tbl |>
  filter(variable %in% names(cacs_acs_default_rates)) |>
  transmute(
    site_id,
    drive_time_min,
    variable,
    estimate,
    moe_90                = moe,
    estimate_pct          = 100 * estimate,
    moe_90_pct            = 100 * moe,
    failure_origin,
    moe_formula_effective,
    moe_fallback
  ) |>
  arrange(site_id, drive_time_min, variable)

head(export_tbl, 15)

#> # A tibble: 15 × 10
#>    site_id    drive_time_min variable   estimate  moe_90 estimate_pct moe_90_pct
#>    <chr>               <int> <chr>         <dbl>   <dbl>        <dbl>      <dbl>
#>  1 AL_SITE_07              5 labor_for…   0.851  0.115          85.1      11.5  
#>  2 AL_SITE_07              5 poverty_r…   0.618  0.104          61.8      10.4  
#>  3 AL_SITE_07              5 snap_rate    0.0399 0.00553         3.99      0.553
#>  4 AL_SITE_07              5 ssi_rate     0.167  0.0373         16.7       3.73 
#>  5 AL_SITE_07              5 unemp_rate   0.0620 0.0151          6.20      1.51 
#>  6 AL_SITE_07             10 labor_for…   0.846  0.112          84.6      11.2  
#>  7 AL_SITE_07             10 poverty_r…   0.602  0.0984         60.2       9.84 
#>  8 AL_SITE_07             10 snap_rate    0.0426 0.00557         4.26      0.557
#>  9 AL_SITE_07             10 ssi_rate     0.165  0.0360         16.5       3.60 
#> 10 AL_SITE_07             10 unemp_rate   0.0630 0.0149          6.30      1.49 
#> 11 AL_SITE_07             15 labor_for…   0.726  0.0958         72.6       9.58 
#> 12 AL_SITE_07             15 poverty_r…   0.360  0.0472         36.0       4.72 
#> 13 AL_SITE_07             15 snap_rate    0.111  0.0120         11.1       1.20 
#> 14 AL_SITE_07             15 ssi_rate     0.124  0.0169         12.4       1.69 
#> 15 AL_SITE_07             15 unemp_rate   0.0922 0.0118          9.22      1.18 
#> # ℹ 3 more variables: failure_origin <chr>, moe_formula_effective <chr>,
#> #   moe_fallback <lgl>

write.csv(export_tbl, "rates.csv", row.names = FALSE) saves it as a CSV file.