plumbers.R

# -------------------------------------------------------------------------
# What do plumbers actually earn?
#
# Fact-check of the viral claim that plumbers are "clearing $250k-$300k/year",
# using the 2024 1-year American Community Survey Public Use Microdata Sample
# (PUMS), accessed with tidycensus.
#
# Universe: employed plumbers, pipefitters, and steamfitters (SOC 47-2152)
# in the United States. Estimates use the PUMS person weight (PWGTP).
# (For standard errors, add rep_weights = "person" to the get_pums() call
# and build a design with to_survey() -- slower download, 80 extra columns.)
#
# Requires a Census API key:
#   census_api_key("YOUR KEY", install = TRUE)
#   Get one at https://api.census.gov/data/key_signup.html
# -------------------------------------------------------------------------

library(tidycensus)
library(dplyr)
library(ggplot2)
library(scales)

# --- 1. Download microdata for plumbers only --------------------------------
# variables_filter subsets to SOC 47-2152 on the API side, so we only
# download plumber records rather than the full national sample.
# Verify the SOCP code with View(pums_variables)
plumbers_raw <- get_pums(
  variables = c("WAGP", "PERNP", "ADJINC", "ESR", "COW", "WKHP", "WKWN"),
  state = "all",
  survey = "acs1",
  year = 2024,
  variables_filter = list(SOCP = "472152")
)

# --- 2. Define the universe and adjust dollars ------------------------------
# ADJINC converts the rolling 12-month income reference periods to
# calendar-year 2024 dollars. PERNP is total earnings (wages + self-employment
# income), which is the fairest reading of "clearing $250k"; WAGP is wage and
# salary income alone.
plumbers <- plumbers_raw |>
  mutate(
    adj = as.numeric(ADJINC),
    adj = if_else(adj > 100, adj / 1e6, adj), # 6 implied decimals in some years
    earnings = PERNP * adj,
    wages = WAGP * adj,
    # Census definition of full-time, year-round: usually 35+ hours per week,
    # 50-52 weeks worked in the past 12 months
    full_time = WKHP >= 35 & WKWN >= 50
  ) |>
  filter(ESR %in% c("1", "2"), !is.na(earnings)) # civilians currently employed

plumbers_ft <- filter(plumbers, full_time)

# --- 3. The exact share earning more than $250,000 --------------------------
weighted_median <- function(x, w) {
  o <- order(x)
  x <- x[o]
  w <- w[o]
  x[which(cumsum(w) >= sum(w) / 2)[1]]
}

summarize_universe <- function(df) {
  summarize(
    df,
    n_sample = n(),
    n_workers = sum(PWGTP),
    median_earnings = weighted_median(earnings, PWGTP),
    share_250k_earnings = weighted.mean(earnings > 250000, PWGTP),
    share_250k_wages = weighted.mean(wages > 250000, PWGTP)
  )
}

stats_all <- summarize_universe(plumbers)
stats_ft <- summarize_universe(plumbers_ft)

print_universe <- function(label, stats) {
  cat(sprintf(
    paste0(
      "\n%s\n",
      "Unweighted sample:                  %s records\n",
      "Estimated workers:                  %s\n",
      "Median annual earnings:             %s\n",
      "Share earning > $250k (earnings):   %.2f%%\n",
      "Share earning > $250k (wages only): %.2f%%\n"
    ),
    label,
    comma(stats$n_sample),
    comma(round(stats$n_workers)),
    dollar(round(stats$median_earnings)),
    100 * stats$share_250k_earnings,
    100 * stats$share_250k_wages
  ))
}

cat("\n2024 1-year ACS PUMS, plumbers/pipefitters/steamfitters (SOC 47-2152)\n")
print_universe("All employed (includes part-time and part-year):", stats_all)
print_universe("Full-time, year-round (35+ hrs/week, 50-52 weeks):", stats_ft)
cat("\n")

share_pct <- 100 * stats_ft$share_250k_earnings

# --- 4. Histogram -----------------------------------------------------------
# Full-time, year-round workers only -- the most generous universe for the
# claim. Weighted by PWGTP so bar heights estimate counts of workers, not
# sample records. Earnings above $400k (a fraction of a percent, and
# top-coded in PUMS anyway) are shown in the last bar.

ink_primary <- "#0b0b0b"
ink_secondary <- "#52514e"
ink_muted <- "#898781"
surface <- "#fcfcfb"
gridline <- "#e1e0d9"
bar_blue <- "#2a78d6"
claim_red <- "#d03b3b"

plot_df <- plumbers_ft |>
  mutate(earnings_capped = pmin(pmax(earnings, 0), 395000))

p <- ggplot(plot_df, aes(x = earnings_capped, weight = PWGTP)) +
  geom_histogram(
    binwidth = 10000,
    boundary = 0,
    closed = "left",
    fill = bar_blue,
    color = surface,
    linewidth = 0.4
  ) +
  geom_vline(
    xintercept = 250000,
    color = claim_red,
    linetype = "42",
    linewidth = 0.6
  ) +
  annotate(
    "text",
    x = 258000,
    y = 52000,
    hjust = 0,
    vjust = 1,
    label = 'The claim:\n"clearing $250k-$300k/year"',
    color = claim_red,
    size = 3.4,
    fontface = "bold",
    lineheight = 1.1
  ) +
  annotate(
    "text",
    x = 258000,
    y = 34000,
    hjust = 0,
    vjust = 1,
    label = sprintf(
      "In reality, %.1f%% of full-time U.S.\nplumbers earn more than $250,000.",
      share_pct
    ),
    color = ink_secondary,
    size = 3.4,
    lineheight = 1.15
  ) +
  scale_x_continuous(
    labels = \(x) {
      if_else(x >= 400000, "$400k+", dollar(x, scale = 1e-3, suffix = "k"))
    },
    breaks = seq(0, 400000, 50000),
    expand = expansion(mult = c(0.01, 0.02))
  ) +
  scale_y_continuous(
    labels = label_number(scale = 1e-3, suffix = "K"),
    expand = expansion(mult = c(0, 0.05))
  ) +
  labs(
    title = "What plumbers actually earn",
    subtitle = sprintf(
      "Annual earnings (wages + self-employment income) of the %s full-time, year-round plumbers, pipefitters\n& steamfitters in the United States. Median: %s.",
      comma(round(stats_ft$n_workers, -3)),
      dollar(round(stats_ft$median_earnings, -2))
    ),
    x = "Annual earnings, 2024 dollars",
    y = "Workers per $10k bin",
    caption = "Source: 2024 1-year ACS Public Use Microdata Sample via the R tidycensus package | SOC 47-2152, civilians employed full-time (35+ hrs/week), year-round (50-52 weeks)\nPerson weights (PWGTP); earnings adjusted with ADJINC. Earnings above $400k are shown in the last bar."
  ) +
  theme_minimal(base_size = 12, base_family = "sans") +
  theme(
    plot.background = element_rect(fill = surface, color = NA),
    panel.background = element_rect(fill = surface, color = NA),
    panel.grid.major.x = element_blank(),
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_line(color = gridline, linewidth = 0.35),
    axis.ticks = element_blank(),
    axis.text = element_text(color = ink_muted, size = 9.5),
    axis.title = element_text(color = ink_secondary, size = 10),
    axis.title.x = element_text(margin = margin(t = 8)),
    axis.title.y = element_text(margin = margin(r = 8)),
    plot.title = element_text(color = ink_primary, face = "bold", size = 17),
    plot.subtitle = element_text(
      color = ink_secondary,
      size = 10.5,
      lineheight = 1.2,
      margin = margin(t = 4, b = 14)
    ),
    plot.caption = element_text(
      color = ink_muted,
      size = 7.5,
      hjust = 0,
      lineheight = 1.25,
      margin = margin(t = 12)
    ),
    plot.title.position = "plot",
    plot.caption.position = "plot",
    plot.margin = margin(16, 18, 12, 16)
  )

ggsave(
  "plumber_salaries.png",
  p,
  width = 10,
  height = 5.8,
  dpi = 200,
  bg = surface
)

message("Saved plumber_salaries.png")