doctors_and_plumbers.R

library(tidycensus)
library(tidyverse)
library(showtext)
font_add_google("Montserrat")
showtext_auto()

us_occ_pums <- get_pums(
  variables = c("OCCP", "WAGP"),
  state = "all",
  survey = "acs1",
  year = 2022,
  recode = TRUE
)

us_phys_plumb <- us_occ_pums %>%
  filter(OCCP_label %in% c("MED-Physicians", "CON-Plumbers, Pipefitters, And Steamfitters"),
         WAGP > 0)

ggplot(us_phys_plumb, aes(x = WAGP, fill = OCCP_label,
                          weight = PWGTP, color = OCCP_label)) +
  geom_histogram(alpha = 0.5, bins = 50) +
  facet_wrap(~OCCP_label, nrow = 2) +
  theme_minimal(base_size = 36, base_family = "Montserrat") +
  scale_fill_manual(values = c("navy", "darkred"),
                    guide = "none") +
  scale_color_manual(values = c("navy", "darkred"),
                    guide = "none") +
  scale_x_continuous(labels = scales::label_dollar(), breaks = c(0, 100000, 200000, 300000, 400000, 500000, 600000, 700000, 800000)) +
  scale_y_continuous(labels = scales::label_comma()) +
  labs(title = "Distributions of physician and plumber wages: 2022 1-year ACS PUMS, US Census",
       y = "Estimated count",
       x = "Wages or salary income, last 12 months",
       caption = "tidycensus R package | Individuals with no wages in 2022 removed from dataset")

Thanks a lot for your work. I highly appreciate it!

For displaying salary data, I would use boundary = 0, binwidth = 25000 instead of bins = 50. And in general, I would put percentages on the y-axis instead of counts.

Looking forward to more plots from you on LinkedIn!