jared
d95c9032c5
The Census of Governments is a complete census only in years ending in 2 and
7. Every other year is a sample, and the sample varies enormously. Neither
cog_geographic_rollup() nor cog_peer_compare()/cog_find_peers() had any
concept of "the universe": each summed or labelled whichever govids happened
to have rows and returned that with nothing distinguishing "every government
reported" from "a fifth of them did".
On the bundled fixture, Wisconsin's 608-city universe rolls up 597
governments in FY2012 and 112 in FY2019. The peer side is worse exposure, not
better: a Madison-scale cohort looks stable because Madison is large, while
governments matched to a small target sit in exactly the population band the
sample cycle hits hardest. Chilton's 15-peer cohort reports 15 of 15 in
FY2012 and 3 of 15 in FY2019.
Implements the owner's settled design: coverage = c("all", "census",
"consistent") on all three verbs, defaulting to "all" so nothing currently
calling them changes, PLUS always-on provenance$coverage carrying per-year
n_units_reporting / n_units_expected / is_census_year and
provenance$coverage_mode. cog_explain() prints a "Reporting coverage"
section. The default mode can no longer mislead silently, which is the point
-- using these verbs correctly must not require knowing the survey calendar.
Decisions worth stating:
- n_units_expected is the universe the CALLER named, not the national one.
That is what makes the ratio mean something: "597 of the 608 Wisconsin
cities you asked about". For peers it is the cohort size, counted over
peer rows only -- including the target would inflate every count by one
and make a cohort that has entirely stopped reporting look non-empty.
- The coverage table is built from the REQUESTED years, not the years
present in the result, so a year in which nothing reported still appears
with n_units_reporting = 0. A year that vanishes silently is precisely
the disclosure failure at issue.
- "census" filters years BEFORE the query, and aborts when the range holds
no census year rather than returning an empty result for a query the
caller believes they made.
- "consistent" exempts the peer-comparison target: it is the subject of the
comparison, not a member of the cohort being balanced, and dropping it
would leave nothing to compare. The summary_* quantiles are computed
AFTER the filter so they describe the cohort actually returned.
- is_census_year is documented as a statement about the survey CALENDAR,
never a claim of completeness -- FY1967 is a census year in which only 97
of Wisconsin's 608 cities report (DoD 3). n_units_reporting is the number
that tells the truth.
On cog_find_peers(), where there is no year range, coverage governs the
cohort VINTAGE: "census" snaps to the most recent census year with an
observed population, so a cohort is not built from a sample year in which
most of the candidate universe is absent. "consistent" is a comparison-time
concept and selects like "all" there, carried on the result for
cog_peer_compare().
One fix to the committed test, which was internally inconsistent. It pinned
n_units_reporting == 597 for FY2012 AND asserted that number equals a raw
cross-check that answers 595. Both numbers are right for different questions:
VERNON VILLAGE and WAUKESHA VILLAGE carry type = 3 in `long` (their
as-of-year identity, as townships) while the xwalk lists them as govs_type =
2 (their present identity, as villages) -- schema v6 made the long table's
geography present-harmonized but `type` still reads as-of-year. The rollup
counts against the requested govid set, so 597 answers "how many of the
governments I asked about reported". The cross-check now scopes to that same
universe instead of to long.type/long.fips_state; it still reads raw parquet
rather than going through the verb under test.
Suite: 670 pass / 0 fail / 2 skip (was 658/0/3). rcmdcheck clean.
The two remaining skips are #11 and #12.
366 lines
16 KiB
R
366 lines
16 KiB
R
# R/peers.R
|
||
|
||
#' Find peer governments by similarity criteria
|
||
#'
|
||
#' Selects peer governments by combinations of government type, state, and
|
||
#' population range at a chosen `year`. Peers are ordered by `|log(pop_ratio)|`
|
||
#' ascending (closest to the target's population first).
|
||
#'
|
||
#' @param target_govid Character scalar — `canonical_govid` of the target.
|
||
#' @param year Integer scalar. Cohort vintage. When `NULL` (default), uses the
|
||
#' most recent year for which the target has an observed population in
|
||
#' `gov_population_yearly`.
|
||
#' @param same_type If `TRUE` (default) restrict peers to the target's
|
||
#' `govs_type`.
|
||
#' @param same_state If `TRUE` restrict peers to the target's `fips_state`.
|
||
#' Default `FALSE`.
|
||
#' @param pop_range Length-2 numeric vector giving lower/upper bounds.
|
||
#' @param is_ratio If `TRUE` (default) `pop_range` is multiplied by the
|
||
#' target's population at `year` to produce absolute bounds. If `FALSE`,
|
||
#' `pop_range` is interpreted as absolute population counts.
|
||
#' @param max_peers Integer cap on the number of peers returned.
|
||
#' @param coverage Survey-cycle handling; see [cog_peer_compare()]. Here it
|
||
#' governs the cohort VINTAGE when `year` is `NULL`: `"census"` snaps to the
|
||
#' most recent census year with an observed population, so a cohort is not
|
||
#' built from a sample year in which most of the candidate universe is
|
||
#' absent. `"consistent"` needs a year range, which cohort selection does not
|
||
#' have, so it selects like `"all"` and is carried on the result as
|
||
#' `attr(x, "coverage")` for [cog_peer_compare()].
|
||
#' @return Tibble with columns `canonical_govid`, `gov_name`, `fips_state`,
|
||
#' `population`, `pop_ratio`, `rank`. The cohort year is attached as
|
||
#' `attr(x, "cohort_year")`.
|
||
#' @export
|
||
cog_find_peers <- function(target_govid,
|
||
year = NULL,
|
||
same_type = TRUE,
|
||
same_state = FALSE,
|
||
pop_range = c(0.7, 1.3),
|
||
is_ratio = TRUE,
|
||
max_peers = 10L,
|
||
coverage = c("all", "census", "consistent")) {
|
||
coverage <- .validate_coverage(coverage)
|
||
if (!is.character(target_govid) || length(target_govid) != 1L) {
|
||
cli::cli_abort("`target_govid` must be a length-1 character string.")
|
||
}
|
||
if (!is.numeric(pop_range) || length(pop_range) != 2L ||
|
||
pop_range[1] >= pop_range[2]) {
|
||
cli::cli_abort("`pop_range` must be a length-2 numeric with lo < hi.")
|
||
}
|
||
if (!is.null(year) &&
|
||
(!(is.numeric(year) || is.integer(year)) || length(year) != 1L)) {
|
||
cli::cli_abort("`year` must be NULL or a length-1 integer.")
|
||
}
|
||
|
||
con <- .ensure_session()
|
||
|
||
# Confirm target exists in the xwalk and pull govs_type / fips_state.
|
||
meta_sql <- sprintf(
|
||
"SELECT canonical_govid, gov_name, govs_type, fips_state
|
||
FROM canonical_fips_xwalk
|
||
WHERE canonical_govid = %s",
|
||
.sql_lit_chr(target_govid)
|
||
)
|
||
meta <- DBI::dbGetQuery(con, meta_sql)
|
||
if (nrow(meta) == 0L) {
|
||
cli::cli_abort(c(
|
||
"govid {target_govid} not found in corpus.",
|
||
i = "v0.1 covers types 0-3 only (state/county/city/township); see vignette('coverage-scope')."
|
||
))
|
||
}
|
||
|
||
cohort_year <- .resolve_cohort_year(con, target_govid, year, coverage)
|
||
|
||
pop_sql <- sprintf(
|
||
"SELECT population FROM gov_population_yearly
|
||
WHERE canonical_govid = %s AND year = %d",
|
||
.sql_lit_chr(target_govid), as.integer(cohort_year)
|
||
)
|
||
target_pop <- DBI::dbGetQuery(con, pop_sql)$population
|
||
if (length(target_pop) == 0L || is.na(target_pop) || target_pop <= 0) {
|
||
cli::cli_abort(c(
|
||
"Target {target_govid} has no observed population in {cohort_year}.",
|
||
i = "Use a year for which population is observed; see gov_population_yearly."
|
||
))
|
||
}
|
||
|
||
if (isTRUE(is_ratio)) {
|
||
lo <- target_pop * pop_range[1]
|
||
hi <- target_pop * pop_range[2]
|
||
} else {
|
||
lo <- pop_range[1]; hi <- pop_range[2]
|
||
}
|
||
|
||
preds <- c(
|
||
sprintf("p.canonical_govid != %s", .sql_lit_chr(target_govid)),
|
||
sprintf("p.year = %d", as.integer(cohort_year)),
|
||
sprintf("p.population BETWEEN %.6f AND %.6f", lo, hi)
|
||
)
|
||
if (isTRUE(same_type)) preds <- c(preds, sprintf("x.govs_type = %d", meta$govs_type))
|
||
if (isTRUE(same_state)) preds <- c(preds, sprintf("x.fips_state = %s", .sql_lit_chr(meta$fips_state)))
|
||
|
||
peers_sql <- sprintf(
|
||
"SELECT p.canonical_govid, x.gov_name, x.fips_state, p.population,
|
||
p.population / %.6f AS pop_ratio
|
||
FROM gov_population_yearly p
|
||
JOIN canonical_fips_xwalk x USING (canonical_govid)
|
||
WHERE %s
|
||
ORDER BY ABS(LN(CAST(p.population AS DOUBLE) / %.6f))
|
||
LIMIT %d",
|
||
target_pop,
|
||
paste(preds, collapse = " AND "),
|
||
target_pop,
|
||
as.integer(max_peers)
|
||
)
|
||
peers <- tibble::as_tibble(DBI::dbGetQuery(con, peers_sql))
|
||
peers$rank <- if (nrow(peers) > 0L) seq_len(nrow(peers)) else integer(0)
|
||
attr(peers, "cohort_year") <- as.integer(cohort_year)
|
||
attr(peers, "pop_range") <- as.numeric(pop_range)
|
||
attr(peers, "is_ratio") <- isTRUE(is_ratio)
|
||
attr(peers, "coverage") <- coverage
|
||
attr(peers, "is_census_year") <- .is_census_year(cohort_year)
|
||
peers
|
||
}
|
||
|
||
# `coverage` picks the cohort vintage when the caller did not name one.
|
||
# "census" snaps to the most recent CENSUS year with an observed population,
|
||
# so a cohort is not silently built from a sample year in which most of the
|
||
# candidate universe is absent. "consistent" is a comparison-time concept --
|
||
# it needs a year RANGE, which cohort selection does not have -- so it selects
|
||
# like "all" here and is carried on the result for cog_peer_compare().
|
||
#' @noRd
|
||
.resolve_cohort_year <- function(con, target_govid, year,
|
||
coverage = "all") {
|
||
if (!is.null(year)) return(as.integer(year))
|
||
if (identical(coverage, "census")) {
|
||
sql <- sprintf(
|
||
"SELECT MAX(year) AS y FROM gov_population_yearly
|
||
WHERE canonical_govid = %s AND year %% 10 IN (2, 7)",
|
||
.sql_lit_chr(target_govid)
|
||
)
|
||
y <- DBI::dbGetQuery(con, sql)$y
|
||
if (length(y) > 0L && !is.na(y)) return(as.integer(y))
|
||
cli::cli_abort(c(
|
||
"{.code coverage = \"census\"} found no census year with an observed population for {target_govid}.",
|
||
i = "Pass an explicit {.arg year}, or use {.code coverage = \"all\"}."
|
||
), class = "uscogdata_no_census_years")
|
||
}
|
||
sql <- sprintf(
|
||
"SELECT MAX(year) AS y FROM gov_population_yearly
|
||
WHERE canonical_govid = %s",
|
||
.sql_lit_chr(target_govid)
|
||
)
|
||
y <- DBI::dbGetQuery(con, sql)$y
|
||
if (length(y) == 0L || is.na(y)) {
|
||
cli::cli_abort(
|
||
"Target {target_govid} has no observed population in any year."
|
||
)
|
||
}
|
||
as.integer(y)
|
||
}
|
||
|
||
#' Compare a target government against a peer set
|
||
#'
|
||
#' Pulls spending for the target plus a peer set (either a
|
||
#' [cog_find_peers()] result or a character vector of `canonical_govid`) and
|
||
#' appends peer-distribution summary rows (`summary_p25`, `summary_p50`,
|
||
#' `summary_p75`) so the result can be faceted by `role` in a single ggplot
|
||
#' call. Those summary rows are quantiles **within each category**, not
|
||
#' quantiles of each peer's total — see the `@return` section before summing
|
||
#' them.
|
||
#'
|
||
#' @param target_govid Character scalar.
|
||
#' @param peers A tibble from [cog_find_peers()] or a character vector of
|
||
#' `canonical_govid`s.
|
||
#' @param category Character scalar or vector.
|
||
#' @param years Integer vector.
|
||
#' @param per_capita Default `TRUE` — peer compare usually normalizes by
|
||
#' population.
|
||
#' @param adjust_to_year Integer base year for CPI-U conversion or `NULL`.
|
||
#' @param expenditure_concept `"direct"` (default) or `"total"`. Currently only
|
||
#' `"direct"` is accepted; the `"total"` option exists in [cog_spending()] for
|
||
#' single-government queries but cannot be used here because combining Total
|
||
#' across peer sets counts intergovernmental transfers twice.
|
||
#' @param coverage How to handle the Census of Governments survey cycle,
|
||
#' which is a **complete census only in years ending in 2 and 7** -- every
|
||
#' other year is a sample, and the sample varies enormously (on the bundled
|
||
#' fixture, Wisconsin's 608-city universe reports 597 governments in FY2012
|
||
#' and 112 in FY2019).
|
||
#'
|
||
#' * `"all"` (default) -- every unit that reported that year. Unchanged
|
||
#' behaviour, so existing code keeps working.
|
||
#' * `"census"` -- census years only. Aborts if the requested range holds
|
||
#' none, rather than silently returning nothing.
|
||
#' * `"consistent"` -- only units reporting in *every* requested year, giving
|
||
#' a balanced panel.
|
||
#'
|
||
#' Regardless of mode, `provenance$coverage` always carries per-year
|
||
#' `n_units_reporting`, `n_units_expected` and `is_census_year`, and
|
||
#' `provenance$coverage_mode` records the mode. `is_census_year` is a
|
||
#' statement about the **survey calendar**, never a claim of completeness:
|
||
#' FY1967 is a census year in which only 97 of Wisconsin's 608 cities
|
||
#' report. `n_units_reporting` is the number that tells the truth.
|
||
#'
|
||
#' The comparison target is exempt from `"consistent"` balancing -- it is the
|
||
#' subject of the comparison, not a member of the cohort -- and the
|
||
#' `summary_*` quantiles are computed AFTER the filter, so they describe the
|
||
#' cohort actually returned. `n_units_reporting` counts peers only, against
|
||
#' the cohort size: "3 of your 15 peers reported in FY2019".
|
||
#' @return Tibble matching [cog_spending()]'s columns, plus a `role`
|
||
#' column taking values `"target"`, `"peer"`, `"summary_p25"`,
|
||
#' `"summary_p50"`, or `"summary_p75"`, `target_rank` (target's rank
|
||
#' among target+peers at `max(years)`, NA for other rows), and
|
||
#' `cohort_year` (the year used to build the peer cohort, read from
|
||
#' `attr(peers, "cohort_year")`; `NA` when `peers` was a bare character
|
||
#' vector). Provenance reports `verb = "cog_peer_compare"`, `peer_count`,
|
||
#' `cohort_year`, and `cohort_govids`.
|
||
#'
|
||
#' **The `summary_*` rows are per-category quantiles: they are not additive.**
|
||
#' Each one is computed **within each `(year, spend_subtype,
|
||
#' category)` cell** across the peer set, so a `summary_p50` row is *the
|
||
#' median peer's value in that one category*, not *the value of the median
|
||
#' peer's total*. The median peer for Police and the median peer for Fire
|
||
#' are usually different governments, so summing `summary_*` rows across
|
||
#' categories does not give any peer's total and misstates the band it
|
||
#' appears to describe — measured at −32.7% to +251.0% across 24 years on
|
||
#' one cohort, with a sign flip at FY2012.
|
||
#'
|
||
#' Facet by `role` **and** `category` (the documented use, and what the
|
||
#' rows are built for). For a genuine "median peer's total spending" line,
|
||
#' sum each peer's own categories first and take the quantile of those
|
||
#' per-government totals:
|
||
#'
|
||
#' ```r
|
||
#' library(dplyr)
|
||
#' cmp |>
|
||
#' filter(role %in% c("target", "peer")) |>
|
||
#' group_by(year, role, canonical_govid) |>
|
||
#' summarise(total = sum(amt_per_capita_real, na.rm = TRUE), .groups = "drop") |>
|
||
#' filter(role == "peer") |>
|
||
#' group_by(year) |>
|
||
#' summarise(p50 = quantile(total, 0.5, na.rm = TRUE))
|
||
#' ```
|
||
#' @export
|
||
cog_peer_compare <- function(target_govid, peers, category, years,
|
||
per_capita = TRUE, adjust_to_year = NULL,
|
||
expenditure_concept = c("direct", "total"),
|
||
coverage = c("all", "census", "consistent")) {
|
||
call <- match.call()
|
||
expenditure_concept <- match.arg(expenditure_concept)
|
||
coverage <- .validate_coverage(coverage)
|
||
if (identical(expenditure_concept, "total")) {
|
||
.abort_concept_not_aggregatable("cog_peer_compare")
|
||
}
|
||
if (!is.character(target_govid) || length(target_govid) != 1L) {
|
||
cli::cli_abort("`target_govid` must be a length-1 character string.")
|
||
}
|
||
cohort_year <- if (is.data.frame(peers)) {
|
||
ay <- attr(peers, "cohort_year")
|
||
if (is.null(ay)) NA_integer_ else as.integer(ay)
|
||
} else {
|
||
NA_integer_
|
||
}
|
||
pop_range <- if (is.data.frame(peers)) attr(peers, "pop_range") else NULL
|
||
is_ratio <- if (is.data.frame(peers)) attr(peers, "is_ratio") else NULL
|
||
peer_govids <- if (is.data.frame(peers)) {
|
||
as.character(peers$canonical_govid)
|
||
} else {
|
||
as.character(peers)
|
||
}
|
||
peer_govids <- peer_govids[!is.na(peer_govids) & nzchar(peer_govids)]
|
||
all_govids <- unique(c(target_govid, peer_govids))
|
||
|
||
years <- .apply_census_years(years, coverage, "cog_peer_compare")
|
||
|
||
r <- cog_spending(all_govids, years, category, per_capita, adjust_to_year)
|
||
r$role <- ifelse(r$canonical_govid == target_govid, "target", "peer")
|
||
|
||
# The target is exempt from balancing: it is the subject of the comparison,
|
||
# not a member of the cohort being balanced, and dropping it would leave a
|
||
# peer comparison with nothing to compare. Filtering happens BEFORE the
|
||
# quantiles below, so a "consistent" cohort's summary rows describe that
|
||
# cohort rather than the unbalanced one.
|
||
if (identical(coverage, "consistent")) {
|
||
r <- .filter_consistent(r, years, keep_ids = target_govid)
|
||
}
|
||
|
||
value_col <- .peer_value_col(per_capita, adjust_to_year)
|
||
|
||
summary_rows <- .peer_summary_rows(r, value_col)
|
||
out <- dplyr::bind_rows(r, summary_rows)
|
||
rank_val <- .peer_target_rank(r, target_govid, years, value_col)
|
||
out$target_rank <- ifelse(out$role == "target", rank_val, NA_integer_)
|
||
out$cohort_year <- cohort_year
|
||
|
||
prov <- attr(r, "provenance") %||% list()
|
||
prov$verb <- "cog_peer_compare"
|
||
prov$call <- paste(deparse(call), collapse = " ")
|
||
prov$peer_count <- length(peer_govids)
|
||
prov$cohort_year <- cohort_year
|
||
prov$cohort_govids <- peer_govids
|
||
prov$pop_range <- pop_range
|
||
prov$is_ratio <- is_ratio
|
||
prov$target <- list(
|
||
canonical_govid = target_govid,
|
||
gov_name = unique(r$gov_name[r$role == "target"])
|
||
)
|
||
# Counted over PEER rows only, against the cohort size: "3 of your 15 peers
|
||
# reported in FY2019". Including the target would inflate every count by one
|
||
# and make a cohort that has entirely stopped reporting look non-empty.
|
||
prov$coverage_mode <- coverage
|
||
prov$coverage <- .coverage_table(
|
||
out, years, length(peer_govids),
|
||
rows = r[r$role == "peer", , drop = FALSE]
|
||
)
|
||
attr(out, "provenance") <- prov
|
||
out
|
||
}
|
||
|
||
#' @noRd
|
||
.peer_value_col <- function(per_capita, adjust_to_year) {
|
||
if (!is.null(adjust_to_year)) {
|
||
if (isTRUE(per_capita)) "amt_per_capita_real" else "amt_real"
|
||
} else {
|
||
if (isTRUE(per_capita)) "amt_per_capita_nominal" else "amt_nominal"
|
||
}
|
||
}
|
||
|
||
#' @noRd
|
||
.peer_summary_rows <- function(r, value_col) {
|
||
peer_rows <- r[r$role == "peer", , drop = FALSE]
|
||
if (nrow(peer_rows) == 0L) {
|
||
return(r[integer(0), , drop = FALSE])
|
||
}
|
||
s <- peer_rows |>
|
||
dplyr::group_by(.data$year, .data$spend_subtype, .data$category) |>
|
||
dplyr::reframe(
|
||
q = c("p25", "p50", "p75"),
|
||
value = stats::quantile(
|
||
.data[[value_col]], c(0.25, 0.50, 0.75), na.rm = TRUE
|
||
)
|
||
)
|
||
s$role <- paste0("summary_", s$q)
|
||
s$gov_name <- dplyr::case_when(
|
||
s$q == "p25" ~ "Peer P25",
|
||
s$q == "p50" ~ "Peer median",
|
||
s$q == "p75" ~ "Peer P75",
|
||
TRUE ~ NA_character_
|
||
)
|
||
s$canonical_govid <- NA_character_
|
||
out <- s[, c("year", "canonical_govid", "gov_name",
|
||
"spend_subtype", "category", "role")]
|
||
out[[value_col]] <- s$value
|
||
out
|
||
}
|
||
|
||
#' @noRd
|
||
.peer_target_rank <- function(r, target_govid, years, value_col) {
|
||
if (nrow(r) == 0L) return(NA_integer_)
|
||
latest <- max(as.integer(years))
|
||
latest_rows <- r[r$year == latest &
|
||
r$role %in% c("target", "peer"), , drop = FALSE]
|
||
if (nrow(latest_rows) == 0L) return(NA_integer_)
|
||
ranked <- dplyr::arrange(latest_rows, dplyr::desc(.data[[value_col]]))
|
||
tr <- which(ranked$canonical_govid == target_govid)[1]
|
||
if (length(tr) == 0L || is.na(tr)) NA_integer_ else as.integer(tr)
|
||
}
|