Does state life expectancy correlate with political party voting?

Overview

      Do people in red states live shorter lives?

      This project examines whether state level life expectancy is statistically associated with each state's political climate in the 2024 presidential election. Political climate is measured using the popular vote margin, defined as (Trump votes - Harris votes) / (Trump votes + Harris votes). A positive value indicates a Republican advantage, and a negative value indicates a Democratic advantage.

      I'm not predicting elections, and I'm not making a claim of causality. There are many reasons why two states differ in life expectancies - differences in average income, availability of medical care, occupations with differing job hazards, etc. I'm just asking whether these two measurable state‑level quantities move together.

      A scatterplot of life expectancy versus vote margin shows a clearly downward trend: states with higher Republican vote margins tend to have lower life expectancy, and states with higher Democratic margins tend to have higher life expectancy.

      The correlation coefficient r is -0.50. This represents a moderate negative relationship - neither weak nor strong, but unmistakenly present. A significance test yields t = -4.041, p = .00019. With 51 observations, this correlation is statistically significant at the .05 level.

Details

      I obtained 2024 presidential percentage popular vote data from the Federal Election Commission. Trump won 31 of 51 (50 plus DC) states in 2024. To visualize the distribution of states, I grouped the vote margin varaible into four bins:

• Dem Majority: margin ≤ -10%
• Dem Narrow: -10% < margin ≤ 0
• Rep Narrow: 0 < margin ≤ 10%
• Rep Majority: margin > 10%

24 of Trump's 31 states were by margins greater than 10%. A bar chart and a US map show how many states fall into each group and where they are located geographically.

      The CDC (Centers for Disease Control and Prevention) publishes life expectancy tables by state. These are period life tables, showing the life expectancy of a newborn under today’s mortality rates, assuming the age‑specific death rates observed in that year (e.g., 2022) stay fixed for the newborn’s entire lifetime. Hawaii has the highest life expectancy at 80.0 years, and Mississippi has the lowest at 70.9 years.

      The following map shows life expectancies by state. I allocated the states by quartile (shortest life expectancy, shorter, longer, longest). I believe there is a relationship especially with southern states having short life expectancies in this map, compared with Republican margins in the prior map.

      The actual correlation coefficient is r = -0.50, which is moderate, but statistically significant.

      Incidentally, I had a little challenge drawing the maps with R library usmap. That library includes Puerto Rico which was not in the voter or life expectancy data, and the map would show Puerto Rico as an NA until I excluded it in the plot_usmap statement.

R code

library(readxl)
pres <- read_excel("C:/Users/Jerry/Desktop/R_files/2024presgeresults.xlsx", n_max=51)
pres$TRUMP_PERCENT <- pres$TRUMP/(pres$TRUMP + pres$HARRIS)   # ratio of popular votes
pres$HARRIS_PERCENT <- pres$HARRIS/(pres$TRUMP + pres$HARRIS)
pres$TRUMP_MARGIN <- round(pres$TRUMP_PERCENT - pres$HARRIS_PERCENT,3)
pres <- pres[, c("STATE", "TRUMP_MARGIN")]   # states are 2 letter abbrevs
colnames(pres)[1] <- "state"  # usmap requires state Column Name to be lowercase "state"

# CDC life expectancies by state:   https://www.cdc.gov/nchs/data/nvsr/nvsr74/nvsr74-12.pdf

library(pdftools)   # extract text from pdf file
library(tidyverse)
raw_text <- pdf_text("C:/Users/Jerry/Desktop/R_files/nvsr74-12.pdf")
page_text <- raw_text[3]   # page 3 only
lines <- read_lines(page_text)

# data cleaning of life expectancy file:
clean_lines <- lines %>%
  str_trim() %>%                      # Remove leading/trailing spaces
  keep(~ .x != "")                    # Drop empty rows

# Extract column headers (row 6 contains headers):
headers <- str_split(clean_lines[6], "\\s{2,}")[[1]]  # these are partial headers

# Process the data rows (Rows 7 to the end):
data_rows <- clean_lines[7:(length(clean_lines)-3)]   # delete footnotes

# Convert text lines into data frame:
life_exp <- data_rows %>%
  # Split columns whenever there are 2 or more spaces
  str_split_fixed("\\s{2,}", n = length(headers)) %>%
  as_tibble(.name_repair = "minimal")

colnames(life_exp) <- c("State", "Tot_Rank", "Tot_LE", "Tot_SE", "Male_Rank", "Male_LE", "Male_SE",
     "Fem_Rank", "Fem_LE", "Fem_SE")
life_exp$State <- gsub("\\.", "", life_exp$State)    # delete periods
life_exp$State <- sub("\\s+$", "", life_exp$State)   # delete spaces after last char 
life_exp <- subset(life_exp, State != "United States")
# convert states from names to abbreviations; District of Columbia will be NA without next line:
life_exp$State <- c(state.abb, "DC")[match(life_exp$State, c(state.name, "District of Columbia"))]   # Convert full name to 2-letter abbreviation
life_exp <- life_exp %>%
  mutate(across(where(is.character) & -1, as.numeric))   # converts all character columns in a data frame into numeric columns, except for the very first column
life_exp <- life_exp[, c("State", "Tot_LE")]   # states are 2 letter abbrevs
colnames(life_exp)[1] <- "state"  # usmap requires state Column Name to be lowercase "state"
print(life_exp)

df <- merge(pres, life_exp, by = "state")

########  Display summaries:  ########

library(ggplot2)

common_theme <- theme(
        plot.title = element_text(size=15, face="bold"),
        plot.subtitle = element_text(size=12.5, face="bold"),
        axis.title = element_text(size=15, face="bold"),
        axis.text = element_text(size=15, face="bold"),
        legend.title = element_text(size=15, face="bold"),
        legend.text = element_text(size=15, face="bold"))

df <- df %>%
  mutate(TRUMP_MARGIN_RANGE = case_when(  
    TRUMP_MARGIN <= -.10 ~ "Dem Majority",
    TRUMP_MARGIN > -.10 & TRUMP_MARGIN <= 0 ~ "Dem Narrow",
    TRUMP_MARGIN > 0 & TRUMP_MARGIN <= .10 ~ "Rep Narrow",
    TRUMP_MARGIN > .10 ~ "Rep Majority",
    TRUE ~ NA_character_ 
  )
)

percent_colors <- c("Dem Majority" = "#883068", "Dem Narrow" = "#4292C6", 
                    "Rep Narrow" = "#FB6A4A", "Rep Majority" = "#CB181D")
df$TRUMP_MARGIN_RANGE <- factor(
  df$TRUMP_MARGIN_RANGE, 
  levels = c("Dem Majority", "Dem Narrow", "Rep Narrow", "Rep Majority")
)

ggplot(df, aes(x = TRUMP_MARGIN_RANGE)) +
  geom_bar(fill = percent_colors) +
  geom_text(
    stat = "count", 
    aes(label = after_stat(count)),
    fontface = "bold", 
    vjust = -0.5
  ) +  
  labs(title="2024 Presidential Election Results by Vote Margin",
       y = "Number of States", x = "% Popular Vote Margin") +
  guides(fill = guide_legend(title = NULL)) + 
  scale_x_discrete(
    labels = c(
      "Dem Majority" = "Dem + 10% or more",
      "Dem Narrow" = "Dem 0 - 10%",
      "Rep Narrow" = "Rep 0 - 10%",
      "Rep Majority" = "Rep + 10% or more")) +
  common_theme

colSums(is.na(df))  # 0
colnames(df)[1] <- "state"   # usmap requires state Column Name to be lowercase "state"
length(df$state)   # 51
df$state <- trimws(toupper(df$state))   #51 states including DC

library(usmap)
unique(usmap::us_map(regions = "states")$full)  # Includes Puerto Rico
plot_usmap(data = df, , regions = "states", values = "TRUMP_MARGIN_RANGE", exclude = "Puerto Rico") +
  labs(title="2024 Presidential Election Results by Vote Margin") +
  scale_fill_manual(values=percent_colors,
    guide = guide_legend(title = NULL, direction="vertical"),
    labels = c(
      "Dem Majority" = "Dem + 10% or more",
      "Dem Narrow" = "Dem 0 - 10%",
      "Rep Narrow" = "Rep 0 - 10%",
      "Rep Majority" = "Rep + 10% or more")
    ) +
  theme(
    legend.position = "bottom",
    legend.box = "vertical",
    plot.title = element_text(size=15, face="bold"),
    legend.title = element_text(size=12, face="bold"),
    legend.text = element_text(size=12, face="bold")
  )

df <- df %>%
  mutate(
    LE_RANGE = case_when(
      ntile(Tot_LE, 4) == 1 ~ "Shortest LE",
      ntile(Tot_LE, 4) == 2 ~ "Shorter LE",
      ntile(Tot_LE, 4) == 3 ~ "Longer LE",
      ntile(Tot_LE, 4) == 4 ~ "Longest LE"
    )
  )

table(df$LE_RANGE)

LE_colors <- c(
  "Shortest LE" = "#D1E5F0",
  "Shorter LE"    = "#92C5DE",
  "Longer LE"     = "#4393C3",
  "Longest LE"  = "#8B3C59"
)

df$LE_RANGE <- factor(
  df$LE_RANGE, 
  levels = c("Shortest LE", "Shorter LE", "Longer LE", "Longest LE")
)

plot_usmap(data = df, regions = "states", values = "LE_RANGE", exclude = "Puerto Rico") +
  labs(title="Life Expectancies by State") +
  scale_fill_manual(values=LE_colors,
     guide = guide_legend(title = NULL, direction="vertical")) +
  theme(
    legend.position = "bottom",
    legend.box = "vertical",
    plot.title = element_text(size=15, face="bold"),
    legend.title = element_text(size=12, face="bold"),
    legend.text = element_text(size=12, face="bold")
  )


########  Corr coeff:  ########

library(ggrepel)
ggplot(data = df, mapping = aes(x = TRUMP_MARGIN, y = Tot_LE)) +
  geom_point(color = "#4D4D4D") +
  geom_smooth(method = "lm", color = "steelblue", se = FALSE, linewidth = 1) +
  geom_text_repel(
    data = df[df$state %in% c("HI", "WV"), ],
    aes(x = TRUMP_MARGIN, y = Tot_LE, label = state),
    size = 4, color = "black", fontface = "bold", nudge_y = 0.25
  ) +
  ggtitle("Life Expectancy vs % Popular Vote Margin by State") +
  xlab("Pop Vote Margin %:  Rep +, Dem - ") +
  ylab("Life Expectancy") + 
  annotate(
    "text",
    x = min(df$TRUMP_MARGIN) + .02,
    y = max(df$Tot_LE) - .02,
    label = "r = -0.50",
    hjust = 0, vjust = 1,
    color = "gray20", fontface = "bold", size = 4) +
  common_theme

r <- round(cor(df[, c("Tot_LE", "TRUMP_MARGIN")])[1, 2], 3)  # -.50 is moderate (neither weak nor strong)

# statistical significance of correlation coefficient:
n <- nrow(df)
dof <- n - 2
t <- round(r*sqrt(n - 2) / sqrt(1 - r^2),3)   # - 4.041
alpha <- .05
p_value <- round(2 * pt(abs(t), df = dof, lower.tail = FALSE), 5)    # .00019
stat_signif <- if(p_value < alpha, "is statistically significant", "is not statistically significant")
cat("r = ", r, "p-value = ", p_value, stat_signif)

End

Babe Ruth: How good was he? Mantle versus Ruth

      As a kid, my baseball hero was Mickey Mantle. In his 18-year Major League career, he amassed tremendous statistics, and surely he was one of the game's all-time greats. When Mantle retired, he was third on the all-time career home run list with 536 home runs, trailing only Babe Ruth (714) and Willie Mays (587). Unfortunately, he suffered numerous injuries, leaving fans to wonder what his numbers would look like if he had enjoyed a healthier career.

      For decades many people considered (and still consider) Babe Ruth to be the greatest all-around baseball player, because he was a star pitcher and a star hitter. The counter-argument to the Mantle "what-if" is: what if Ruth had not spent his first six years as a pitcher and instead had been a full-time hitter? (Similarly, what if Ted Williams had not spent some of his prime years in the military.)

      Nowadays Shohei Ohtani is certainly a once-in-a-generation two-way star. However, he spent his early prime years in Japan's professional baseball, and he may not have enough time left in US baseball to amass these historic lifetime stats. Aaron Judge is another exceptional modern hitter; hopefully, he stays healthy so we can watch his career numbers climb.

      To compare Mantle and Ruth’s lifetime stats objectively, I turned to data. The Sean Lahman baseball dataset contains Major League player stats back to 1871, and is available in the R library lahman. This was a good opportunity to practice data manipulation using dplyr.

      Ruth played 102 more games than Mantle over a 22-year career (2,503 versus 2,401). His lifetime statistics eclipse Mantle's in every offensive category except stolen bases — though Ruth did steal home 10 times and hit nearly twice as many triples, suggesting he was faster than most realize. His lifetime batting average was an impressive .342, sitting just behind Tris Speaker and Ted Williams.

      Lineup protection played a massive role for both men. Ruth generally batted directly before Lou Gehrig, and Mantle batted directly before Yogi Berra. While both Ruth and Mantle drew plenty of walks, presumably pitchers rarely chose to walk them intentionally just to face Gehrig or Berra.

      I was also curious about their defensive metrics. For a fair comparison, I filtered the data to isolate only their outfield appearances. There is no value in comparing an outfield throw to an assist Ruth made on a comeback grounder while pitching. Similarly, I excluded Mantle's infield appearances; he primarily played first base in his final two seasons and filled in briefly at other infield spots early in his career.

      Ruth played 222 more outfield games than Mantle (2,241 versus 2,019). Their total outfield putouts were nearly identical, likely because Mantle played centerfield and covered more ground. However, Ruth recorded nearly twice as many outfield assists (204 versus 117), which aligns with the arm strength expected of a former pitcher. While Ruth's 204 assists don't match stars such as Roberto Clemente's 266, it remains a highly respectable number.

      Ultimately, the data shows that even if you completely ignore his pitching career, Babe Ruth built an outstanding, standalone career as both a hitter and a fielder.

              
          MANTLE  RUTH
games       2401  2503
at_bats     8102  8398
runs        1677  2174
hits        2415  2873
doubles      344   506
triples       72   136
home_runs    536   714
BA         0.298 0.342
rbi         1509  2217
sb           153   123
bb          1733  2062
so          1710  1330
games_of    2019  2241
putouts     4438  4444
assists      117   204

      Here is my R code:

library(Lahman)
library(tidyverse)

data(People)
which(People$nameLast == "Mantle")   # 11773,  playerID = mantlmi01
which(People$nameLast == "Ruth")   # 16674 , playerID = ruthba01 
data(Batting)    
data(Fielding)   # for fielding, want games in outfield POS == 'OF'

df_fielding <- Fielding %>% 
  filter(playerID %in% c('mantlmi01', 'ruthba01'), POS == 'OF') %>% 
  group_by(playerID) %>% 
  summarize(games_of = sum(G, na.rm = TRUE), putouts = sum(PO, na.rm = TRUE),
    assists = sum(A, na.rm = TRUE)
  )

df <- Batting %>% 
  filter(playerID %in% c('mantlmi01', 'ruthba01')) %>% 
   group_by(playerID) %>% 
   summarize(games = sum(G, na.rm = TRUE), at_bats = sum(AB, na.rm = TRUE),
                runs = sum(R, na.rm = TRUE), hits = sum(H, na.rm = TRUE),
                doubles = sum(X2B, na.rm = TRUE), 
                triples = sum(X3B, na.rm = TRUE), home_runs = sum(HR, na.rm = TRUE),
                BA = round(hits / at_bats,3), 
                rbi = sum(RBI, na.rm = TRUE),  sb = sum(SB, na.rm = TRUE),
                bb = sum(BB, na.rm = TRUE),  so = sum(SO, na.rm = TRUE),
                )  %>%
  merge(df_fielding, by = "playerID", all.x = TRUE) 

df_transposed <- as.data.frame(t(df))
colnames(df_transposed) <- c("MANTLE", "RUTH")
df_transposed <- df_transposed[-1, ]
df_transposed

End

How much have prices increased?

      The Consumer Price Index (CPI) is a widely used measure of the prices of goods and services purchased by households. It’s the primary tool for tracking inflation and changes in the cost of living over time. The index is built from monthly price collections on a “basket” of goods and services from a sample of retail and service establishments. Historical CPI data is easy to download.

      A key feature of the CPI is that prices are adjusted for quality changes. If the price of a car’s side mirror rises by $200, but $120 of that increase reflects the mirror becoming “smart” rather than “dumb,” only the remaining $80 is counted as inflation. Similarly, if the price of a medical procedure rises because new equipment improves the quality of care, the portion attributable to improved quality is removed from the inflation calculation. Consumers still pay for these quality improvements, whether they want them or not, so in many cases the CPI understates pure cost increases.

      There are eight major categories of the CPI, and each category has its own index: Food & Beverages, Housing, Apparel, Transportation, Medical Care, Recreation, Education & Communication, and Other. These are weighted to form the overall CPI, with the largest weights as Housing at about 44% of the total, Transportation at 17%, Food & Beverages at 14%, and Medical Care at 8%. (Each of these is further sub-divided into its own index; for example, Other includes Personal Care, and Personal Care has seprate indices for Cosmetics, Perfume, Bath, and Nail Preparations.) Of course the weights will not reflect your percentages of what you buy.

      Downloading historical CPI data from FRED (Federal Reserve Bank of St. Louis) was easier than I expected. You need an API key which you can get from https://fredaccount.stlouisfed.org/login

      Here is a line graph showing cumulative CPI growth for the overall CPI and each of the eight major categories through June 2026, indexed to December 2016 = 1.00. The overall CPI has risen 37.1% since December 2016. Housing (rent, insurance, energy, etc.) has increased the most at 44.8%. Transportation (vehicle purchases, fuel, maintenance, insurance, public transit fares) is next at 43.5%. Food (groceries and restaurants) is up 39.9%. Medical Care is lower at 25.9%. Medical Care includes out‑of‑pocket spending on providers, hospitals, and insurance, but excludes employer‑paid and government‑paid health insurance premiums.

      The CPIs exclude a lot of things like the quality changes I mentioned above, and other items that you may pay but that the government does not classify as personal expenses. We all feel the cost of our groceries going up - I discussed this previously in groceries . Each category has its own reasons why it is increasing; I will leave that discussion to the economists, except to say that the cost of energy affects a lot of items in the cost of production and delivery.

      Like many broad measures, the CPI is an attempt to estimate the overall cost of goods and services. But what ultimately matters to you is the actual cost of the things you buy.

      Here is my R code:

library(tidyquant)
library(dplyr)
library(tidyr)
library(lubridate)
library(ggplot2)
library(ggrepel)   # repel overlapping text labels

# Set your API environment variable
Sys.setenv(FRED_API_KEY = "xxxx")


# Define the official FRED database tracking codes
cpi_series <- c(
  "CPIAUCSL", "CPIFABSL", "CPIHOSSL", "CPIAPPSL", 
  "CPITRNSL", "CPIMEDSL", "CPIRECNS", "CPIEDUNS", "CPIOGSNS"
)

# Download and clean data vectors
raw_data <- tq_get(cpi_series, get = "economic.data")   # get from FRED

# 2. Clean, Filter, and Perform Group-Indexing
cpi_processed <- raw_data %>%
  mutate(
    Year  = year(date),
    Month = month(date)
  ) %>%
  # Keep all Decembers from 2015 onward OR strictly isolate June 2026
  filter((Month == 12 & Year >= 2015) | (Year == 2026 & Month == 6)) %>%
  # Convert raw tracking codes into readable titles
  mutate(Category = case_when(
    symbol == "CPIAUCSL" ~ "Overall CPI",
    symbol == "CPIFABSL" ~ "1. Food & Bev",
    symbol == "CPIHOSSL" ~ "2. Housing",
    symbol == "CPIAPPSL" ~ "3. Apparel",
    symbol == "CPITRNSL" ~ "4. Transportation",
    symbol == "CPIMEDSL" ~ "5. Medical Care",
    symbol == "CPIRECNS" ~ "6. Recreation",
    symbol == "CPIEDUNS" ~ "7. Education & Comm",
    symbol == "CPIOGSNS" ~ "8. Other Goods"
  )) %>%
  # Chronologically sort each group, then anchor base-100 to the first row (Dec 2015)
  group_by(Category) %>%
  arrange(date, .by_group = TRUE) %>%
  mutate(Indexed_Value = (price / first(price)) * 100) %>% 
  ungroup() %>%
  # Convert Timeline to ordered categories for a clean discrete X-Axis
  mutate(Period = if_else(Month == 6, paste0(Year, " (June)"), as.character(Year))) %>%
  mutate(Period = factor(Period, levels = unique(Period[order(date)])))

# 3. Isolate final data point rows for the text tags
label_data <- cpi_processed %>%
  group_by(Category) %>%
  filter(date == max(date)) %>%
  ungroup()

# 4. Generate the Chart with the Categorical String Axis Baseline
ggplot(cpi_processed, aes(x = Period, y = Indexed_Value, color = Category, group = Category)) +
  geom_line(aes(linewidth = ifelse(Category == "Overall CPI", 1.5, 0.8))) +
  geom_point(size = 2) +
  
  # Non-overlapping direct text labels
  geom_text_repel(
    data = label_data,
    aes(label = paste0(Category, " (", round(Indexed_Value, 1), ")")),
    nudge_x = 0.5,             
    direction = "y",          
    hjust = 0,                
    segment.color = "grey50", 
    segment.size = 0.4,
    force = 2,
    fontface = "bold",
    size = 4   # millimters
  ) +
  
  # High-contrast visual color mapping matrix
  scale_color_manual(values = c(
    "Overall CPI"         = "#000000", 
    "1. Food & Bev"       = "#E64B35", 
    "2. Housing"          = "#56B4E9", 
    "3. Apparel"          = "#009E73", 
    "4. Transportation"   = "#4D8805", 
    "5. Medical Care"     = "#0072B2", 
    "6. Recreation"       = "#D55E00", 
    "7. Education & Comm" = "#CC79A7", 
    "8. Other Goods"      = "#999999"  
  )) +
  
  # Format plot margins to prevent label cropping
  scale_x_discrete(expand = expansion(mult = c(0.05, 0.35))) +
  scale_linewidth_identity() + 
  labs(
    title = "10.5-Year Cumulative CPI Growth Comparison",
    subtitle = "Base Index: December 2016 = 100",
    x = "Reporting Period",
    y = "Indexed Value (Relative to 100)"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    legend.position = "none",          
    panel.grid.minor = element_blank(),
    text = element_text(face = "bold"),
    plot.title = element_text(face = "bold", size = 14),
    axis.text = element_text(face = "bold")
  )

# 5. Save the final graphic output file
# ggsave("cpi_growth_comparison.png", width = 12, height = 7, dpi = 300, bg = "white")

# Extract and print the final 10.5-year cumulative values
final_column_summary <- cpi_processed %>%
  filter(date == max(date)) %>%
  select(Category, Indexed_Value) %>%
  mutate(Indexed_Value = round(Indexed_Value, 2)) %>%
  arrange(desc(Indexed_Value)) # Sorts from highest inflation to lowest

print(as.data.frame(final_column_summary))

End

Testing racial predictions with BISG

      Regulators expect insurance companies, banks, and others not to discriminate by race, but businesses are not allowed to collect the race of their customers. One solution for testing discrimination in the aggregate is to use U.S. Census data to estimate the probable race of each customer based on the customer’s name and geographic location. (The Census race data is self identified and excludes a multi racial category; I will leave it to the sociologists to discuss issues with these things and whether Hispanic is a race.)

      The common algorithm for estimating race from names and geography is BISG (Bayesian Improved Surname Geocoding). As an example, suppose there is a customer Mary Johnson who lives in Essex County, NJ. If we want to predict her race only from her surname, we use nationwide Census probabilities such as P(Black ∣ Johnson) = 0.3441, P(White ∣ Johnson) = 0.5438, etc.

      To improve the prediction, we add her county, using geographic weights such as the relative density of each race in Essex compared with nationwide. The resulting probability P(Black | Johnson & Essex) = (Black national prob × Black geographic weight) / ∑_race (race national prob × race geographic weight) = 0.7393.

      Finally, we add her first name. Under a Naive Bayes assumption that surname, location, and first name are conditionally independent given race, Bayes’ rule allows us to multiply the components: P(race ∣ surname & geo & first ) ≈ P(race) x P(surname ∣ race) x P(geo ∣ race) x P(first ∣ race). The probability P(Black ∣ Johnson & Mary & Essex) becomes 0.7158.

      A summary of these probabilities is:

Model Type	White	Black	Hispanic	Asian	Other
Surname Only	.5438	.3441	.0272	.0074	.0774
Surname+Geo	.1766	.7393	.0244	.0046	.0551
Surname+First+Geo	.2481	.7158	.0042	.0013	.0306

      These calculations can be done with the R package wru (Who are you?). You will need a Census API key to download the Census data, available from https://api.census.gov/data/key_ssignup.html .

      I was interested in measuring the accuracy of the BISG predictions. It is not easy to find real (not simulated) data linking name and race. However, several states make their voter registration records public. I requested and downloaded Florida voting registration data from Harvard Dataverse ( https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/UBIG3F ). I limited my study to a single Florida county, Lee, which had over 500,000 registered voters.

      The Florida voting data required considerable data cleaning. Among the amusing data issues is that some Lee County voters registered with a specialty post office box zip code, or a zip code with a typo, or an out-of-state zip code. The first two were mapped to county level, and the third category was dropped.

      There were 545,187 registered Lee County voters, distributed by race as follows:

Asian	Black	Hispanic	Other	White	NA
6912	32407	66677	9435	422333	7423

      I performed the same BISG algorithm as above, against the Florida voter registration data for Lee County. The overall accuracy is 88%, but this is misleading because the Florida data is highly imbalanced toward the White race. A better accuracy measure is one that is calculated separately for each race, and I chose Balanced Accuracy = (Sensitivity + Specificity)/2. Sensitivity measures how well the algorithm finds people of that specific race, and specificity measures how well the algorithm avoids misclassifying people of other races into that specific race.

      My Balanced Accuracy results of comparing BISG race predictions versus Florida self-identified race data is as follows:

White	Black	Hispanic	Asian	Other
.8226	.7179	.9007	.7251	.5041

      I think these results are pretty good for Hispanic and White, not as good for Black and Asian, and of course pretty poor for Other. Other includes Native American, Native Hawaiian, and others, and are a small overall percentage of the county population.

The citation for the Florida voter data is: Sood, Gaurav, 2017, "Florida Voter Registration Data (2017 and 2022)", https://doi.org/10.7910/DVN/UBIG3F, Harvard Dataverse, V2.

      Here is my R code to predict Mary Johnson.

library(wru)
library(dplyr)

# 1. Create a 1-row dataset for Mary Johnson in Essex County, NJ (FIPS "013")
test_voter <- data.frame(
  surname = "JOHNSON",
  first = "MARY",
  state = "NJ",
  county = "013", # Essex County FIPS
  stringsAsFactors = FALSE
)

# 2. Grab the live Census geographic data for New Jersey
nj_census <- get_census_data(
  key = "xxxx",
  states = "NJ",
  age = FALSE,
  sex = FALSE
)

# 3. Generate the 3 prediction variations
pred_baseline  <- predict_race(voter.file = test_voter, surname.only = TRUE)
pred_geo       <- predict_race(voter.file = test_voter, census.data = nj_census, census.geo = "county")
pred_first_geo <- predict_race(voter.file = test_voter, census.data = nj_census, census.geo = "county", names.to.use = "surname, first")

# 4. Format and clean output to perfectly fit the RStudio layout
combined_predictions <- bind_rows(
  pred_baseline  %>% mutate(model_type = "Surname Only"),
  pred_geo       %>% mutate(model_type = "Surname+Geo"),
  pred_first_geo %>% mutate(model_type = "Surname+First+Geo")
) %>%
  rename(
    White = pred.whi,
    Black = pred.bla,
    Hispan = pred.his,
    Asian = pred.asi,
    Other = pred.oth
  ) %>%
  mutate(across(c(White, Black, Hispan, Asian, Other), ~ round(., 4))) %>%
  select(model_type, surname, first, White, Black, Hispan, Asian, Other)

# 5. Display the output summary
print(combined_predictions, row.names = FALSE)
# End of Mary Johnson prediction
##################################################################################

      Here is my R code to measure the accuracy of the BISG predictions against Lee County voters. A considerable portion of the code is for data cleaning of the Florida voter registration data.

library(data.table)
library(dplyr)
library(stringr)
library(caret)
library(wru)

voter_file_path <- "xxxx/LEE_20220621.txt"
Sys.setenv(CENSUS_API_KEY = "xxxx")

raw_lines <- fread(
  file = voter_file_path,
  sep = NULL, # Read the row as a single string to lock layout shifting
  header = FALSE,
  col.names = "raw_row"
)

wru_ready_data <- raw_lines %>%
  mutate(
    # Extract structural components securely out of text blocks
    surname = str_match(raw_row, "\\d{9}\\s+(\\w+)")[,2],
    first = str_match(raw_row, "\\d{9}\\s+\\w+\\s+(\\w+)")[,2],
    zipcode = str_extract(raw_row, "\\b3\\d{4}\\b"),
    raw_race = str_match(raw_row, "\\b([MF])\\s+(\\d)\\b")[,3]
  ) %>%
  # Filter out empty records or text extraction errors
  filter(!is.na(surname) & !is.na(raw_race)) %>%
  # Map Florida's state administrative codes to wru categories
  mutate(
    true_race = case_when(
      raw_race == "5" ~ "white",
      raw_race == "3" ~ "black",
      raw_race == "4" ~ "hispanic",
      raw_race == "2" ~ "asian",
      raw_race %in% c("1", "6", "7") ~ "other",
      TRUE ~ NA_character_ # Automatically standardizes missing responses to NA
    )
  )

# Clean data of invalid zip codes and reroute PO Boxes to residential ZCTAs
wru_ready_data <- wru_ready_data %>%
  # Exclude out-of-state ZIP codes (Georgia 30xxx/31xxx, Alabama 35xxx/36xxx)
  filter(!substr(zipcode, 1, 2) %in% c("30", "31", "35", "36")) %>%
  # Map missing/typos and PO boxes directly to matching valid residential ZCTAs
  mutate(
    zipcode = case_when(
      # Defunct / Typos / Missing data map to Central Lee County Baseline
      zipcode %in% c("fl-NA", "33929", "33932", "33945", "33970") ~ "33901",
     
      # PO Box explicit routing to residential census counterparts
      zipcode == "33902" ~ "33901",  # Fort Myers PO Box -> Fort Myers Residential
      zipcode == "33906" ~ "33907",  # Fort Myers PO Box -> South Fort Myers Residential
      zipcode == "33910" ~ "33904",  # Cape Coral PO Box -> Cape Coral Residential
      zipcode == "33915" ~ "33919",  # Fort Myers PO Box -> Cypress Lake Residential
      zipcode == "33918" ~ "33903",  # N. Fort Myers PO Box -> N. Fort Myers Residential
      zipcode == "33994" ~ "33928",  # Bonita Springs PO Box -> Estero/Bonita Residential
      zipcode == "34106" ~ "34102",  # Naples PO Box -> Naples Residential
      zipcode == "34133" ~ "34135",  # Bonita Springs PO Box -> Bonita Residential
      zipcode == "34136" ~ "34135",  # Bonita Springs PO Box -> Bonita Residential
     
      TRUE ~ zipcode # Keep all other valid residential ZCTAs as they are
    ),
    # Ensure wru recognizes the geographic county boundary for the fallback imputation
    county = "12071" # FIPS code for Lee County, FL
  )

nrow(wru_ready_data)
table(wru_ready_data$true_race, useNA = "always")

lee_county_test <- as.data.table(wru_ready_data) %>%
  filter(!is.na(true_race)) %>%
  mutate(
    state = "fl",
    surname = as.character(surname),
    first = as.character(first),
    zcta = as.character(zipcode) # Required column label mapping for ZCTA use
  )

florida_zcta_layers <- wru::get_census_data(
  key = Sys.getenv("CENSUS_API_KEY"),
  state = "FL",
  age = FALSE,
  sex = FALSE,
  census.geo = "zcta",
  county.list = NULL # Required parameter boundary choice for independent ZCTAs
)

predicted_lee_county <- wru::predict_race(
  voter.file = lee_county_test,
  census.surname = TRUE,
  surname.only = FALSE,
  census.geo = "zcta",
  census.data = florida_zcta_layers,
  impute.missing = TRUE, # Forces wru to fall back to County priors if needed
  skip_bad_geos = TRUE
)

# "53233 (10.1%) individuals' last names were not matched, but ZCTA's baseline geographic racial demographics to calculate the prediction.

evaluation_ready_data <- predicted_lee_county %>%
  rowwise() %>%
  mutate(
    # Isolate which of the five columns held the highest posterior probability string
    max_prob_col = c("white", "black", "hispanic", "asian", "other")[
      which.max(c(pred.whi, pred.bla, pred.his, pred.asi, pred.oth))
    ]
  ) %>%
  ungroup() %>%
  # Convert fields into factor arrays to prevent tracking sequence errors
  mutate(
    true_race = factor(true_race, levels = c("white", "black", "hispanic", "asian", "other")),
    predicted_map = factor(max_prob_col, levels = c("white", "black", "hispanic", "asian", "other"))
  )

# Compute the final metric confusion matrix tracking loop
accuracy_report <- confusionMatrix(
  data = evaluation_ready_data$predicted_map,
  reference = evaluation_ready_data$true_race
)

print(accuracy_report$overall["Accuracy"]) # Overall model precision
print(accuracy_report$byClass[, "Balanced Accuracy"]) # Success broken down per group

End

RSA Encryption

      When I was in high school, students would pass notes to each other in class. I even had an English teacher who was rumored to be having an affair with another teacher, and she would pass notes to him, using a student to deliver the note. Nowadays, students would use their phones and text each other.

      I recently learned a small amount about encryption for a discrete math class. Encryption is the process of scrambling readable data into unreadable code using a secret mathematical key. Cryptography is the general subject of securing information by converting it into unreadable code. During World War II, mathematician Alan Turing played a crucial role in deciphering German messages. Also during World War II, Native American Code Talkers created unbreakable code that enemy cryptographers could never decipher. Today encryption is used in protecting usernames, passwords, and personal data from interception, in many areas such as financial transactions.

      Cryptography is a vast subject. I will limit this essay to the RSA algorithm which is reasonably simple, although it becomes computationally complex for large files. I am also intentionally using small numbers; real RSA uses enormous numbers that make it hard to hack. In fact, RSA is being phased out starting in 2030 due to the threat of future quantum computers.

      Briefly, if you and I agree to use the RSA encryption algorithm, and you have the public key number (n,e) where n is the product of two prime numbers and you have the private key number d, then I can give you an encrypted message and you can decrypt it.

      Suppose the public key is (n,e) = (33,3) and the private key is d = 7. I have a message, and I encrypt it. Suppose my encrypted message is 3 32 12 3 11 26 32 16 9 21 15 32 28 1 5 31 24 1 32 8 21 24 5 28 . Prior to encryption, 1 is a, 2 is b, etc., 27 is comma, and 32 is space. If you can decrypt this, you will know a secret of my youth :)

      Recall the modulo operator "mod": if x and y are integers, x mod y is the remainder after x/y. For example, 17 mod 12 = 5, because 17/12 = 1 + remainder 5.

      The following is a brief summary of RSA encryption with small numbers:

The recipient has a public key (n,e) that is available to the sender, eg (33,3); and a private key, eg d = 7. (The recipient first chooses primes p and q, eg 3 and 11, and computes n = p*q. The recipient computes Totient φ = (p-1)(q-1)=20, and selects an e that is coprime to φ. Then the recipient performs a modular inverse to solve for d in e*d = 1 mod φ, eg 3*d = 1 mod 20; this can be done using the Extended Euclidean Algorithm, although R package numbers does this as d <- modinv(e,φ) . Here d is 7.)

The sender converts the message from letter to plaintext m, eg m=13; and then from m to encrypted ciphertext = (m^e)(mod n), eg c = (13^3)(mod 33) = 19. (The R expression is (m^e) %% n . However, this may fail if numbers get too large. One alternative is modpower(13,3,33) from the numbers package.)

The recipient decrypts c to recover m, m = (c^d)mod n, eg m = (19^7)(mod 33) = 13.

As an example, if the message is “math”, the plaintext m is (13 1 20 8), the encrypted ciphertext c is (19 1 14 17), the decrypted plaintext m is (13 1 20 8), and the decrypted message is (m a t h).

      Here is my R code.

library(numbers)

# --- RECIPIENT KEYS ---
p <- 3
q <- 11
n <- p * q
phi <- (p - 1) * (q - 1)
e <- 3
d <- modinv(e, phi)

# --- SENDER SIDE ---
m_text <- "MATH"
m_text <- tolower(m_text)
print(paste("The unencrypted message m is:", m_text))

char_vector <- unlist(strsplit(m_text, split = ""))

# Mapping: letters → 1–26, punctuation → custom codes
special_lookup <- c("," = 27, "." = 28, " " = 32)

m_tokens <- sapply(char_vector, function(char) {
  if (char %in% names(special_lookup)) {
    special_lookup[char]
  } else if (char %in% letters) {
    match(char, letters)
  } else {
    stop(paste("Unsupported character:", char))
  }
})

# Encryption
c <- sapply(m_tokens, function(msg_num) modpower(msg_num, e, n))
print(paste("The encrypted ciphertext c is:", paste(c, collapse = " ")))

# Decryption
m_decrypted <- sapply(c, function(cipher_num) modpower(cipher_num, d, n))

# Reverse mapping
result <- sapply(m_decrypted, function(num) {
  if (num %in% special_lookup) {
    names(special_lookup)[special_lookup == num]
  } else {
    letters[num]
  }
})

print(paste("The decrypted ciphertext integers are:", paste(m_decrypted, collapse = " ")))
print(paste("The decrypted message is:", paste(result, collapse = " ")))

End

Transgender Day of Visibility

      March 31 is Transgender Day of Visibility. I’m not transgender myself, but I have friends, acquaintances, and family members who are. Chances are you do too, whether you realize it or not.

      My understanding is that transgender is an umbrella term for people whose gender identity or expression differs from the sex they were assigned at birth. It reflects a deeply held internal sense of self—something not defined by appearance, clothing, or medical procedures. Being transgender is about identity, not sexual orientation.

      What I’ve learned over time is that many transgender people face challenges most of us never have to think about. These include discrimination, gaps in legal protection, denial of medical care, and even physical violence. There are also everyday barriers that rarely make headlines, like difficulty obtaining accurate driver’s licenses or passports—documents most of us take for granted.

      I found this resource helpful: Understanding the Transgender Community .

      I wish every transgender person could live their life openly, safely, and without being hassled for who they are.

      The five color Transgender Pride Flag was designed by Monica Helms in 1999. I made this flag in computer language R. Here is the R code.

########################

library(ggplot2)

# Define the colors in order: Blue, Pink, White, Pink, Blue
trans_colors <- c("#5BCEFA", "#F5A9B8", "#FFFFFF", "#F5A9B8", "#5BCEFA")

# Create a data frame for the 5 stripes
flag_data <- data.frame(
  stripe = factor(1:5),
  height = rep(1, 5)
)

# Plot the flag
ggplot(flag_data, aes(x = 1, y = height, fill = stripe)) +
  geom_bar(stat = "identity", width = 1, color = NA) +
  scale_fill_manual(values = rev(trans_colors)) + # Reverse to stack correctly
  theme_void() + # Remove axes and labels
  theme(legend.position = "none") +
  coord_cartesian(expand = FALSE)

End

Does every finite string of numbers appear within π ?

      In honor of Pi Day (March 14), I offer the following: Does every finite string of numbers like your social security number eventually appear somewhere within the decimal expansion of π ?

      So said Harold Finch, the computer genius on TV show "Person of Interest" ("You are being watched ... ") in a 2013 episode where he poses as a substitute high school math teacher.

      "Pi... keeps on going, forever, without ever repeating. Which means that contained within this string of decimals, is every single other number. Your birthdate, combination to your locker, your social security number, it's all in there, somewhere. And if you convert these decimals into letters, you would have every word that ever existed... all of the world's infinite possibilities rest within this one simple circle."

      We know π is irrational; it cannot be expressed as the quotient of two integers.   Video

      An irrational number has an infinite, non-repeating decimal expansion.

      A number is defined as normal if any finite string of numbers like your social security number will eventually appear somewhere in its decimal expansion. (The is a simpler version of a more complicated mathematical definition , but it is sufficient for this discussion.) Mr. Finch is claiming π is a normal number. However, mathematicians have not yet proved this. (I once attempted to add this last statement and the fact that Finch was wrong, to Wikipedia's Person of Interest page, but the Wikipedia editors declined it.)

      Although π is defined geometrically as the ratio of a circle's circumference to its diameter, mathematicians calculate long sequences of π using infinite series, rather than physical measurements. Various websites such as angio provide files containing up to 1 million digits of π.

      I wrote some R code that would take a number (or a word, converted to a number) and search for where that number could be found within π. I decided the first 100,000 digits of π were enough for me. Of course, the fact that a number can not be found within the first 100,000 digits, does not preclude that it can be found in the next 100,000 digits, or the next 1 million digits. And so on. Infinitely many digits goes on for a long time.

      I tested these examples:

• "Eggs". e = 5, g = 7, g = 7, s = 19.   57719 appears starting at the 6026 th digit of π.
• "0123" appears starting at the 27846 th digit of π.
• "01234567" was not found in the first 100k digits. But according to The Pi Search Page it can be found starting at the 112,099,767 digit.

I did not try "0123456789", but another source says it was not found in the first 2 billion digits.

      Here is my R code. Happy Pi Day!

# Symbol for pi is \u03c0
library(stringr)
library(readr)

find_in_pi <- function(input_value) {
  # 1. Load the Pi data
  data.raw <- readr::read_file("https://assets.angio.net/100000.txt")
  data.vec <- unlist(str_split(data.raw, pattern = ""))
  data.vec <- data.vec[-c(1, 2)] # remove '3.'
  pi_string <- paste(as.character(data.vec), collapse = "")
  
  # 2. Process Input: Determine if it's a word or a numeric string
  # Convert to character and lowercase for uniformity
  clean_input <- tolower(as.character(input_value))
  
  if (grepl("[a-z]", clean_input)) {
    # It's a word: Convert letters to numbers (1-26)
    word_only <- gsub("[^a-z]", "", clean_input)
    z <- unlist(str_split(word_only, ""))
    q <- match(z, letters)
    search_target <- paste(q, collapse = "")
  } else {
    # It's already a number: Just strip non-digits (like decimals or spaces)
    search_target <- gsub("[^0-9]", "", clean_input)
  }
  
  # 3. Search Pi
  pos <- regexpr(search_target, pi_string)
  
  # 4. Handle Results
  if (pos[1] == -1) {
    return(paste0("The sequence '", search_target, "' was not found in the first 100k digits."))
  } else {
    start <- pos[1]
    end <- start + nchar(search_target) - 1
    return(list(
      found = search_target,
      start_index = start,
      end_index = end,
      index_sequence = seq(from = start, to = end)
    ))
  }
}

# Search for a word
find_in_pi("Eggs")

# Search for a number;  enclose in quotes if a leading zero
find_in_pi("01234")
find_in_pi("01234567")

End

Years having a month with exactly four weeks

      A year having a month with exactly four weeks will occur in 2026, 2037, 2043, 2054.

      Such a year has February 1 on a Sunday, and is a non-Leap year so February has exactly 28 days.

      Here is the R code:

library(lubridate)

date_sequence <- seq(
  as.Date("2026-02-01"),
  as.Date("2056-02-01"),
  by = "years"
)

df <- data.frame(
  Feb_1 = date_sequence,
  day_of_week = weekdays(date_sequence), 
  is_leap = leap_year(date_sequence)
)

df <- subset(df, day_of_week == "Sunday" & is_leap == FALSE)

cat("The following years have a month (Feb) with exactly four weeks:", (year(df$Feb_1)))

# 2026, 2037, 2043, 2054 

End

My other girlfriend, Irmaa. (Not a typo: two a's)

     

      She is lovely, isn't she?

      Irmaa sat alone at a corner table in the kind of restaurant where the silverware had weight and the waiters moved like well trained ghosts. A single glass of Sancerre caught the light, throwing a pale shimmer across the white tablecloth. She liked that — the quiet glow, the sense that even photons paused to admire her.

      Her simple but elegant hairstyle, with a few intentional strands of silver, framed a face that could shift from stern to playful in a heartbeat. Tonight it hovered somewhere in between. She tapped a manicured finger against the stem of her glass, scanning the room with the calm confidence of someone who knew she brought in more revenue than most small nations.

      The maître d’ approached with a deferential smile. “Will anyone be joining you this evening, Ms. Irmaa?”

      She returned the smile, subtle and knowing. “Oh, people join me every year,” she said, “but never at the same table.”

So who is Irmaa?

      My real significant other girlfriend knows one of my life philosophies is "You can only have one Honey." When I am watching TV or a movie, and a character who already has a partner starts looking for a second partner, I predict bad things will happen. Male or female, straight or gay, it doesn't matter. "You can only have one Honey."

      And true to my philosophy, I have, reluctantly, added Irmaa to my life, and bad things have happened. Sigh.

      IRMAA stands for Income-Related Monthly Adjustment Amount. It is a surcharge directly deducted from your Social Security check, and added to Medicare Part B and Part D premiums for beneficiaries with higher incomes. IRMAA generates tens of billions of dollars each year. IRMAA premiums are slated to rise 30 percent from 2026 to 2030, according to the 2025 Medicare Trustees Report. In comparison the Social Security cost of living adjustment was set at 2.8 percent for 2026. See IRMAA.

      My IRMAA surcharge is quite large relative to my standard Medicare premium.

      If you imagine IRMAA as a single individual collecting all those surcharges herself, she’d be pulling in revenue on the scale of tech billionaires. IRMAA as a person would be staggeringly wealthy.

      I think of her as my girlfriend named Irmaa. I am one of those who pays for her extravagant lifestyle. She takes and takes, and does not give much back. Of course, she has quite a few other boyfriends. And girlfriends.

      Irma (single a) is a quite uncommon name nowadays. Irma is currently #2609 in U.S. births according to The Bump's Data. The World Meteorological Organization retired Irma as an Atlantic hurricane name after its disasterous 2017 hurricane, feeling the future use of the name would be insensitive.

.       Meanwhile Irma calmly enjoys her expensive meals, confidently knowing over 5 million Medicare beneficiaries - men and women - support her lifestyle. She's not your girlfriend yet? Just wait a few years.

      (Photo from Google Gemini, description of Irmaa by me in consultation with Microsoft Copilot.)

Do prices feel like they are rising more than the CPI says?

      Do prices feel like they are rising more than the CPI (Consumer Price Index) says? Yes.

      For 2025 the overall CPI increased by a pretty moderate 2.7% (CPI-U for all urban consumers, December to December). The CPI measures the price change of a basket of goods over eight major categories: food & beverages, housing, apparel, transportation, medical care, recreation, education & communication, and other, over a large number of metropolitan statitical areas (MSA). Each month Bureau of Labor Statistics data collectors collect about 94,000 prices. For some items, prices are adjusted to remove the effect of a change of quality (for example today's Toyota Corolla is a much different car than ten years ago). There are separate indices for over 200 items countrywide, and separate indices for the eight major categories by MSA.

      Of course what matters is price changes on the items YOU buy, and especially the items you buy frequently like food. Each individual item has its own reasons for its price changes - egg prices are affected by bird flu, coffee prices are affected by weather and geopolitical conditions in other countries, etc. These factors are beyond our control.

      I was interested in comparing food price changes for common frequently purchased items. This is based on CPIs for countrywide data in categories like eggs. Of course I am not measuring YOUR specific egg purchase (such as store brand, organic, medium size, white, 12 count, at CTown in Intercourse, Pennsylvania).

      I began with 11 years of December values of ten separate food categories and also the overall CPI-U, seasonally adjusted. The values I used are as follows (although I understand values change as there are changes to seasonal adjustment factors, weightings, etc.)

      For comparison purposes I set all eleven indices to a 100 index value at 2015, and I plotted them on a line graph. Each year's value is now relative to 2015.

      On a ten-year basis, ground beef, eggs, and lettuce are running higher than the overall CPI. Those three food indices have experienced a lot of volatility. The pandemic year 2022 was the year of the largest increase for many foods. Interestingly, apple prices have gone down!

      Another way to look at this is to look at bar graphs with percentage changes. Here are graphs with percentage changes 2015 to 2025, 2024 to 2025, and 2021 to 2022:

      For 2025 compared with 2015, I was not aware how much the prices had risen for ground beef, lettuce, and eggs.

      For 2025 compared to 2024, the increase in prices for ground beef, eggs, and coffee have been much larger than the 2.7% overall increase in the CPI. This is probably no surprise if you are the household member doing the grocery shopping.

      For 2022 versus 2021, the increase in all foods except ground beef exceeded the overall CPI. We probably remember this during the Covid pandemic years.

      In the insurance world, the goods and services that insurance pays for have also increased by more than the overall CPI. For example, medical care costs, auto maintenance & repair costs, and building construction & repair costs have consistently outpaced the overall CPI. So don't be surprised to see health insurance, auto insurance, and home insurance prices to continue to increase.

      Here is the R code:

library(readxl)
library(dplyr)
library(tidyr)
library(ggplot2)

df <- read_excel("C:/Users/Jerry/Desktop/R_files/CPI_10_years.xlsx")
print(df, n = Inf, width = Inf)

# 1. Convert data to 'long' format for easier plotting
df_long <- df %>%
  pivot_longer(cols = -Year, names_to = "Category", values_to = "Index")

# 2. Re-base every category so 2015 = 100
df_normalized <- df_long %>%
  group_by(Category) %>%
  mutate(Index_100 = (Index / Index[Year == 2015]) * 100) %>%
  ungroup()

# 3. Plot line graph "Race to the Top"
ggplot(df_normalized, aes(x = Year, y = Index_100, color = Category)) +
  geom_line(aes(size = Category == "Overall CPI"), show.legend = FALSE) +
  # Use subset to only label the last point (2025)
  geom_text(data = subset(df_normalized, Year == 2025), 
            aes(label = Category), 
            hjust = -0.1,  # Push text to the right of the point
            size = 4, fontface="bold",
            show.legend = FALSE) +
  # Fix X-Axis: No decimals, every year labeled
  scale_x_continuous(breaks = seq(2015, 2025, 1), limits = c(2015, 2027)) +
  # Highlight the Overall CPI in black/thicker line
  scale_size_manual(values = c("TRUE" = 2.0, "FALSE" = 0.8)) +
  scale_color_manual(values = c("Overall CPI" = "black", 
                                "Eggs" = "red", 
                                "Ground Beef" = "darkred",
                                "Apples" = "forestgreen",
                                "Milk" = "blue",
                                "Bread" = "orange",
                                "Coffee" = "brown",
                                "Lettuce" = "lightgreen",
                                "Soup" = "#008080",
                                "Breakfast Cereal" = "purple",
                                "Other Fresh Veg" = "gray")) +
  theme_minimal() +
  theme(
    legend.position = "none",        # Remove legend
    plot.title = element_text(face = "bold", size = 16), 
    axis.title = element_text(face = "bold"), 
    axis.text = element_text(face = "bold"),
    plot.margin = margin(r = 100),   # Add space on the right for labels
    panel.grid.minor = element_blank()
  ) +
  coord_cartesian(clip = 'off') +    # Prevent labels from being cut off
  labs(title = "Cumulative Food Price Changes (2015 = 100)",
       subtitle = "Food categories vs. Overall CPI Baseline",
       y = "Index (2015 = 100)",
       x = "Year")

#4 Plot bar graph with % changes
# filter data for start and end years 

plot_cpi_change <- function(data, start_yr, end_yr) {
  
  # 1. Prepare dynamic column names for the mutation step
  start_col <- paste0("Yr", start_yr)
  end_col <- paste0("Yr", end_yr)
  
  # 2. Data Processing
  df_bar <- data %>% 
    filter(Year %in% c(start_yr, end_yr)) %>%
    pivot_longer(cols = -Year, names_to = "Category", values_to = "Index") %>% 
    pivot_wider(names_from = Year, names_prefix = "Yr", values_from = Index) %>%
    # Use .data[[]] to dynamically reference the columns created by pivot_wider
    mutate(Pct_Change = ((.data[[end_col]] - .data[[start_col]]) / .data[[start_col]]) * 100) %>%
    mutate(Category = reorder(Category, -Pct_Change))
  
  # 3. Visualization
  ggplot(df_bar, aes(x = Category, y = Pct_Change, fill = Category)) +
    geom_bar(stat = "identity", color = "black", size = 0.5) +
    # Logic for text placement based on positive/negative change
    geom_text(aes(label = paste0(round(Pct_Change, 1), "%"),
                  vjust = ifelse(Pct_Change >= 0, -0.5, 1.5)), 
              fontface = "bold", size = 4.5) +
    scale_fill_manual(values = c("Overall CPI" = "black", "Soup" = "#008080", "Eggs" = "red", 
                                 "Ground Beef" = "darkred", "Apples" = "forestgreen", "Milk" = "blue", 
                                 "Bread" = "orange", "Coffee" = "brown", "Lettuce" = "lightgreen", 
                                 "Breakfast Cereal" = "purple", "Other Fresh Veg" = "gray")) +
    theme_minimal() + 
    theme(legend.position = "none", 
          plot.title = element_text(face = "bold", size = 18, hjust = 0.5), 
          axis.title = element_text(face = "bold", size = 14), 
          axis.text.x = element_text(face = "bold", size = 11, angle = 45, hjust = 1), 
          axis.text.y = element_text(face = "bold", size = 11), 
          panel.grid.major.x = element_blank()) +
    coord_cartesian(clip = 'off') +
    labs(title = paste("CPI FOOD PRICE CHANGES (", start_yr, "TO", end_yr, ")"), 
         y = "PERCENT CHANGE (%)", x = "CATEGORY")
}

plot_cpi_change(df, 2015, 2025)
plot_cpi_change(df, 2024, 2025)
plot_cpi_change(df, 2021, 2022)

End