diff --git a/.gitignore b/.gitignore index 7c7683a..cbd2f9b 100644 --- a/.gitignore +++ b/.gitignore @@ -25,3 +25,5 @@ claude.md /Meta/ docs /doc/ +cpsvote.Rcheck/ +cpsvote_0.2.0.tar.gz diff --git a/DESCRIPTION b/DESCRIPTION index 7d3ee9b..6daa767 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -41,7 +41,11 @@ Suggests: ggplot2, usmap, ggthemes, - tweenr + tweenr, + kableExtra, + ggrepel, + googlesheets4, + readxl VignetteBuilder: knitr Roxygen: list(markdown = TRUE) RoxygenNote: 7.3.3 diff --git a/NAMESPACE b/NAMESPACE index dc34252..821c623 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,6 +1,7 @@ # Generated by roxygen2: do not edit by hand export("%>%") +export(cps_check_codes) export(cps_data_dir) export(cps_docs_dir) export(cps_download_data) diff --git a/NEWS.md b/NEWS.md index cccb841..8b5f1b4 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,6 @@ # cpsvote 0.2 +- Add `cps_check_codes()`, which warns when a response documented in the codebook (e.g. the `-9` "No Response" code for `PES1`/`VRS_VOTE`) is entirely absent from a given year's data. `cps_recode_vote()` now calls this automatically. This catches a real discontinuity: the Census Bureau's public-use file stops including any `-9` respondents starting in 2022, without documenting the change (see `vignette("voting")` for details, h/t Michael Hanmer) - Include new 2020, 2022, and 2024 CPS VRS data - Note: 2024 microdata downloaded from Census Bureau (not NBER); VEP reweighting data from University of Florida Election Lab - Update maintainer and copyright info diff --git a/R/cps_check_codes.R b/R/cps_check_codes.R new file mode 100644 index 0000000..c9ecdbc --- /dev/null +++ b/R/cps_check_codes.R @@ -0,0 +1,71 @@ +#' check that expected response codes actually appear in the data +#' +#' @description The Census Bureau occasionally changes how it codes +#' nonresponse between survey years without documenting the change. For +#' example, the `-9` ("No Response") code for `PES1`/`VRS_VOTE` appears in +#' every CPS VRS release from 1994-2020, but is entirely absent starting in +#' 2022 (the codebook still lists it, but no respondent is actually coded +#' that way -- see `vignette("voting")` for details). Because +#' [cps_label()] and [cps_recode_vote()] do not error when an expected +#' response is simply absent, this kind of change can silently shift what +#' population is excluded from `cps_turnout` and `hurachen_turnout` without +#' any warning. This function compares the responses documented in `factors` +#' for a given year against what is actually observed in `data`, and warns +#' when a documented response is completely unobserved for a year where it +#' should exist. +#' @param data labelled CPS data (i.e. the output of [cps_label()]), +#' containing `vote_col` and a `YEAR` column +#' @param vote_col which column to check +#' @param check_labels which labels (as they appear in `factors$value`, +#' case-insensitive) to check for. Defaults to the three categories the +#' Census treats as nonvoters: "Don't Know", "Refused", and "No Response" +#' @param factors A data frame containing the label codes to be applied, +#' as in [cps_label()] +#' @param names_col Which column of `factors` contains the column names of +#' `data` +#' @return `data`, invisibly. This function is called for its side effect of +#' issuing warnings. +#' @examples cps_check_codes(cps_label(cps_2024_10k)) +#' +#' @export +cps_check_codes <- function(data, + vote_col = "VRS_VOTE", + check_labels = c("DON'T KNOW", "REFUSED", "NO RESPONSE"), + factors = cpsvote::cps_factors, + names_col = "new_name") { + + if (!(vote_col %in% colnames(data)) || !("YEAR" %in% colnames(data))) { + return(invisible(data)) + } + + # numeric (pre-`cps_label()`) data isn't labelled yet, so there's nothing + # to compare against `factors$value` + if (is.numeric(data[[vote_col]])) { + return(invisible(data)) + } + + check_labels <- toupper(check_labels) + + for (yr in sort(unique(data$YEAR))) { + documented <- toupper(factors$value[factors[[names_col]] == vote_col & factors$year == yr]) + expected <- intersect(check_labels, documented) + if (length(expected) == 0) next + + observed <- toupper(as.character(data[[vote_col]][data$YEAR == yr])) + missing <- expected[!(expected %in% observed)] + + if (length(missing) > 0) { + warning( + "cps_check_codes: for YEAR = ", yr, ", column '", vote_col, "', the codebook ", + "documents the response(s) ", paste(sQuote(missing), collapse = ", "), + " but none of them appear in the data. This usually means the Census Bureau ", + "changed its coding of nonresponse for this year (see vignette(\"voting\") for ", + "the 2022+ `-9`/No Response discontinuity), and turnout estimates from ", + "`cps_recode_vote()` / `cps_reweight_turnout()` may be biased as a result.", + call. = FALSE + ) + } + } + + invisible(data) +} diff --git a/R/refactor.R b/R/refactor.R index 42748ad..2f23c48 100644 --- a/R/refactor.R +++ b/R/refactor.R @@ -148,20 +148,26 @@ cps_refactor <- function(data, move_levels = TRUE) { #' where multiple categories map to "No", and one which follows the guidelines #' from Hur & Achen (2013) of setting these categories to `NA`. See the Vignette #' for more information on this process. +#' @details Before recoding, this function calls [cps_check_codes()], which +#' warns if any of `items` is documented in the codebook for a given year but +#' never actually observed in the data -- as happens with the `-9` ("No +#' Response") code starting in 2022. See `vignette("voting")` for details. #' @param data the input data set #' @param vote_col which column contains the voting variable -#' @param items which items should be "No" in the CPS coding and `NA` in the +#' @param items which items should be "No" in the CPS coding and `NA` in the #' Hur & Achen coding -#' +#' #' @return `data` with two columns attached, `cps_turnout` and `hurachen_turnout`, #' voting variables recoded according to the process above #' @importFrom rlang .data #' @examples cps_recode_vote(cps_refactor(cps_label(cps_2016_10k))) -#' +#' #' @export -cps_recode_vote <- function(data, +cps_recode_vote <- function(data, vote_col = "VRS_VOTE", items = c("DON'T KNOW", "REFUSED", "NO RESPONSE")) { + cps_check_codes(data, vote_col = vote_col, check_labels = items) + if (is.numeric(data[[vote_col]])) { output <- dplyr::mutate(data, cps_turnout = dplyr::case_when( diff --git a/man/cps_check_codes.Rd b/man/cps_check_codes.Rd new file mode 100644 index 0000000..3329893 --- /dev/null +++ b/man/cps_check_codes.Rd @@ -0,0 +1,53 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/cps_check_codes.R +\name{cps_check_codes} +\alias{cps_check_codes} +\title{check that expected response codes actually appear in the data} +\usage{ +cps_check_codes( + data, + vote_col = "VRS_VOTE", + check_labels = c("DON'T KNOW", "REFUSED", "NO RESPONSE"), + factors = cpsvote::cps_factors, + names_col = "new_name" +) +} +\arguments{ +\item{data}{labelled CPS data (i.e. the output of \code{\link[=cps_label]{cps_label()}}), +containing \code{vote_col} and a \code{YEAR} column} + +\item{vote_col}{which column to check} + +\item{check_labels}{which labels (as they appear in \code{factors$value}, +case-insensitive) to check for. Defaults to the three categories the +Census treats as nonvoters: "Don't Know", "Refused", and "No Response"} + +\item{factors}{A data frame containing the label codes to be applied, +as in \code{\link[=cps_label]{cps_label()}}} + +\item{names_col}{Which column of \code{factors} contains the column names of +\code{data}} +} +\value{ +\code{data}, invisibly. This function is called for its side effect of +issuing warnings. +} +\description{ +The Census Bureau occasionally changes how it codes +nonresponse between survey years without documenting the change. For +example, the \code{-9} ("No Response") code for \code{PES1}/\code{VRS_VOTE} appears in +every CPS VRS release from 1994-2020, but is entirely absent starting in +2022 (the codebook still lists it, but no respondent is actually coded +that way -- see \code{vignette("voting")} for details). Because +\code{\link[=cps_label]{cps_label()}} and \code{\link[=cps_recode_vote]{cps_recode_vote()}} do not error when an expected +response is simply absent, this kind of change can silently shift what +population is excluded from \code{cps_turnout} and \code{hurachen_turnout} without +any warning. This function compares the responses documented in \code{factors} +for a given year against what is actually observed in \code{data}, and warns +when a documented response is completely unobserved for a year where it +should exist. +} +\examples{ +cps_check_codes(cps_label(cps_2024_10k)) + +} diff --git a/man/cps_recode_vote.Rd b/man/cps_recode_vote.Rd index 7708db8..63dd6ca 100644 --- a/man/cps_recode_vote.Rd +++ b/man/cps_recode_vote.Rd @@ -33,6 +33,12 @@ where multiple categories map to "No", and one which follows the guidelines from Hur & Achen (2013) of setting these categories to \code{NA}. See the Vignette for more information on this process. } +\details{ +Before recoding, this function calls \code{\link[=cps_check_codes]{cps_check_codes()}}, which +warns if any of \code{items} is documented in the codebook for a given year but +never actually observed in the data -- as happens with the \code{-9} ("No +Response") code starting in 2022. See \code{vignette("voting")} for details. +} \examples{ cps_recode_vote(cps_refactor(cps_label(cps_2016_10k))) diff --git a/sandbox/make_massachusetts_gif.R b/sandbox/make_massachusetts_gif.R new file mode 100644 index 0000000..9ce2f47 --- /dev/null +++ b/sandbox/make_massachusetts_gif.R @@ -0,0 +1,25 @@ +# Animated ternary plot of Massachusetts's vote mode share over time +# (1996-2024), all years. A single dot labeled "MA", with a trail showing +# its path through the ternary space over time. + +library(here) +library(dplyr) +source(here("sandbox", "ternary_gif_lib.R")) + +vote_mode <- readRDS(here("sandbox", "vote_mode_data.rds")) %>% + filter(STATE == "MA") + +years_seq <- seq(1996, 2024, 2) + +make_vote_mode_gif( + vote_mode = vote_mode, + years_seq = years_seq, + output_gif = here("sandbox", "massachusetts_vote_mode.gif"), + title = "The Move Away From Election Day Voting in Massachusetts", + subtitle = "Share of votes cast by mode, 1996-2024", + highlight_states = c("MA"), + highlight_color = "red", + trail = TRUE +) + +cat("Wrote", here("sandbox", "massachusetts_vote_mode.gif"), "\n") diff --git a/sandbox/make_midterm_gif.R b/sandbox/make_midterm_gif.R new file mode 100644 index 0000000..ef82dbc --- /dev/null +++ b/sandbox/make_midterm_gif.R @@ -0,0 +1,19 @@ +# Animated ternary plot of vote mode share by state, midterm election years +# only (1998-2022). One dot per state per year. + +library(here) +source(here("sandbox", "ternary_gif_lib.R")) + +vote_mode <- readRDS(here("sandbox", "vote_mode_data.rds")) + +years_seq <- seq(1998, 2022, 4) + +make_vote_mode_gif( + vote_mode = vote_mode, + years_seq = years_seq, + output_gif = here("sandbox", "midterm_vote_mode.gif"), + title = "The Move Away From Election Day Voting in America", + subtitle = "Midterm elections: share of votes cast by mode, 1998-2022" +) + +cat("Wrote", here("sandbox", "midterm_vote_mode.gif"), "\n") diff --git a/sandbox/make_presidential_gif.R b/sandbox/make_presidential_gif.R new file mode 100644 index 0000000..85d62ef --- /dev/null +++ b/sandbox/make_presidential_gif.R @@ -0,0 +1,19 @@ +# Animated ternary plot of vote mode share by state, presidential election +# years only (1996-2024). One dot per state per year. + +library(here) +source(here("sandbox", "ternary_gif_lib.R")) + +vote_mode <- readRDS(here("sandbox", "vote_mode_data.rds")) + +years_seq <- seq(1996, 2024, 4) + +make_vote_mode_gif( + vote_mode = vote_mode, + years_seq = years_seq, + output_gif = here("sandbox", "presidential_vote_mode.gif"), + title = "The Move Away From Election Day Voting in America", + subtitle = "Presidential elections: share of votes cast by mode, 1996-2024" +) + +cat("Wrote", here("sandbox", "presidential_vote_mode.gif"), "\n") diff --git a/sandbox/massachusetts_vote_mode.gif b/sandbox/massachusetts_vote_mode.gif new file mode 100644 index 0000000..3dd8d75 Binary files /dev/null and b/sandbox/massachusetts_vote_mode.gif differ diff --git a/sandbox/midterm_vote_mode.gif b/sandbox/midterm_vote_mode.gif new file mode 100644 index 0000000..ce365b6 Binary files /dev/null and b/sandbox/midterm_vote_mode.gif differ diff --git a/sandbox/prep_vote_mode_data.R b/sandbox/prep_vote_mode_data.R new file mode 100644 index 0000000..fa19b29 --- /dev/null +++ b/sandbox/prep_vote_mode_data.R @@ -0,0 +1,32 @@ +# Purpose: Build the weighted YEAR x STATE vote-mode share table used by the +# three ternary-plot GIF scripts (presidential years, midterm years, MA-only). +# This is the same data-prep pipeline as sandbox/animate_snowglobe.R, factored +# out so it only needs to run once. + +library(dplyr) +library(tidyr) +library(srvyr) +library(here) + +devtools::load_all(here(), quiet = TRUE) + +options(cpsvote.datadir = here("cps_data")) + +cps <- cps_load_basic(years = seq(1996, 2024, 2)) + +cps_weighted <- cps %>% + filter(YEAR > 1995) %>% + srvyr::as_survey_design(weights = turnout_weight) + +vote_mode <- cps_weighted %>% + select(YEAR, STATE, VRS_VOTEMETHOD_CON) %>% + filter(if_all(everything(), ~ !is.na(.x))) %>% + group_by(YEAR, STATE, VRS_VOTEMETHOD_CON) %>% + summarize(survey_mean(na.rm = TRUE)) %>% + select(-ends_with('_se')) %>% + pivot_wider(id_cols = c("YEAR", "STATE"), names_from = "VRS_VOTEMETHOD_CON", values_from = "coef", + values_fill = 0) + +saveRDS(vote_mode, here("sandbox", "vote_mode_data.rds")) + +cat("Saved", nrow(vote_mode), "rows across years:", paste(sort(unique(vote_mode$YEAR)), collapse = ", "), "\n") diff --git a/sandbox/presidential_vote_mode.gif b/sandbox/presidential_vote_mode.gif new file mode 100644 index 0000000..e74cfdd Binary files /dev/null and b/sandbox/presidential_vote_mode.gif differ diff --git a/sandbox/ternary_gif_lib.R b/sandbox/ternary_gif_lib.R new file mode 100644 index 0000000..da9c979 --- /dev/null +++ b/sandbox/ternary_gif_lib.R @@ -0,0 +1,112 @@ +# Shared helpers for rendering animated ternary plots of vote mode share +# (Election Day / Early In-Person / Mail-VBM), in the same visual style as +# img/MDvote_mode.gif. Used by make_presidential_gif.R, make_midterm_gif.R, +# and make_massachusetts_gif.R. + +library(dplyr) +library(tweenr) +library(ggtern) +library(magick) +library(here) + +#' Tween a YEAR x STATE vote-mode table across a sequence of anchor years, +#' holding on each anchor and smoothly interpolating between them. +build_tweened_frames <- function(vote_mode, years_seq, nframes_hold = 10, nframes_transition = 10) { + vote_mode <- ungroup(vote_mode) + years_seq <- sort(unique(years_seq)) + gap <- min(diff(years_seq)) + + frames <- filter(vote_mode, YEAR == years_seq[1]) %>% + keep_state(nframes_hold) + + for (yr in years_seq[-1]) { + frames <- frames %>% + tween_state(filter(vote_mode, YEAR == yr), 'linear', id = STATE, nframes = nframes_transition) %>% + keep_state(nframes_hold) + } + + frames %>% + mutate(YEAR = floor((YEAR - years_seq[1]) / gap) * gap + years_seq[1]) +} + +#' Render one ternary-plot frame per unique .frame value in `tweened_data` and +#' write the resulting PNGs to `frame_dir`. +render_frames <- function(tweened_data, frame_dir, title, subtitle, + highlight_states = c(), highlight_color = "red", + trail = FALSE, width = 4.25, height = 3.5) { + dir.create(frame_dir, showWarnings = FALSE, recursive = TRUE) + + tweened_data <- tweened_data %>% + ungroup() %>% + mutate(highlight = STATE %in% highlight_states) + + trail_data <- if (trail && length(highlight_states) == 1) { + tweened_data %>% + filter(STATE == highlight_states[1]) %>% + select(.frame, `ELECTION DAY`, `BY MAIL`, EARLY) + } else { + NULL + } + + for (frame in unique(tweened_data$.frame)) { + frame_data <- filter(tweened_data, .frame == frame) + yr <- unique(frame_data$YEAR) + + p <- ggplot(frame_data, + aes(y = `ELECTION DAY`, x = `BY MAIL`, z = EARLY, label = STATE, colour = highlight)) + + geom_text(vjust = 0.5, hjust = 0.5, size = 1.5, position = "identity") + + scale_colour_manual(values = c("FALSE" = "black", "TRUE" = highlight_color), guide = "none") + + coord_tern() + + labs(x = "Mail", + y = "Election Day", + z = "Early", + title = title, + subtitle = paste0(subtitle, "\nYear: ", yr)) + + theme_classic() + + theme(plot.title = element_text(hjust = 0.5, size = 8), + plot.subtitle = element_text(hjust = 0.5, size = 6), + axis.title = element_text(size = 5), + axis.text = element_text(size = 4)) + + if (!is.null(trail_data)) { + p <- p + geom_path(data = filter(trail_data, .frame <= frame), + aes(y = `ELECTION DAY`, x = `BY MAIL`, z = EARLY, group = 1), + colour = highlight_color, inherit.aes = FALSE, linewidth = 0.5) + } + + ggsave(plot = p, + filename = file.path(frame_dir, paste0('frame', stringr::str_pad(frame, width = 4, pad = "0"), '.png')), + width = width, height = height) + } +} + +#' Stitch the PNGs in `frame_dir` into a GIF at `output_gif`. +write_gif <- function(frame_dir, output_gif, fps = 10) { + list.files(path = frame_dir, pattern = "\\.png$", full.names = TRUE) %>% + sort() %>% + purrr::map(image_read) %>% + image_join() %>% + image_animate(fps = fps) %>% + image_write(output_gif) +} + +#' End-to-end: tween, render frames, write gif, clean up frame_dir. +make_vote_mode_gif <- function(vote_mode, years_seq, output_gif, title, subtitle, + highlight_states = c(), highlight_color = "red", + trail = FALSE, fps = 10, + nframes_hold = 10, nframes_transition = 10, + frame_dir = NULL, cleanup = TRUE) { + if (is.null(frame_dir)) { + frame_dir <- file.path(tempdir(), paste0("frames_", tools::file_path_sans_ext(basename(output_gif)))) + } + + tweened <- build_tweened_frames(vote_mode, years_seq, nframes_hold, nframes_transition) + render_frames(tweened, frame_dir, title, subtitle, + highlight_states = highlight_states, highlight_color = highlight_color, + trail = trail) + write_gif(frame_dir, output_gif, fps = fps) + + if (cleanup) unlink(frame_dir, recursive = TRUE) + + invisible(output_gif) +} diff --git a/vignettes/articles/validation_2022_turnout.Rmd b/vignettes/articles/validation_2022_turnout.Rmd new file mode 100644 index 0000000..8a8b22e --- /dev/null +++ b/vignettes/articles/validation_2022_turnout.Rmd @@ -0,0 +1,402 @@ +--- +title: "Validating cpsvote's 2022 Turnout Estimates" +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + echo = TRUE, + message = FALSE, + warning = FALSE, + fig.width = 7, + fig.height = 5 +) +# resolved relative to this file's own location (vignettes/articles/), which +# knitr guarantees as the working directory during rendering -- deliberately +# not using here::here(), which resolves inconsistently under pkgdown's build +# subprocess +cps_data_dir <- normalizePath(file.path("..", "..", "cps_data"), mustWork = FALSE) + +# This article downloads full raw CPS microdata (several years, tens of +# thousands of rows each) plus external McDonald VEP and Census Table 4a +# data over the network. None of that is available on the pkgdown CI +# runner -- there's no cached cps_data/ there, and pulling it fresh from +# Census/NBER/Google Sheets on every CI build would be slow and fragile +# (the pre-2020 NBER mirror cps_read() falls back to is no longer reliably +# reachable). So: skip the data-heavy chunks entirely when running in CI, +# or when the local raw-data cache isn't present, and say so explicitly +# rather than silently producing an empty article. +IN_CI <- identical(Sys.getenv("CI"), "true") +HAS_CPS_DATA <- dir.exists(cps_data_dir) && + length(list.files(cps_data_dir, pattern = "\\.zip$")) > 0 +RUN_LIVE <- !IN_CI && HAS_CPS_DATA + +knitr::opts_chunk$set(eval = RUN_LIVE) +``` + +```{r, echo = FALSE, eval = !RUN_LIVE, results = "asis"} +cat( + "> **Note:** this article computes its tables and figures from full raw", + "CPS microdata files and live downloads (Census Bureau, Google Sheets)", + "that aren't available in the environment used to build this site", + "automatically. The code below is shown for reference; render this", + "article locally with a populated `cps_data/` cache (see `cps_download_data()`)", + "to reproduce its output.\n" +) +``` + +```{r setup} +library(cpsvote) +library(dplyr) +library(srvyr) +library(ggplot2) +library(knitr) +library(kableExtra) +``` + +## Overview + +This article validates `cpsvote`'s 2022 turnout estimates against two independent, authoritative sources, and documents a data issue discovered while doing so. It's an internal validation exercise (originally run as `sandbox/validate_2022_turnout.Rmd`), published here as a pkgdown article rather than a package vignette because it downloads external data, takes a while to run, and produces output (long state-by-state tables, several figures) that's more useful as a standing reference than as CRAN-checked documentation. + +Two comparisons are run: + +1. **CPS standard weight vs. McDonald VEP turnout** — checks that `cpsvote`'s numbers reproduce the well-documented CPS overreporting bias, and shows what the Hur-Achen reweighting (`turnout_weight`) buys you. +2. **`cpsvote` vs. Census Table 4a** — a stricter consistency check. If `cpsvote` reproduces the Census Bureau's own published methodology (standard weight, `cps_turnout` coding), state estimates should match Table 4a's numbers almost exactly. They don't, out of the box, and tracking down why turned up a real (if narrow) limitation in using `PWSSWGT` for citizen-only subgroup analysis. + +A third section documents a related, more consequential issue found afterward: starting in 2022, the Census Bureau's public files stop using the `-9` ("No Response") code for the vote question at all, which silently affects `cps_turnout` for 2022 and 2024. + +Throughout, `hurachen_turnout` is used for Comparison 1 and `cps_turnout` for Comparison 2, since Comparison 2 is explicitly reproducing the Census's own reporting convention (which codes Don't Know/Refused/No Response as "No" rather than dropping them). + +```{r load-cps, cache = TRUE} +cps_2022 <- cps_load_basic(years = 2022, datadir = cps_data_dir) %>% + as_survey_design(weights = WEIGHT) # standard Census weight, NOT turnout_weight +``` + +## Comparison 1: CPS Standard Weight vs. McDonald VEP + +The CPS is well known to overreport turnout relative to the actual count of ballots cast, because non-voters are more likely than voters to refuse to participate in (or to complete) the survey. Michael McDonald's Voting-Eligible Population (VEP) turnout figures, built from official state vote counts, are the standard benchmark against which this overreporting is measured. This is precisely the bias that the Hur & Achen (2013) reweighting procedure — implemented here as `cps_reweight_turnout()` — is designed to correct. + +```{r load-mcdonald, cache = TRUE} +googlesheets4::gs4_deauth() + +gid_2022 <- "17iSZIBPP6jj0C2wcqpym9wNuNCtTiqM_tMweMjmvces" + +mcdonald_2022 <- googlesheets4::read_sheet( + gid_2022, + range = "A3:N54", + col_names = c( + "state_name", "highestoffice", "status", "results_source", + "pct_highestoffice_vep", "vep", "vap", "pct_noncitizen", + "prison", "probation", "parole", "total_ineligible_felon", + "overseas_eligible", "state_abb" + ) +) %>% + filter(state_abb != "US") %>% + transmute( + STATE = state_abb, + state_name = stringr::str_remove_all(state_name, "\\*$"), + mcd_turnout_pct = round(as.numeric(pct_highestoffice_vep) * 100, 1) + ) + +cps_turnout <- cps_2022 %>% + group_by(STATE) %>% + summarize(cps_turnout_pct = survey_mean(hurachen_turnout == "YES", na.rm = TRUE)) %>% + select(STATE, cps_turnout_pct) %>% + mutate(cps_turnout_pct = round(cps_turnout_pct * 100, 1)) + +comparison <- left_join(mcdonald_2022, cps_turnout, by = "STATE") %>% + mutate(difference = cps_turnout_pct - mcd_turnout_pct) %>% + arrange(state_name) +``` + +### State-by-state comparison + +```{r table} +comparison %>% + select( + State = state_name, + `McDonald VEP (%)` = mcd_turnout_pct, + `CPS Std. Weight (%)` = cps_turnout_pct, + `Difference (pp)` = difference + ) %>% + kable(align = c("l", "r", "r", "r"), + caption = "2022 State Voter Turnout: CPS Standard Weight vs. McDonald VEP Estimates") %>% + kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) %>% + scroll_box(height = "400px") %>% + column_spec(4, bold = TRUE, + color = ifelse(abs(comparison$difference) > 5, "red", "black")) +``` + +```{r summary} +summary_stats <- comparison %>% + summarize( + `Mean difference (pp)` = round(mean(difference, na.rm = TRUE), 2), + `Median difference (pp)` = round(median(difference, na.rm = TRUE), 2), + `SD of difference (pp)` = round(sd(difference, na.rm = TRUE), 2), + `States overestimating` = sum(difference > 0, na.rm = TRUE), + `States > 5pp off` = sum(abs(difference) > 5, na.rm = TRUE), + `States (n)` = n() + ) %>% + tidyr::pivot_longer(everything(), names_to = "Statistic", values_to = "Value") + +kable(summary_stats, caption = "Summary of CPS (standard weight) vs. McDonald differences (percentage points)") %>% + kable_styling(bootstrap_options = c("striped", "condensed"), full_width = FALSE) +``` + +Every one of the 51 states/DC overestimates McDonald's VEP turnout using the standard weight, with a mean gap of about 12.5 percentage points (median 12.2) and 50 of 51 states off by more than 5 points. This lands squarely in the range of the CPS overreporting bias documented elsewhere (Hur & Achen report the No Response category alone knocking nearly 10 points off the 2008 estimate once properly coded); the exact size of the gap varies by year and by whether Don't Know/Refused/No Response are coded as missing (as here) or as "No" (Comparison 2 below), but the direction and rough magnitude are consistent with the literature. + +```{r scatter, fig.cap = "CPS (standard weight) vs. McDonald VEP turnout by state, 2022. The dashed line is the 45-degree line of perfect agreement. Points above the line indicate CPS overestimation."} +comparison %>% + mutate(label = ifelse(abs(difference) > 15, STATE, NA)) %>% + ggplot(aes(x = mcd_turnout_pct, y = cps_turnout_pct)) + + geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray50") + + geom_point(aes(color = abs(difference) > 5), size = 2) + + ggrepel::geom_text_repel(aes(label = label), size = 2.5, color = "red", + na.rm = TRUE, max.overlaps = 20) + + scale_color_manual(values = c("FALSE" = "steelblue", "TRUE" = "red"), guide = "none") + + scale_x_continuous(labels = scales::percent_format(scale = 1), limits = c(20, 70)) + + scale_y_continuous(labels = scales::percent_format(scale = 1), limits = c(20, 80)) + + labs( + x = "McDonald VEP Turnout (%)", y = "CPS Standard Weight Turnout (%)", + title = "CPS vs. McDonald Turnout Estimates, 2022", + subtitle = "Using standard Census survey weight; points above the line = CPS overestimates" + ) + + theme_bw() +``` + +Figure 1 shows virtually every state above the 45-degree line — the standard Census weight (`WEIGHT`, i.e. `PWSSWGT`) systematically overstates turnout because it does not adjust for differential non-response between voters and non-voters. This is exactly the gap the Hur-Achen procedure targets. + +```{r ha-cps, cache = TRUE} +cps_ha <- cps_load_basic(years = 2022, datadir = cps_data_dir) %>% + as_survey_design(weights = turnout_weight) + +cps_ha_turnout <- cps_ha %>% + group_by(STATE) %>% + summarize(ha_turnout_pct = survey_mean(hurachen_turnout == "YES", na.rm = TRUE)) %>% + select(STATE, ha_turnout_pct) %>% + mutate(ha_turnout_pct = round(ha_turnout_pct * 100, 1)) + +comparison_ha <- left_join(mcdonald_2022, cps_ha_turnout, by = "STATE") %>% + mutate(difference_ha = ha_turnout_pct - mcd_turnout_pct) +``` + +```{r ha-scatter, fig.cap = "CPS (Hur-Achen reweighted) vs. McDonald VEP turnout by state, 2022. Reweighting pulls points markedly closer to the 45-degree line compared with Figure 1."} +comparison_ha %>% + mutate(label = ifelse(abs(difference_ha) > 5, STATE, NA)) %>% + ggplot(aes(x = mcd_turnout_pct, y = ha_turnout_pct)) + + geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray50") + + geom_point(aes(color = abs(difference_ha) > 5), size = 2) + + ggrepel::geom_text_repel(aes(label = label), size = 2.5, color = "red", + na.rm = TRUE, max.overlaps = 20) + + scale_color_manual(values = c("FALSE" = "steelblue", "TRUE" = "red"), guide = "none") + + scale_x_continuous(labels = scales::percent_format(scale = 1), limits = c(20, 70)) + + scale_y_continuous(labels = scales::percent_format(scale = 1), limits = c(20, 80)) + + labs( + x = "McDonald VEP Turnout (%)", y = "CPS Hur-Achen Weight Turnout (%)", + title = "CPS (Hur-Achen weight) vs. McDonald Turnout Estimates, 2022", + subtitle = "Using Hur-Achen reweighted survey weight; points near the line = close agreement" + ) + + theme_bw() +``` + +Comparing Figures 1 and 2 is the core demonstration of what `cpsvote` buys you: `turnout_weight` (built by `cps_reweight_turnout()`) post-stratifies each state so that the reweighted CPS sample turnout matches McDonald's VEP figure by construction — which is why this comparison, on its own, can't be used to validate the package's *implementation*. That's what Comparison 2 is for. + +## Comparison 2: cpsvote vs. Census Table 4a + +This is a stricter test. If `cpsvote` correctly reproduces the Census Bureau's *own* methodology — standard `WEIGHT`, plus the Census's own coding of Don't Know/Refused/No Response as "No" (`cps_turnout`, not `hurachen_turnout`) — then state-level estimates should match the Bureau's published [Table 4a](https://www.census.gov/data/tables/time-series/demo/voting-and-registration/p20-586.html) ("Reported Voting and Registration of the Total Voting-Age Population, for States: November 2022") almost exactly. Any real discrepancy points to a problem in the package. + +```{r census-data, cache = TRUE} +census_url <- "https://www2.census.gov/programs-surveys/cps/tables/p20/586/vote04a_2022.xlsx" +census_file <- tempfile(fileext = ".xlsx") +download.file(census_url, census_file, quiet = TRUE, mode = "wb") + +census_raw <- suppressMessages(readxl::read_excel(census_file, col_names = FALSE, skip = 4)) + +state_lookup <- tibble::tibble( + state_name = c(toupper(state.name), "DISTRICT OF COLUMBIA"), + state_abb = c(state.abb, "DC") +) + +# column ...12 = "Voted - Percent of citizen population" +census_table4a <- census_raw %>% + select(state_name = `...1`, voted_pct_citizen = `...12`) %>% + filter(!is.na(state_name), state_name %in% state_lookup$state_name) %>% + left_join(state_lookup, by = "state_name") %>% + mutate(census_pct = as.numeric(voted_pct_citizen)) %>% + select(STATE = state_abb, state_name, census_pct) + +# column ...3 = citizen voting-age population (thousands) +citizen_vap_targets <- census_raw %>% + select(state_name = `...1`, citizen_pop_thou = `...3`) %>% + filter(!is.na(state_name), state_name %in% state_lookup$state_name) %>% + left_join(state_lookup, by = "state_name") %>% + transmute(STATE = state_abb, census_vap_persons = as.numeric(citizen_pop_thou) * 1000) +``` + +### The naive comparison, and the denominator problem + +A first pass, computing `cps_turnout` state proportions directly from the survey design (the same way Comparison 1 computed `hurachen_turnout` proportions), shows every state overestimating Table 4a too — which is suspicious, because this comparison is supposed to reproduce the Census's own numbers. + +The tempting hypothesis is that `cpsvote` uses the wrong population denominator (e.g. including non-citizens). Tracing through `cps_label()` rules this out: code `-1` for `PES1` is labelled `"Not in Universe"`, which `cps_label()` upper-cases and sets to `NA` by default — this correctly drops non-citizens and respondents under 18 from the denominator of `survey_mean(..., na.rm = TRUE)`. + +The actual cause is the weight itself. `WEIGHT` (i.e. `PWSSWGT`) is calibrated by the Census Bureau to the **total** civilian noninstitutional population, including non-citizens. When the analysis is restricted to VRS respondents (citizens 18+) and their weights are summed, that total falls short of the Bureau's own, independently-derived citizen voting-age population control totals (Table 4a, column 3) — because the weight was never calibrated to that subgroup in the first place. + +```{r denom-check} +vrs_raw <- cps_2022$variables %>% + filter(!is.na(cps_turnout)) %>% + mutate(STATE = as.character(STATE)) + +denom_before <- vrs_raw %>% + group_by(STATE) %>% + summarize(before_persons = round(sum(WEIGHT)), .groups = "drop") + +denom_verify <- left_join(denom_before, citizen_vap_targets, by = "STATE") %>% + mutate(shortfall_k = round((before_persons - census_vap_persons) / 1000), + shortfall_pct = round((before_persons - census_vap_persons) / census_vap_persons * 100, 1)) %>% + arrange(STATE) + +denom_verify %>% + select(STATE, + `cpsvote Denom. (persons)` = before_persons, + `Census Citizen VAP (persons)` = census_vap_persons, + `Shortfall (000s)` = shortfall_k, + `Shortfall (%)` = shortfall_pct) %>% + kable(align = c("l", "r", "r", "r", "r"), + caption = "Denominator before adjustment: cpsvote sum-of-weights vs. Census citizen VAP") %>% + kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) %>% + scroll_box(height = "400px") +``` + +Across all 51 states/DC, `cpsvote`'s implicit denominator falls short of the Census's citizen VAP total, by 11.2% on average. Since the numerator (weighted count of "Yes" respondents) is the same either way, a smaller denominator mechanically inflates the turnout estimate relative to Table 4a in every single state. This is a genuine limitation of using `PWSSWGT` for a citizen-only subgroup analysis, not a coding bug — but it means naive `cps_turnout` estimates from `cpsvote` should not be directly compared to Table 4a without correction. + +### The fix: use the Census's own denominator + +A uniform rescaling of weights (e.g. `survey::postStratify()`) cannot fix this: multiplying every weight in a stratum by a constant `k` leaves any proportion `sum(k*w*x) / sum(k*w)` unchanged, since `k` cancels. The actual fix is to compute the estimate the way Table 4a itself is built — numerator from the CPS survey, but the *denominator* replaced by the Bureau's independent citizen-VAP control total: + +$$\text{Turnout} = \frac{\sum_i w_i \cdot \mathbf{1}[\text{voted}_i]}{\text{Census citizen VAP}}$$ + +```{r post-stratify} +yes_totals <- vrs_raw %>% + group_by(STATE) %>% + summarize(yes_wt = sum(WEIGHT[cps_turnout == "YES"]), .groups = "drop") + +our_census_ps <- left_join(yes_totals, citizen_vap_targets, by = "STATE") %>% + mutate(our_pct = round(yes_wt / census_vap_persons * 100, 1)) %>% + select(STATE, our_pct) + +census_vs_ours <- left_join(census_table4a, our_census_ps, by = "STATE") %>% + mutate(diff = our_pct - census_pct, + label = ifelse(abs(diff) > 3, STATE, NA)) +``` + +```{r census-comparison-table} +census_vs_ours %>% + select( + State = state_name, + `Census Table 4a (%)` = census_pct, + `cpsvote, corrected denom. (%)` = our_pct, + `Difference (pp)` = diff + ) %>% + arrange(State) %>% + kable(align = c("l", "r", "r", "r"), + caption = "2022 State Turnout: cpsvote (Census citizen-VAP denominator) vs. Census Table 4a") %>% + kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) %>% + scroll_box(height = "400px") %>% + column_spec(4, bold = TRUE, color = ifelse(abs(census_vs_ours$diff) > 3, "red", "black")) +``` + +```{r census-scatter, fig.cap = "Figure 2: cpsvote (Census citizen-VAP denominator) vs. Census Table 4a by state, 2022. After substituting the correct denominator, every state falls essentially on the 45-degree line."} +ggplot(census_vs_ours, aes(x = census_pct, y = our_pct)) + + geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray50") + + geom_point(aes(color = abs(diff) > 3), size = 2) + + ggrepel::geom_text_repel(aes(label = label), size = 2.5, color = "red", + na.rm = TRUE, max.overlaps = 20) + + scale_color_manual(values = c("FALSE" = "steelblue", "TRUE" = "red"), guide = "none") + + scale_x_continuous(labels = scales::percent_format(scale = 1), limits = c(30, 75)) + + scale_y_continuous(labels = scales::percent_format(scale = 1), limits = c(30, 75)) + + labs( + x = "Census Table 4a: Voted % of Citizen Population", + y = "cpsvote: Voted % (corrected denominator, Census coding)", + title = "cpsvote vs. Census Table 4a, 2022 (Figure 2)", + subtitle = "Points now lie on the dashed 45° line" + ) + + theme_bw() +``` + +With the corrected denominator, the mean difference across all 51 states/DC drops from 6.6 percentage points (naive `cpsvote` denominator) to **0.0 percentage points**, with a maximum discrepancy of just 0.1pp anywhere in the country — Figure 2 shows every state essentially on the 45-degree line. This includes Colorado, which showed the largest gap (14.3pp) before the fix and looked like it might be a state-specific data quality issue; it turns out to be fully explained by the same denominator effect as every other state, with no residual discrepancy once the citizen-VAP denominator is used. The remaining takeaway is narrow but important: **`cps_turnout`/`hurachen_turnout` state proportions computed directly from `WEIGHT` should not be compared to Table 4a without substituting Table 4a's own citizen-VAP denominator**; the underlying turnout coding in `cpsvote` is correct. + +## The 2022+ No Response Discontinuity + +While reviewing 2022 and 2024 data in response to a [user-reported issue](https://github.com/Reed-EVIC/cpsvote/issues/69) from Dr. Michael Hanmer, a separate and more consequential problem turned up: starting in 2022, the Census Bureau's public-use file stops using the `-9` ("No Response") code for `PES1`/`VRS_VOTE` entirely. + +```{r no-response-codes, cache = TRUE} +years_to_check <- c(2016, 2018, 2020, 2022, 2024) + +code_pcts <- lapply(years_to_check, function(yr) { + d <- cps_read(years = yr, dir = cps_data_dir) + tab <- table(d$VRS_VOTE) + tibble::tibble(year = yr, code = names(tab), pct = round(100 * as.integer(tab) / sum(tab), 1)) +}) %>% + bind_rows() %>% + tidyr::pivot_wider(id_cols = year, names_from = code, values_from = pct, values_fill = 0) + +age_pcts <- lapply(years_to_check, function(yr) { + d <- cps_read(years = yr, dir = cps_data_dir) + d %>% + filter(VRS_VOTE == -1) %>% + summarize(year = yr, n = n(), pct_18plus = round(100 * mean(AGE >= 18), 1)) +}) %>% + bind_rows() +``` + +```{r no-response-table} +code_pcts %>% + select(Year = year, `-9 (No Response)` = `-9`, `-1 (Not in Universe)` = `-1`, + `1 (Yes)` = `1`, `2 (No)` = `2`) %>% + kable(align = "r", caption = "Percent of PES1/VRS_VOTE respondents by code, by year (full raw files)") %>% + kable_styling(bootstrap_options = c("striped", "condensed"), full_width = FALSE) +``` + +`-9` is present in every year through 2020 and entirely absent in 2022 and 2024, coinciding with a roughly 10-point jump in the share of respondents coded `-1`. Because genuine "Not in Universe" cases should overwhelmingly be people who were never eligible to be asked the question (minors, mainly), the age composition of the `-1` group is informative: + +```{r age-table} +age_pcts %>% + select(Year = year, `-1 count` = n, `% of -1 group aged 18+` = pct_18plus) %>% + kable(align = "r", caption = "Age composition of the '-1' (Not in Universe) group, by year") %>% + kable_styling(bootstrap_options = c("striped", "condensed"), full_width = FALSE) +``` + +The adult share of the `-1` group roughly doubles right when `-9` disappears, and the size of the increase is in the same ballpark as the number of `-9` respondents seen in 2020 (about 10,000, median age 49, 100% adult). This strongly suggests the Census Bureau folded "No Response" adults into the generic "Not in Universe" bucket starting in 2022, rather than that these respondents actually vanished from the survey. + +**Why this matters:** `cps_label()` treats `-1` as missing (`NA`) *before* `cps_recode_vote()` ever runs. Through 2020, a `-9` respondent was a real observation, counted as "No" for `cps_turnout` (matching the Census's own convention) and correctly dropped to `NA` for `hurachen_turnout`. Since 2022, that same person (now coded `-1`) is wiped to `NA` at the labelling step and never reaches the recoding logic at all — `cps_turnout`'s denominator silently shrinks, inflating the naive turnout estimate, with no error or warning. `hurachen_turnout` is comparatively less affected, since Don't Know/Refused/No Response were always excluded from it regardless of the exact numeric code, but the population used to calibrate the Hur-Achen reweighting could be subtly affected too. + +As of this package version, `cps_recode_vote()` calls `cps_check_codes()` automatically, which compares the codebook's documented responses for a given year against what's actually observed and warns when one is documented but completely missing — exactly this situation: + +```{r check-codes-demo} +withCallingHandlers( + invisible( + cps_read(years = 2022, dir = cps_data_dir) %>% + cps_label() %>% + cps_refactor() %>% + cps_recode_vote() + ), + warning = function(w) { + cat("Warning:", conditionMessage(w), "\n") + invokeRestart("muffleWarning") + } +) +``` + +See `vignette("voting")` for the full writeup. `cpsvote` does not currently attempt to algorithmically recover the lost `-9` responses (e.g. by using age to split the `-1` bucket back into "Not in Universe" vs. "No Response") — that would be a substantive change to how turnout is coded and needs further discussion before it's adopted. + +## Conclusions and Recommendations + +- **Use `turnout_weight` (Hur-Achen reweighted), not `WEIGHT`, for headline turnout estimates.** The standard weight systematically overreports turnout by around 12 points on average in 2022; the reweighting exists specifically to correct this. +- **Don't compare `cpsvote`'s `WEIGHT`-based state proportions directly to Census Table 4a.** `PWSSWGT` is calibrated to the full civilian population, not the citizen voting-age population Table 4a reports against; the fix is to substitute Table 4a's own citizen-VAP totals as the denominator, not to rescale the survey weight. Once that's done, `cpsvote`'s turnout coding matches Table 4a almost exactly (mean difference 0.0pp, max 0.1pp, across all 51 states/DC) — including for Colorado, which looked like an outlier before the denominator was corrected. +- **Treat 2022 and 2024 `cps_turnout` estimates with caution.** The disappearance of the `-9` "No Response" code starting in 2022 means people who would previously have been counted as non-voters are now silently excluded from the denominator instead, inflating naive turnout estimates for those years specifically. `hurachen_turnout` is less affected. `cps_recode_vote()` will now warn when this happens — don't ignore that warning. +- **When in doubt, validate.** Both of the substantive findings here (the Table 4a denominator issue and the 2022+ No Response discontinuity) were caught by comparing `cpsvote` output against external, authoritative benchmarks rather than by inspecting the package's code in isolation. The same approach is recommended before relying on `cpsvote` estimates for any year not covered by this validation. + +## Data Sources + +- McDonald's 2022 VEP turnout data: `https://docs.google.com/spreadsheets/d/17iSZIBPP6jj0C2wcqpym9wNuNCtTiqM_tMweMjmvces` +- Census Bureau Table 4a for 2022: `https://www2.census.gov/programs-surveys/cps/tables/p20/586/vote04a_2022.xlsx` +- GitHub issue reporting the 2022+ No Response discontinuity: `https://github.com/Reed-EVIC/cpsvote/issues/69` diff --git a/vignettes/voting.Rmd b/vignettes/voting.Rmd index 18a3dc5..d2ee37c 100644 --- a/vignettes/voting.Rmd +++ b/vignettes/voting.Rmd @@ -90,6 +90,18 @@ Professor Michael McDonald of the University of Florida helpfully provides guida **The Hur and Achen corrections have been integrated into the `cpsvote` package.** +### A 2022+ Discontinuity in the `-9` "No Response" Code {#noresponse2022} + +Starting with the 2022 VRS, the Census Bureau's public-use file no longer contains any respondents coded `-9` ("No Response") for `PES1`/`VRS_VOTE`. This code is present in every release from 1994 through 2020 (2020 alone has over 10,000 such cases), and the CPS's own codebook still documents `-9` as a valid value for 2022 and 2024 -- but no respondent actually receives that code in either year's data. + +Looking at the age composition of respondents coded `-1` ("Not in Universe") is informative here: in 2016-2020, only about 12-13% of the `-1` group is 18 or older, consistent with that code overwhelmingly capturing people who were never asked the question (e.g. minors). In 2022 and 2024, that share roughly doubles to 22-23%, while the `-1` group grows by about the number of respondents who used to be coded `-9`. This strongly suggests that the Census Bureau folded "No Response" adults into the generic "Not in Universe" bucket starting in 2022, rather than that these respondents actually vanished from the survey. + +This matters because `cps_label()` treats `-1` as missing data by default (via its `na_vals` argument) *before* `cps_recode_vote()` ever runs. In years through 2020, a `-9` respondent was a real observation that got coded "No" for `cps_turnout` (per the Census's own convention) and correctly dropped to `NA` for `hurachen_turnout`. Starting in 2022, that same person -- now coded `-1` -- is wiped out to `NA` at the labelling stage and never reaches `cps_recode_vote()` at all. The practical effect is that `cps_turnout`'s denominator silently loses these respondents rather than counting them as nonvoters, which inflates the naive, uncorrected turnout estimate for 2022 and 2024 (pushing it toward the "topline"/overreported rate rather than the Census's reported rate). `hurachen_turnout` is comparatively less affected, since Don't Know/Refused/No Response respondents were already excluded from it regardless of the exact numeric code, but the composition of the state-level reweighting calibration may be subtly affected too. + +Neither `cps_label()` nor `cps_recode_vote()` will error or otherwise flag this -- they simply proceed with fewer respondents. To help catch this kind of silent coding change, `cps_recode_vote()` now calls [cps_check_codes()] automatically, which compares the codebook's documented responses for a given year against what's actually observed in the data and issues a warning when one is entirely missing (as with `-9` in 2022 and 2024). You can also call `cps_check_codes()` directly to audit a data set on its own. + +As of this writing, `cpsvote` does not attempt to algorithmically recover the lost `-9` responses (e.g. by using respondent age to split the `-1` bucket back into "Not in Universe" versus "No Response"); that would be a methodological change to the package's core turnout coding and needs further discussion before being adopted. If you are working with 2022 or later data and rely on `cps_turnout`, treat the resulting turnout estimate for those years with caution. + ## Example: Estimating Voter Turnout in 2020 {#sample2020} ```{r setup2, echo = TRUE}