Area weighting of tract estimates
Source:vignettes/theory-spatial-aggregation.Rmd
theory-spatial-aggregation.RmdA drive-time area is drawn from travel times, and its edge cuts
across census tracts: the area around a site takes in some tracts whole
and parts of others. catchmentACS estimates characteristics of the
people and households in the area, such as the number below the poverty
level, from the American Community Survey (ACS) estimates for these
tracts. The weight of each tract is computed from the area it shares
with the drive-time area. The margins of error of these estimates,
half-widths of confidence intervals at the 90 percent level by default,
are derived in
vignette("theory-moe-propagation", package = "catchmentACS"),
and the rates are described in
vignette("theory-derived-rates", package = "catchmentACS").
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.36The 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.000 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.000000The 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.000 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.225033The 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 coverageThe 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.