Area weighting of tract estimates

Area weighting and its assumption

As in vignette("methodology", package = "catchmentACS"), we write $I_s$ for the drive-time area of site $s$ at one drive time and $T_j$ for census tract $j$. The part of the tract inside the drive-time area is $I_s \cap T_j$, and $\lvert\,\cdot\,\rvert$ denotes area. An ACS estimate describes a whole tract, while the drive-time area may contain only part of it.

Area weighting is a simple and widely used way to move counts from one set of areas to another. It assumes that whatever a variable counts is spread evenly over the tract’s area, which rarely holds in practice (Comber and Zeng 2019, 8). Let $\rho_j$ be the number per unit area, in tract $T_j$, of whatever the variable counts: people, households, or a subgroup such as the people below the poverty level. If $\rho_j$ is the same everywhere in the tract, the number in any part $S$ of the tract is proportional to the area of that part:

$$\text{number in } S = \rho_j \, \lvert S \rvert = \big(\text{number in } T_j\big) \cdot \frac{\lvert S \rvert}{\lvert T_j \rvert}.$$

With $S = I_s \cap T_j$, the part of a tract count that lies in the drive-time area is the count multiplied by the share of the tract’s area inside the drive-time area. This share is the coverage weight of the next section.

Counts, such as the number of people below the poverty level or the number of households, can be split this way. Medians and per-person values, such as median household income (B19013_001) and per capita income (B19301_001), do not add up over area: half of a tract does not have half of its per capita income. For these, the package averages the tract values with weights proportional to the area each tract shares with the drive-time area. This average stands in for the median or per-person value of the drive-time area, but its weights follow area and ignore how many people live in each tract. It can therefore be far from that value when the overlapping tracts differ in population density, even if people are spread evenly within each tract; the next section gives an example. The margin of error reported with such an average describes the sampling error of the average, not the difference between it and the drive-time area’s own value, and that difference does not shrink as the ACS margins of error do. A rate, such as the poverty rate, is the ratio of two counts, and each of the two counts is split by area.

The assumption fails where the people of a tract live on a small part of its area. Suppose that 80 percent of a tract is uninhabited woodland and that its residents live on the remaining 20 percent. A drive-time area that covers all of the woodland and none of the settled part gives the tract a coverage weight of 0.8, so it is credited with 80 percent of the tract’s residents, although none of them live inside. The error grows with how unevenly people are spread within the tracts, as in a large rural tract that contains a small town. Methods that use other data, such as land cover, to place people within tracts relax the assumption (Comber and Zeng 2019, 3–5). catchmentACS uses area weighting only: weight_method = "population" is not implemented yet, and using it gives an error.

Coverage weights and area shares

Both weights are computed from the overlap area $\lvert I_s \cap T_j \rvert$ and differ in what it is divided by. The coverage weight of tract $j$ divides it by the area of the tract,

$$w^{cov}_{sj} = \frac{\lvert I_s \cap T_j \rvert}{\lvert T_j \rvert} \;\in\; [0, 1],$$

which is 1 when the tract lies wholly inside the drive-time area and 0 when the two do not overlap. The coverage weights of the overlapping tracts need not sum to one, since each is a share of a different tract. Their sum adds up the covered fractions of the tracts: two whole tracts and half of a third give 2.5, and so do five tracts that are each half covered. The column weight_sum of the result reports the sum for each count, median, and per-person value.

The area share of tract $j$ divides the overlap area by the total overlap area of the tracts,

$$w^{mean}_{sj} = \frac{\lvert I_s \cap T_j \rvert}{\sum_k \lvert I_s \cap T_k \rvert},$$

so the area shares sum to one. These sums, like all sums over tracts below, run over the tracts kept for the site, drive time, and variable. A tract is kept if its coverage weight is above the argument min_weight (see the section on the calculation) and the ACS data have a row for it and the variable. The area shares are proportional to the coverage weights multiplied by the tract areas, $w^{mean}_{sj} \propto w^{cov}_{sj}\, \lvert T_j \rvert$. When the tracts have the same area, the area shares are therefore the coverage weights divided by their sum.

A count is estimated by the sum of the tract estimates $Y_j$ weighted by coverage weights,

$$\widehat{Y}_s = \sum_j w^{cov}_{sj}\, Y_j ,$$

which grows with the part of each tract inside the drive-time area. With area shares in place of coverage weights, the result would be an average of the tract counts, on the scale of one tract.

A median or a per-person value is estimated by the average of the tract estimates $X_j$ weighted by area shares,

$$\widehat{\bar X}_s = \sum_j w^{mean}_{sj}\, X_j .$$

Because the area shares are nonnegative and sum to one, $\widehat{\bar X}_s$ is at least the smallest and at most the largest of the tract values. A sum of the same values weighted by coverage weights, whose sum can be above or below one, can fall outside that range.

The area shares weight each tract by its overlap area and ignore how many people live there. Two made-up tracts show how much this can matter. Tract 1 has an area of 1 km², 1,000 residents, and a per capita income of 10,000 dollars, and lies wholly inside the drive-time area. Tract 2 has an area of 100 km², 1,000 residents, and a per capita income of 50,000 dollars, and 10 km² of it lie inside the drive-time area. The code below applies the package’s formula for a per-person value to these numbers by hand. For comparison, it also computes the per capita income of the people inside, on the assumption that within each tract both the residents and their income per person are the same everywhere.

tracts <- data.frame(
  tract      = c(1, 2),
  area_km2   = c(1, 100),       # area of the tract
  inside_km2 = c(1, 10),        # area of its part inside the drive-time area
  residents  = c(1000, 1000),
  income_pc  = c(10000, 50000)  # per capita income, dollars
)
tracts$w_cov  <- tracts$inside_km2 / tracts$area_km2         # coverage weight
tracts$w_mean <- tracts$inside_km2 / sum(tracts$inside_km2)  # area share
tracts

#>   tract area_km2 inside_km2 residents income_pc w_cov     w_mean
#> 1     1        1          1      1000     10000   1.0 0.09090909
#> 2     2      100         10      1000     50000   0.1 0.90909091

# The package's formula: the average of the tract values weighted by area shares
area_share_average <- sum(tracts$w_mean * tracts$income_pc)

# The people inside and their per capita income (both assumed even within each tract)
residents_inside <- sum(tracts$w_cov * tracts$residents)
income_inside    <- sum(tracts$w_cov * tracts$residents * tracts$income_pc) /
  residents_inside

c(area_share_average = area_share_average,
  residents_inside   = residents_inside,
  income_inside      = income_inside)

#> area_share_average   residents_inside      income_inside 
#>           46363.64            1100.00           13636.36

The package’s formula gives 46,364 dollars, while the 1,100 people inside have a per capita income of 13,636 dollars. Most of the overlap area lies in tract 2, but most of the people inside live in tract 1. The per capita income of the people inside is their total income divided by their number, a ratio of two counts. The area-share average equals it when the overlapping tracts have the same number of residents per unit area.

A rate, such as the poverty rate, is estimated by the ratio of two coverage-weighted sums of tract estimates, $\widehat{A}_s = \sum_j w^{cov}_{sj}\, A_j$ for its numerator and $\widehat{B}_s = \sum_j w^{cov}_{sj}\, B_j$ for its denominator. vignette("theory-derived-rates", package = "catchmentACS") shows when this ratio is an average of the tract rates, weighted by the part of each tract’s denominator inside the drive-time area, and describes the margins of error of rates.

How the package chooses the weights

Each row of the results records the kind of estimate in the column estimand_family and the weights used in weight_basis: coverage weights for counts and rates, and area shares for medians and per-person values. The kind is read from the ACS code alone, as the help page of cacs_intersect_weight() and vignette("methodology", package = "catchmentACS") describe. Median age (B01002_001), for instance, is a median from a table that the package does not list, so it is added up like a count, without a warning.

For a median, the area-share average is a looser stand-in than for a per-person value. Even when the tracts have the same population density, an average of tract medians is in general not the median of the combined population, which depends on the income distribution in each tract. The ACS publishes its 5-year detailed tables for all areas down to block groups, census tracts included, and one of these tables, B19001, gives the number of households in 16 income brackets. The package does not use these brackets: the area-share average of the tract medians stands in for the median of the drive-time area. The brackets themselves can be aggregated, because they are counts: asking for B19001_017 in variables gives the number of households in that bracket for the drive-time area, with a coverage-weighted sum like any other count. What the package has no function for is turning brackets into a median.

Measuring the areas

cacs_intersect_weight() computes the overlaps and their areas on a flat map. Before measuring any area, it projects the drive-time areas and the tracts to EPSG:5070 (NAD83 / Conus Albers) and measures all areas there, in square meters. This Albers projection is equal-area: the area of a shape on the projected map is its area on the reference ellipsoid. Every area in a run is measured on this one map.

A coverage weight is the ratio of two areas within the same tract, so it changes little with the way area is measured. EPSG defines this projection for the contiguous 48 states. The package checks its own box for that scope, drawn around the contiguous United States and the District of Columbia with a quarter of a degree to spare. The bounding boxes of both inputs must lie between 24.14 and 49.64 degrees north and between -125.25 and -66.68 degrees east. Data for Alaska, Hawaii, or Puerto Rico give an error.

The calculation for each site

Before it combines anything, cacs_intersect_weight() repairs invalid geometries in both inputs (see below) and skips, with a warning, tracts whose area is zero or not finite. The skipped tracts are listed in the skipped_geoids attribute of the result. Tracts with a coverage weight at or below min_weight (by default 1e-6), such as tracts that only touch the edge of the area, are dropped before the area shares and the estimates are computed.

The help page of cacs_intersect_weight() lists the steps of the calculation. It also describes the single row of NA values, with failure_origin = "intersection", that a site and drive time gets when no tract is left or the calculation fails.

Invalid geometries, such as polygons whose edges cross, are repaired with sf::st_make_valid() or, if some remain invalid, with a buffer of zero width (sf::st_buffer()), each time with a warning. This is done for both inputs and again for the overlaps of each pair. Repair can change a shape and its area. A geometry that is still invalid after both attempts is replaced by an empty one. An empty tract has zero area and is skipped, an empty drive-time area gives its pair the row described above, and an empty overlap adds nothing. If every geometry of an input is still invalid, the function stops with an error.

Recomputing the weights for one site

The weights and one count for one site can be recomputed from the two areas that cacs_intersect_weight() records for each tract when keep_tract_audit = TRUE. The data for this example ship with the package. In these data, squares of nearly equal size, scattered with gaps between them and carrying random ACS values, take the place of census tracts. Each site sits at the center of one of the squares, and its drive-time areas are circles around it. The ACS values are made up, so the numbers below show the arithmetic and not a real place. At a drive time of 10 minutes, the circle of site AL_SITE_17 contains the site’s own square and cuts two other squares at its edge. Because of the gaps, the squares cover only a small part of the circle; real census tracts tile a state, so for an area of this size the sum of the coverage weights there would be far larger.

# Data bundled with the package
iso <- readRDS(system.file(
  "extdata", "legacy_2025_isochrones.rds", package = "catchmentACS"
))
acs <- readRDS(system.file(
  "extdata", "sample_alabama_subset.rds", package = "catchmentACS"
))

site_id    <- "AL_SITE_17"
drive_time <- 10L

iso_one <- iso[
  iso$site_id == site_id & iso$drive_time_min == drive_time, ,
  drop = FALSE
]

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

The cacs_tract_audit attribute of the result lists, for each tract that the drive-time area overlaps, the area of the overlap (int_area_m2) and the area of the tract (tract_area_m2), in square meters on EPSG:5070:

areas <- attr(weighted, "cacs_tract_audit") |>
  select(GEOID, int_area_m2, tract_area_m2) |>
  arrange(desc(int_area_m2))

areas

#> # A tibble: 3 × 3
#>   GEOID       int_area_m2 tract_area_m2
#>   <chr>             <dbl>         <dbl>
#> 1 01125001001    9233961.      9233961.
#> 2 01127001001    1038977.      9233961.
#> 3 01123001001    1038971.      9233961.

The coverage weights divide the overlap area by the tract area, and the area shares divide it by the total overlap area:

hand <- areas |>
  mutate(
    w_cov  = int_area_m2 / tract_area_m2,        # |I_s n T_j| / |T_j|
    w_mean = int_area_m2 / sum(int_area_m2)      # |I_s n T_j| / sum_k |I_s n T_k|
  )

hand |> select(GEOID, w_cov, w_mean)

#> # A tibble: 3 × 3
#>   GEOID       w_cov w_mean
#>   <chr>       <dbl>  <dbl>
#> 1 01125001001 1.00  0.816 
#> 2 01127001001 0.113 0.0918
#> 3 01123001001 0.113 0.0918

c(
  sum_w_cov  = sum(hand$w_cov),    # not 1 in general: each has its own denominator
  sum_w_mean = sum(hand$w_mean)    # 1 by construction
)

#>  sum_w_cov sum_w_mean 
#>   1.225033   1.000000

The coverage weights are 1, 0.113, and 0.113: the site’s own square lies wholly inside the circle, and the other two only partly. Their sum, 1.225, is the overlap area counted in squares, because these three squares have the same area; across the whole data set the squares differ in area by about five percent.

The number of people below the poverty level (B17001_002) is a count, so the package estimates it by $\widehat{Y}_s = \sum_j w^{cov}_{sj}\, Y_j$. The code below takes the three tract estimates from the ACS data and forms this sum:

counts <- acs |>
  sf::st_drop_geometry() |>
  filter(GEOID %in% hand$GEOID, variable == "B17001_002") |>
  select(GEOID, Y = estimate)

hand_total <- hand |>
  select(GEOID, w_cov) |>
  left_join(counts, by = "GEOID")

hand_total

#> # A tibble: 3 × 3
#>   GEOID       w_cov     Y
#>   <chr>       <dbl> <dbl>
#> 1 01125001001 1.00   1020
#> 2 01127001001 0.113  1104
#> 3 01123001001 0.113   828

hand_Y_hat       <- sum(hand_total$w_cov * hand_total$Y)   # sum_j w_cov * Y_j
hand_weight_sum  <- sum(hand_total$w_cov)                  # sum_j w_cov

c(hand_Y_hat = hand_Y_hat, hand_weight_sum = hand_weight_sum)

#>      hand_Y_hat hand_weight_sum 
#>     1237.382230        1.225033

The row of the package’s result for the same variable has estimand_family = "spatial_total" and weight_basis = "coverage":

pkg_row <- weighted |>
  filter(variable == "B17001_002") |>
  select(variable, estimate, weight_sum, n_tracts,
         estimand_family, weight_basis)

pkg_row

#> # A tibble: 1 × 6
#>   variable   estimate weight_sum n_tracts estimand_family weight_basis
#>   <chr>         <dbl>      <dbl>    <int> <chr>           <chr>       
#> 1 B17001_002    1237.       1.23        3 spatial_total   coverage

The hand-computed sum, 1,237.4 people, is the estimate in the row above, and the sum of the coverage weights, 1.225, is its weight_sum. The same two areas are recorded for every tract of every site and drive time in the result, so any row can be recomputed this way.

References