The CDC publishes firearm mortality for each State per 100,000 persons
https://www.cdc.gov/nchs/state-stats/deaths/firearms.html?CDC_AAref_Val=https://www.cdc.gov/nchs/pressroom/sosmap/firearm_mortality/firearm.htm. Each State’ firearm control laws can be categorized as very strict to very lax. The purpose of this Story is to answer the question, ” Do stricter firearm control laws help reduce firearm mortality?”
For this assignment you will need to:
Access the firearm mortality data from the CDC using an available API (https://open.cdc.gov/apis.html)
Create a 5 point Likert scale categorizing gun control laws from most lax to strictest and assign each state to the most appropriate Likert bin.
Determine wether stricter gun control laws result in reduced gun violence deaths
Present your story using heat maps
# load libraries
library(tidyr)
library(janitor)
Attaching package: 'janitor'The following objects are masked from 'package:stats':
chisq.test, fisher.testlibrary(ggplot2)
library(robotstxt)
library(rvest)
library(dplyr)
Attaching package: 'dplyr'The following objects are masked from 'package:stats':
filter, lagThe following objects are masked from 'package:base':
intersect, setdiff, setequal, unionlibrary(stringr)
library(usmap)
library(patchwork)#load data
mortality_data <- read.csv("mortality_data-table.csv") %>%
clean_names()
# Script to Scrape Historical Giffords Gun Law Scorecard Data (2020-2024)
# This will extract state grades and rankings from multiple years
# Define the URLs for each year
urls <- list(
"2020" = "https://giffords.org/lawcenter/resources/scorecard2020/",
"2021" = "https://giffords.org/lawcenter/resources/scorecard2021/",
"2022" = "https://giffords.org/lawcenter/resources/scorecard2022/",
"2023" = "https://giffords.org/lawcenter/resources/scorecard2023/",
"2024" = "https://giffords.org/lawcenter/resources/scorecard/"
)
# Function to scrape data from a single year
scrape_year <- function(url, year) {
cat("Scraping", year, "data from:", url, "\n")
tryCatch({
# Read the webpage
page <- read_html(url)
# Try to find the table
tables <- page %>% html_table(fill = TRUE)
if (length(tables) > 0) {
cat(" Found", length(tables), "table(s)\n")
# Look for the main data table
for (i in seq_along(tables)) {
table <- tables[[i]]
# Check if this looks like the state data table
# Should have columns for State, Grade, etc.
if (ncol(table) >= 3 && nrow(table) >= 50) {
cat(" Using table", i, "with", nrow(table), "rows and", ncol(table), "columns\n")
# Clean up the table
# Column names vary by year, so we need to be flexible
names(table) <- make.names(names(table))
# Add year column
table$Year <- year
return(table)
}
}
}
# If table method didn't work, try alternative parsing
cat(" Table method didn't work, trying alternative parsing...\n")
# Look for state data in the HTML structure
# Find all table rows
rows <- page %>% html_nodes("tr")
if (length(rows) > 0) {
cat(" Found", length(rows), "table rows\n")
# Extract data from rows
state_data <- lapply(rows, function(row) {
cells <- row %>% html_nodes("td, th") %>% html_text(trim = TRUE)
if (length(cells) >= 3) {
return(cells)
}
return(NULL)
})
# Remove NULL entries
state_data <- Filter(Negate(is.null), state_data)
if (length(state_data) > 10) { # Should have at least 10 states
cat(" Extracted data for", length(state_data), "rows\n")
# Convert to data frame
max_cols <- max(sapply(state_data, length))
state_data <- lapply(state_data, function(x) {
length(x) <- max_cols
return(x)
})
df <- as.data.frame(do.call(rbind, state_data), stringsAsFactors = FALSE)
df$Year <- year
return(df)
}
}
cat(" WARNING: Could not extract data for", year, "\n")
return(NULL)
}, error = function(e) {
cat(" ERROR scraping", year, ":", e$message, "\n")
return(NULL)
})
}
# Scrape all years
all_data <- list()
for (year in names(urls)) {
cat("\n")
result <- scrape_year(urls[[year]], year)
if (!is.null(result)) {
all_data[[year]] <- result
}
Sys.sleep(2)
}
Scraping 2020 data from: https://giffords.org/lawcenter/resources/scorecard2020/
Found 1 table(s)
Using table 1 with 50 rows and 5 columns
Scraping 2021 data from: https://giffords.org/lawcenter/resources/scorecard2021/
Found 1 table(s)
Using table 1 with 51 rows and 5 columns
Scraping 2022 data from: https://giffords.org/lawcenter/resources/scorecard2022/
Found 1 table(s)
Using table 1 with 51 rows and 5 columns
Scraping 2023 data from: https://giffords.org/lawcenter/resources/scorecard2023/
Found 1 table(s)
Using table 1 with 51 rows and 5 columns
Scraping 2024 data from: https://giffords.org/lawcenter/resources/scorecard/
Found 1 table(s)
Using table 1 with 51 rows and 5 columns# Save each year separately
cat("=== Saving Data ===\n")=== Saving Data ===for (year in names(all_data)) {
filename <- paste0("giffords_", year, "_raw.csv")
write.csv(all_data[[year]], filename, row.names = FALSE)
cat("Saved:", filename, "\n")
}Saved: giffords_2020_raw.csv
Saved: giffords_2021_raw.csv
Saved: giffords_2022_raw.csv
Saved: giffords_2023_raw.csv
Saved: giffords_2024_raw.csv # Try to combine all years if data structure is similar
if (length(all_data) > 1) {
}NULL# Return the data for further use
all_data$`2020`
# A tibble: 50 × 6
Gun.LawStrength.Ranked. State X2020Grade Gun.DeathRate.Ranked.
<int> <chr> <chr> <int>
1 37 Alabama F 5
2 42 Alaska F 1
3 45 Arizona F 16
4 39 Arkansas F 9
5 1 California A 44
6 15 Colorado C+ 18
7 3 Connecticut A- 45
8 11 Delaware B 40
9 24 Florida C- 26
10 32 Georgia F 14
# ℹ 40 more rows
# ℹ 2 more variables: Gun.DeathRate.Per.100K. <dbl>, Year <chr>
$`2021`
# A tibble: 51 × 6
Gun.Law.Strength.............Ranked. State Grade Gun.Death.Rate......…¹
<chr> <chr> <chr> <chr>
1 31 Alabama F 5
2 41 Alaska F 6
3 42 Arizona F 20
4 50 Arkansas F 8
5 1 California A 44
6 13 Colorado B 22
7 3 Connecticut A- 45
8 12 Delaware B 25
9 24 Florida C- 29
10 28 Georgia F 15
# ℹ 41 more rows
# ℹ abbreviated name: ¹Gun.Death.Rate.............Ranked.
# ℹ 2 more variables: Gun.Death.Rate.............per.100K. <chr>, Year <chr>
$`2022`
# A tibble: 51 × 6
Gun.Law.Strength........................…¹ State Grade Gun.Death.Rate......…²
<chr> <chr> <chr> <chr>
1 38 Alab… F 4
2 41 Alas… F 6
3 42 Ariz… F 17
4 50 Arka… F 8
5 1 Cali… A 43
6 14 Colo… B 18
7 3 Conn… A- 45
8 13 Dela… B 23
9 23 Flor… C- 34
10 34 Geor… F 14
# ℹ 41 more rows
# ℹ abbreviated names:
# ¹Gun.Law.Strength..................................................Ranked.,
# ²Gun.Death.Rate..................................................Ranked.
# ℹ 2 more variables:
# Gun.Death.Rate..................................................per.100K. <chr>,
# Year <chr>
$`2023`
# A tibble: 51 × 6
Gun.Law.Strength........................…¹ State Grade Gun.Death.Rate......…²
<chr> <chr> <chr> <chr>
1 35 Alab… F 4
2 40 Alas… F 7
3 41 Ariz… F 12
4 48 Arka… F 8
5 1 Cali… A 44
6 10 Colo… A- 19
7 3 Conn… A 45
8 13 Dela… B+ 39
9 24 Flor… D+ 31
10 33 Geor… F 14
# ℹ 41 more rows
# ℹ abbreviated names:
# ¹Gun.Law.Strength..................................................Ranked.,
# ²Gun.Death.Rate..................................................Ranked.
# ℹ 2 more variables:
# Gun.Death.Rate..................................................per.100K. <chr>,
# Year <chr>
$`2024`
# A tibble: 51 × 6
Gun.Law.Strength........................…¹ State Grade Gun.Death.Rate......…²
<chr> <chr> <chr> <chr>
1 34 Alab… F 3
2 40 Alas… F 5
3 41 Ariz… F 14
4 48 Arka… F 7
5 1 Cali… A 44
6 10 Colo… A- 20
7 3 Conn… A 45
8 13 Dela… A- 39
9 24 Flor… C- 30
10 31 Geor… F 13
# ℹ 41 more rows
# ℹ abbreviated names:
# ¹Gun.Law.Strength..................................................Ranked.,
# ²Gun.Death.Rate..................................................Ranked.
# ℹ 2 more variables:
# Gun.Death.Rate..................................................per.100K. <chr>,
# Year <chr># Extract each year as its own data frame
data_2020 <- as.data.frame(all_data[["2020"]])
data_2021 <- as.data.frame(all_data[["2021"]])
data_2022 <- as.data.frame(all_data[["2022"]])
data_2023 <- as.data.frame(all_data[["2023"]])
data_2024 <- as.data.frame(all_data[["2024"]])# See what column names each year has
lapply(all_data, names)$`2020`
[1] "Gun.LawStrength.Ranked." "State"
[3] "X2020Grade" "Gun.DeathRate.Ranked."
[5] "Gun.DeathRate.Per.100K." "Year"
$`2021`
[1] "Gun.Law.Strength.............Ranked."
[2] "State"
[3] "Grade"
[4] "Gun.Death.Rate.............Ranked."
[5] "Gun.Death.Rate.............per.100K."
[6] "Year"
$`2022`
[1] "Gun.Law.Strength..................................................Ranked."
[2] "State"
[3] "Grade"
[4] "Gun.Death.Rate..................................................Ranked."
[5] "Gun.Death.Rate..................................................per.100K."
[6] "Year"
$`2023`
[1] "Gun.Law.Strength..................................................Ranked."
[2] "State"
[3] "Grade"
[4] "Gun.Death.Rate..................................................Ranked."
[5] "Gun.Death.Rate..................................................per.100K."
[6] "Year"
$`2024`
[1] "Gun.Law.Strength..................................................Ranked."
[2] "State"
[3] "Grade"
[4] "Gun.Death.Rate..................................................Ranked."
[5] "Gun.Death.Rate..................................................per.100K."
[6] "Year" # Function to standardize column names AND data types
standardize_columns <- function(df, year) {
# Get current column names
cols <- names(df)
# Create new standardized names
new_names <- cols
# Standardize each column
for (i in seq_along(cols)) {
col <- cols[i]
# Gun Law Rank
if (grepl("Gun.*Law.*Strength.*Ranked", col, ignore.case = TRUE)) {
new_names[i] <- "Gun_Law_Rank"
}
# State (keep as is)
else if (grepl("^State$", col)) {
new_names[i] <- "State"
}
# Grade
else if (grepl("Grade", col, ignore.case = TRUE)) {
new_names[i] <- "Grade"
}
# Gun Death Rank
else if (grepl("Gun.*Death.*Rate.*Ranked", col, ignore.case = TRUE)) {
new_names[i] <- "Gun_Death_Rank"
}
# Gun Death Rate (per 100K)
else if (grepl("Gun.*Death.*Rate.*100K", col, ignore.case = TRUE)) {
new_names[i] <- "Gun_Death_Rate"
}
# Year
else if (grepl("^Year$", col)) {
new_names[i] <- "Year"
}
}
# Apply new names
names(df) <- new_names
# Select columns and FIX DATA TYPES
df <- df %>%
select(State, Grade, Gun_Law_Rank, Gun_Death_Rank, Gun_Death_Rate, Year) %>%
mutate(
State = as.character(State),
Grade = as.character(Grade),
Gun_Law_Rank = as.numeric(as.character(Gun_Law_Rank)), # Convert to numeric
Gun_Death_Rank = as.numeric(as.character(Gun_Death_Rank)), # Convert to numeric
Gun_Death_Rate = as.numeric(as.character(Gun_Death_Rate)), # Convert to numeric
Year = as.numeric(as.character(Year)) # Convert to numeric
) %>%
# Remove whitespace
mutate(across(where(is.character), trimws)) %>%
# Remove any header rows that snuck in
filter(!is.na(Gun_Law_Rank) & State != "State" & State != "")
return(df)
}
# Apply to all years
cleaned_data <- list()
for (year in names(all_data)) {
cat("Cleaning", year, "data...\n")
cleaned_data[[year]] <- standardize_columns(all_data[[year]], year)
}Cleaning 2020 data...
Cleaning 2021 data...Warning: There were 3 warnings in `mutate()`.
The first warning was:
ℹ In argument: `Gun_Law_Rank = as.numeric(as.character(Gun_Law_Rank))`.
Caused by warning:
! NAs introduced by coercion
ℹ Run `dplyr::last_dplyr_warnings()` to see the 2 remaining warnings.Cleaning 2022 data...Warning: There were 3 warnings in `mutate()`.
The first warning was:
ℹ In argument: `Gun_Law_Rank = as.numeric(as.character(Gun_Law_Rank))`.
Caused by warning:
! NAs introduced by coercion
ℹ Run `dplyr::last_dplyr_warnings()` to see the 2 remaining warnings.Cleaning 2023 data...Warning: There were 3 warnings in `mutate()`.
The first warning was:
ℹ In argument: `Gun_Law_Rank = as.numeric(as.character(Gun_Law_Rank))`.
Caused by warning:
! NAs introduced by coercion
ℹ Run `dplyr::last_dplyr_warnings()` to see the 2 remaining warnings.Cleaning 2024 data...Warning: There were 3 warnings in `mutate()`.
The first warning was:
ℹ In argument: `Gun_Law_Rank = as.numeric(as.character(Gun_Law_Rank))`.
Caused by warning:
! NAs introduced by coercion
ℹ Run `dplyr::last_dplyr_warnings()` to see the 2 remaining warnings.# Now combine all years
combined_data <- bind_rows(cleaned_data)
combined_data <- combined_data %>%
clean_names()# Add Likert scale (1-5) based on rankings
combined_data <- combined_data %>%
mutate(
likert_scale = cut(gun_law_rank,
breaks = c(0, 10, 20, 30, 40, 50),
labels = c(5, 4, 3, 2, 1),
include.lowest = TRUE),
likert_scale = as.numeric(as.character(likert_scale)),
category = cut(gun_law_rank,
breaks = c(0, 10, 20, 30, 40, 50),
labels = c("Very Strict", "Strict", "Moderate", "Lax", "Very Lax"),
include.lowest = TRUE)
)
# Create a lookup table
state_lookup <- data.frame(
state = c("Alabama", "Alaska", "Arizona", "Arkansas", "California",
"Colorado", "Connecticut", "Delaware", "Florida", "Georgia",
"Hawaii", "Idaho", "Illinois", "Indiana", "Iowa",
"Kansas", "Kentucky", "Louisiana", "Maine", "Maryland",
"Massachusetts", "Michigan", "Minnesota", "Mississippi", "Missouri",
"Montana", "Nebraska", "Nevada", "New Hampshire", "New Jersey",
"New Mexico", "New York", "North Carolina", "North Dakota", "Ohio",
"Oklahoma", "Oregon", "Pennsylvania", "Rhode Island", "South Carolina",
"South Dakota", "Tennessee", "Texas", "Utah", "Vermont",
"Virginia", "Washington", "West Virginia", "Wisconsin", "Wyoming"),
state_code = c("AL", "AK", "AZ", "AR", "CA", "CO", "CT", "DE", "FL", "GA",
"HI", "ID", "IL", "IN", "IA", "KS", "KY", "LA", "ME", "MD",
"MA", "MI", "MN", "MS", "MO", "MT", "NE", "NV", "NH", "NJ",
"NM", "NY", "NC", "ND", "OH", "OK", "OR", "PA", "RI", "SC",
"SD", "TN", "TX", "UT", "VT", "VA", "WA", "WV", "WI", "WY")
)
# Add state codes to combined_data
combined_data <- combined_data %>%
left_join(state_lookup, by = "state")#joining datasets for the
complete_df <- combined_data %>%
inner_join(mortality_data, by = c("state_code" = "state", "year"))# Calculate trend data
trend_plot <- complete_df %>%
filter(year %in% 2020:2023) %>%
group_by(year, category) %>%
summarize(avg_rate = mean(rate, na.rm = TRUE), .groups = "drop") %>%
mutate(category = factor(category,
levels = c("Very Strict", "Strict", "Moderate", "Lax", "Very Lax")))
# Create improved plot
ggplot(trend_plot, aes(x = year, y = avg_rate, color = category, group = category)) +
geom_line(size = 1.2) +
geom_point(size = 3) +
geom_text(
data = trend_plot %>% filter(year == 2023),
aes(label = round(avg_rate, 1)),
hjust = -0.3,
size = 3.5,
fontface = "bold",
show.legend = FALSE
) +
scale_color_manual(
name = "Gun Law Category",
values = c(
"Very Strict" = "#08519c",
"Strict" = "#3182bd",
"Moderate" = "#9ecae1",
"Lax" = "#fc9272",
"Very Lax" = "#a50f15"
),
breaks = c("Very Strict", "Strict", "Moderate", "Lax", "Very Lax")
) +
scale_x_continuous(breaks = 2020:2023, limits = c(2020, 2023.5)) +
scale_y_continuous(limits = c(0, 25)) +
labs(
title = "Firearm Mortality Trends by Gun Law Strictness",
subtitle = "States with stricter gun laws consistently show lower and more stable mortality rates (2020-2023)",
x = "Year",
y = "Average Mortality Rate (per 100,000)"
) +
theme_minimal() +
theme(
plot.title = element_text(size = 14, face = "bold", hjust = 0),
plot.subtitle = element_text(size = 11, color = "gray40", hjust = 0),
legend.position = "right",
legend.title = element_text(face = "bold"),
panel.grid.minor = element_blank(),
axis.title = element_text(face = "bold")
)Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
The mortality rate gap between states with very strict laws (~8-9 per 100k) and very lax laws (~20-22 per 100k) remains remarkably consistent across all four years.
# Filter for 2023
plot_data <- complete_df %>% filter(year == 2023)
# Heat map: Gun Law Strictness (1 = Very Lax → 5 = Very Strict)
heat_gun_laws <- plot_usmap(data = plot_data, values = "likert_scale") +
scale_fill_gradientn(
colors = c("#08519c", "#3182bd", "#6baed6", "#bdd7e7", "#eff3ff"),
name = "Gun Law Strictness",
breaks = 1:5,
labels = c("1\n(Least strict)", "2", "3", "4", "5\n(Most strict)")
) +
theme_void() +
theme(
legend.position = "bottom",
legend.title = element_text(size = 9, face = "bold"),
legend.text = element_text(size = 5)
)
# Heat map: Firearm Mortality Rate
heat_mortality <- plot_usmap(data = plot_data, values = "rate") +
scale_fill_gradientn(
colors = c("#fff5f0", "#fcbba1", "#fc9272", "#fb6a4a", "#de2d26", "#a50f15"),
name = "Mortality Rate (per 100k)",
n.breaks = 6
) +
theme_void() +
theme(
legend.position = "bottom",
legend.title = element_text(size = 9, face = "bold"),
legend.text = element_text(size = 5)
)
# Combine side by side with improved title
combined_heatmaps <- heat_gun_laws + heat_mortality +
plot_layout(ncol = 2) +
plot_annotation(
title = "Stricter Gun Laws, Lower Mortality Rates",
subtitle = "2023 Data",
theme = theme(
plot.title = element_text(size = 15, hjust = 0.5, face = "bold"),
plot.subtitle = element_text(size = 11, hjust = 0.5, color = "gray40"),
plot.caption = element_text(size = 9, hjust = 1, color = "gray50")
)
)
# Display
combined_heatmapsWarning: annotation$theme is not a valid theme.
Please use `theme()` to construct themes.
# Top and bottom states
extreme_states <- plot_data %>%
arrange(desc(rate)) %>%
mutate(rank = row_number()) %>%
filter(rank <= 10 | rank > n() - 10) %>%
mutate(category = ifelse(rank <= 10, "Highest Mortality", "Lowest Mortality"))
ggplot(extreme_states, aes(x = reorder(state, rate), y = rate, fill = likert_scale)) +
geom_col() +
geom_text(aes(label = round(rate, 1)), hjust = -0.2, size = 3) +
coord_flip() +
scale_fill_gradientn(
colors = c("#08519c", "#3182bd", "#6baed6", "#bdd7e7", "#eff3ff"),
name = "Gun Law\nStrictness",
breaks = 1:5,
labels = c("1 (Least)", "2", "3", "4", "5 (Most)"),
limits = c(1, 5)
) +
labs(
title = "States with Highest and Lowest Firearm Mortality Rates",
x = NULL,
y = "Mortality Rate (per 100k)"
) +
theme_minimal() +
theme(
plot.title = element_text(size = 14, face = "bold"),
legend.position = "right",
legend.key.height = unit(1, "cm"),
legend.key.width = unit(0.5, "cm"),
legend.title = element_text(size = 10, face = "bold"),
legend.text = element_text(size = )
)
library(ggrepel)
# Define the color palette once
my_colors <- c("#08519c", "#3182bd", "#6baed6", "#bdd7e7", "#eff3ff")
outliers <- plot_data %>%
group_by(likert_scale) %>%
mutate(
Q1 = quantile(rate, 0.25),
Q3 = quantile(rate, 0.75),
IQR = Q3 - Q1,
is_outlier = rate < (Q1 - 1.5 * IQR) | rate > (Q3 + 1.5 * IQR)
) %>%
filter(is_outlier)
ggplot(plot_data, aes(x = factor(likert_scale), y = rate)) +
geom_boxplot(aes(fill = factor(likert_scale)), alpha = 0.7) +
geom_jitter(width = 0.2, alpha = 0.5, size = 2, color = "gray40") +
geom_text_repel(
data = outliers,
aes(label = state),
size = 3,
max.overlaps = 20
) +
scale_fill_manual(values = my_colors) + # CHANGED: Use same colors
labs(
title = "Stricter Gun Laws Associated with Lower Firearm Mortality Rates",
x = "Gun Law Strictness (1 = Least Strict, 5 = Most Strict)",
y = "Mortality Rate (per 100k)"
) +
theme_minimal() +
theme(legend.position = "none")
This analysis examined the relationship between gun law strictness and firearm mortality rates across all 50 US states using data from the CDC and the Giffords Law Center spanning 2020-2023.
The visual evidence from the heat maps reveals a clear geographic pattern: states with the strictest gun laws (rated 4-5 on Likert scale) consistently show lower firearm mortality rates, often below 10 deaths per 100,000 people. In stark contrast, states with the least restrictive gun laws (rated 1-2), concentrated in the South and Mountain West, experience significantly higher mortality rates, frequently exceeding 20 deaths per 100,000 people.
The bar chart analysis reinforces this pattern. Among the 10 states with the highest mortality rates, 9 have gun law strictness ratings of 1 or 2 (Very Lax or Lax). Conversely, among the 10 states with the lowest mortality rates, the majority have strictness ratings of 4 or 5 (Strict or Very Strict). The boxplot distribution further illustrates this relationship, showing a consistent downward trend in median firearm mortality rates as gun law strictness increases from 1 to 5.
The trend analysis from 2020-2023 demonstrates that this relationship is not only strong but also persistent over time. The mortality rate gap between states with very strict laws (averaging 7-9 per 100,000) and very lax laws (averaging 19-22 per 100,000) remains remarkably consistent across all four years. Notably, states with the strictest gun laws also exhibit the least volatility in mortality rates, suggesting these policies may provide more stable public health outcomes. While all categories experienced an increase in 2021, likely reflecting broader societal impacts of the COVID-19 pandemic, the relative ordering remained unchanged.
Do stricter gun laws reduce firearm deaths? Based on this state-level analysis, the evidence strongly suggests yes. States with stricter gun control legislation demonstrate substantially lower firearm mortality rates compared to states with permissive gun laws, and this pattern holds consistently across multiple years.