Mini-Assignment 5
Submitted by Uday Bhanu Singh
Setup
#Loading CSV data
feData <- read_csv("fatalEncounters.csv")
names(feData) <- c(
"id", "name", "age", "gender", "race",
"race_imputed", "imputation_prob", "image_url",
"date", "injury_address", "death_city", "state",
"death_zip", "death_county", "address",
"latitude", "longitude", "agencies",
"force", "uid_temp", "name_temp",
"armed", "weapon", "aggressive",
"fleeing", "description_temp", "url_temp",
"description", "disposition", "intended_force",
"document_url", "empty32", "empty33",
"id_formula", "id_redundant"
)
head(feData)## # A tibble: 6 × 35
## id name age gender race race_imputed imputation_prob image_url date
## <dbl> <chr> <dbl> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 31495 Ashley … 28 Female Afri… African-Ame… Not imputed https://… 12/3…
## 2 31496 Name wi… NA Female Race… <NA> <NA> <NA> 12/3…
## 3 31497 Name wi… NA Male Race… <NA> <NA> <NA> 12/3…
## 4 31491 Johnny … 36 Male Race… <NA> <NA> <NA> 12/3…
## 5 31492 Dennis … 44 Male Euro… <NA> <NA> <NA> 12/3…
## 6 31493 Ny'Dari… 21 Male Race… <NA> <NA> <NA> 12/3…
## # ℹ 26 more variables: injury_address <chr>, death_city <chr>, state <chr>,
## # death_zip <chr>, death_county <chr>, address <chr>, latitude <chr>,
## # longitude <dbl>, agencies <chr>, force <chr>, uid_temp <dbl>,
## # name_temp <chr>, armed <chr>, weapon <chr>, aggressive <chr>,
## # fleeing <chr>, description_temp <chr>, url_temp <chr>, description <chr>,
## # disposition <chr>, intended_force <chr>, document_url <chr>, empty32 <lgl>,
## # empty33 <lgl>, id_formula <dbl>, id_redundant <dbl>ggplot(feData, aes(x=gender)) +
geom_bar() + geom_text(
aes(label = after_stat(count)),
stat = "count",
vjust = -0.3,
)
Making a basic chart to test the data. I first examined the gender distribution of individuals recorded in the Fatal Encounters database. Men account for approximately 90% of recorded fatalities, indicating a significant gender divide in the statistic.
Transforming latitude and longitude columns into spatial date.
fe_geom <- feData |>
mutate(
latitude = as.numeric(latitude),
longitude = as.numeric(longitude)
) |>
filter(
!is.na(latitude),
!is.na(longitude),
between(latitude, -90, 90),
between(longitude, -180, 180)) |>
st_as_sf(
coords = c("longitude", "latitude"),
crs = 4326,
remove = FALSE
)Importing ACS data to examine demographic trends in the data at the county level.
# County-level ACS demographic data
acs_vars <- c(
pop = "B01003_001",
poverty_total = "B17001_001",
poverty_count = "B17001_002",
race_total = "B03002_001",
white_nh = "B03002_003",
black_nh = "B03002_004",
asian_nh = "B03002_006",
hispanic = "B03002_012"
)
countyData <- get_acs(
geography = "county",
variables = acs_vars,
year = 2024,
survey = "acs5",
state = NULL,
output = "wide",
geometry = TRUE
) |>
mutate(
poverty_pct = 100 * poverty_countE / poverty_totalE,
white_pct = 100 * white_nhE / race_totalE,
black_pct = 100 * black_nhE / race_totalE,
asian_pct = 100 * asian_nhE / race_totalE,
hispanic_pct = 100 * hispanicE / race_totalE
) |>
select(
GEOID, NAME, geometry,
pop = popE,
poverty_pct, white_pct, black_pct,
asian_pct, hispanic_pct
) |>
st_transform(4326)Also getting the same data at state level:
stateData <- get_acs(
geography = "state",
variables = acs_vars,
year = 2024,
survey = "acs5",
output = "wide",
geometry = TRUE
) |>
mutate(
poverty_pct = 100 * poverty_countE / poverty_totalE,
white_pct = 100 * white_nhE / race_totalE,
black_pct = 100 * black_nhE / race_totalE,
asian_pct = 100 * asian_nhE / race_totalE,
hispanic_pct = 100 * hispanicE / race_totalE
) |>
select(
GEOID, NAME, geometry,
pop = popE,
poverty_pct, white_pct, black_pct,
asian_pct, hispanic_pct
) |>
st_transform(4326)Performing spatial joins to associate each row in the Fatal Encounters data to a county
# Assign each fatal encounter to a county
fe_county <- st_join(
fe_geom,
countyData |> select(GEOID, NAME),
join = st_intersects,
left = FALSE
)
# Count fatal encounters by county
fe_counts <- fe_county |>
st_drop_geometry() |>
count(GEOID, name = "fe_count")
# Join counts to county demographics
countyFE <- countyData |>
left_join(fe_counts, by = "GEOID") |>
mutate(
fe_count = replace_na(fe_count, 0L),
fe_per_10k = if_else(pop > 0, 10000 * fe_count / pop, NA_real_)
)I am using fatal encounter rate per 10K of the poplation to account for population variation.
Creating a seperate data frame for state level aggregation:
stateFE <- countyFE |>
mutate(
STATEFP = substr(GEOID, 1, 2)
) |>
st_drop_geometry() |>
group_by(STATEFP) |>
summarise(
fe_count = sum(fe_count, na.rm = TRUE),
.groups = "drop"
) |>
right_join(stateData, by = c("STATEFP" = "GEOID")) |>
mutate(
STATE = NAME,
fe_count = replace_na(fe_count, 0L),
fe_per_10k = if_else(pop > 0, 10000 * fe_count / pop, NA_real_)
) |>
st_as_sf()Racial Demographics Analysis
countyFE |>
st_drop_geometry() |>
filter(pop >= 20000) |> #Creating a population filter for meaningful analysis
pivot_longer(
cols = c(white_pct, black_pct, asian_pct, hispanic_pct),
names_to = "race",
values_to = "percentage"
) |>
ggplot(aes(x = percentage, y = fe_per_10k)) +
geom_point(alpha = 0.2, size = 0.8) +
geom_smooth(method = "lm", se = FALSE) +
facet_wrap(~race, scales = "free_x") +
labs(
x = "Population Share (%)",
y = "Fatal Encounters per 10,000",
title = "Racial Composition and Fatal Encounter Rates"
) +
theme_minimal()+
stat_cor()
Counties with larger Black population shares exhibit a moderate positive correlation with fatal encounter rates (R = 0.30), while White (R = −0.12) and Hispanic (R = −0.11) population shares show weak negative correlations. The relationship with Asian population share is negligible and not statistically significant (R = −0.043, p = 0.061). Overall, racial composition is associated with geographic variation in fatal encounters, although these county-level relationships do not establish individual racial disparities or causation.
Time-Series Analysis
feTime <- feData |>
mutate(
date = mdy(date),
year = year(date),
month_year = floor_date(date, unit = "month")
) |>
filter(!is.na(year) & (year != 1999)) |> #Filtering out 1999 as it is only partially present in the data
count(year, month_year, name = "fe_count")
ggplot(feTime, aes(x=month_year, fe_count)) + geom_point() +
geom_line(color = "steelblue", alpha = 0.4, linetype = "dashed") +
labs(
x = "Month",
y = "Fatal Encounters",
title = "Monthly Fatal Encounters in the United States"
) +
geom_smooth(method = "lm", se = TRUE, color = "red")+
stat_cor()+
theme_minimal()
feYearly <- feTime |>
group_by(year) |>
summarise(fe_count = sum(fe_count), .groups = "drop")
ggplot(feYearly, aes(x= year, fe_count)) + geom_point() +
geom_line(color = "steelblue", alpha = 0.4, linetype = "dashed") +
scale_x_continuous(breaks = seq(2000, 2022, by = 2)) +
labs(
x = "Year",
y = "Fatal Encounters",
title = "Yearly Fatal Encounters in the United States"
) +
geom_smooth(method = "lm", se = TRUE, color = "red")+
stat_cor()+
theme_minimal()
There is a definite increase in fatal encounter reports in the past two decades as per the data. The slope shows more than a two-fold increase, as well as statistical significance.
However, I am unsure about if this is an artifact of the dataset getting more reports, or a genuine increase.
To address this concern, I decided to cross check the data with another dataset. Mapping Police Violence is the closest alternative I could find online:
Mapping Police Violence Comparison
#Cleaning MPC data
mpvData <- read_csv("Mapping Police Violence.csv") |>
mutate(
latitude = as.numeric(latitude),
longitude = as.numeric(longitude)
) |>
filter(
!is.na(latitude),
!is.na(longitude),
between(latitude, -90, 90),
between(longitude, -180, 180)) |>
st_as_sf(
coords = c("longitude", "latitude"),
crs = 4326,
remove = FALSE
)
feTimeMPV <- mpvData |>
mutate(
date = mdy(date),
year = year(date),
month_year = floor_date(date, unit = "month")
) |>
filter(!is.na(year) & (year != 2026)) |>
count(year, month_year, name = "fe_count")
ggplot(feTimeMPV, aes(x=month_year, fe_count)) + geom_point() +
geom_line(color = "steelblue", alpha = 0.4, linetype = "dashed") +
labs(
x = "Month",
y = "Fatal Encounters",
title = "Monthly Fatal Encounters in the United States"
) +
geom_smooth(method = "lm", se = TRUE, color = "red")+
stat_cor()+
theme_minimal()
feYearlyMPV <- feTimeMPV |>
group_by(year) |>
summarise(fe_count = sum(fe_count), .groups = "drop")
ggplot(feYearlyMPV, aes(x= year, fe_count)) + geom_point() +
geom_line(color = "steelblue", alpha = 0.4, linetype = "dashed") +
scale_x_continuous(breaks = seq(2000, 2025, by = 2)) +
labs(
x = "Year",
y = "Fatal Encounters",
title = "Yearly Fatal Encounters in the United States"
) +
geom_smooth(method = "lm", se = TRUE, color = "red")+
stat_cor()+
theme_minimal()
Mapping Police Violence also suggests a statistically significant
increase in police violence, although lower in magnitude.
# Restricting both datasets to 2013–2021
feMPV_filtered <- feYearlyMPV |>
filter(between(year, 2013, 2021))
fe_filtered <- feYearly |>
filter(between(year, 2013, 2021))
# Calculating regression slopes
slope_mpv <- coef(lm(fe_count ~ year, data = feMPV_filtered))[2]
slope_fe <- coef(lm(fe_count ~ year, data = fe_filtered))[2]
ggplot() +
geom_point(data = feMPV_filtered,
aes(x = year, y = fe_count), color = "steelblue") +
geom_line(data = feMPV_filtered,
aes(x = year, y = fe_count),
color = "steelblue", alpha = 0.4, linetype = "dashed") +
geom_smooth(data = feMPV_filtered,
aes(x = year, y = fe_count),
method = "lm", se = TRUE, color = "steelblue") +
ggpubr::stat_cor(data = feMPV_filtered,
aes(x = year, y = fe_count),
color = "steelblue",
label.x.npc = 0.05, label.y.npc = 0.95) +
geom_point(data = fe_filtered,
aes(x = year, y = fe_count), color = "firebrick") +
geom_line(data = fe_filtered,
aes(x = year, y = fe_count),
color = "firebrick", alpha = 0.4, linetype = "dashed") +
geom_smooth(data = fe_filtered,
aes(x = year, y = fe_count),
method = "lm", se = TRUE, color = "firebrick") +
ggpubr::stat_cor(data = fe_filtered,
aes(x = year, y = fe_count),
color = "firebrick",
label.x.npc = 0.05, label.y.npc = 0.85) +
scale_x_continuous(breaks = 2013:2021) +
labs(
x = "Year",
y = "Annual Deaths",
title = "Police-Related Deaths in the United States (2013–2021)",
subtitle = "Blue: Mapping Police Violence | Red: Fatal Encounters"
) +
annotate(
"label",
x = 2013.5, y = 1300,
label = paste0("MPV slope: ", round(slope_mpv, 2), " deaths/year"),
color = "steelblue",
hjust = 0, vjust = 1.5
) +
annotate(
"label",
x = 2013.5, y = 1700,
label = paste0("FE slope: ", round(slope_fe, 2), " deaths/year"),
color = "firebrick",
hjust = 0, vjust = 3.5
) +
theme_minimal()
# Calculating regression slopes
slope_mpv <- coef(lm(fe_count ~ year, data = feYearlyMPV))[2]
slope_fe <- coef(lm(fe_count ~ year, data = feYearly))[2]
ggplot() +
geom_point(data = feYearlyMPV,
aes(x = year, y = fe_count), color = "steelblue") +
geom_line(data = feYearlyMPV,
aes(x = year, y = fe_count),
color = "steelblue", alpha = 0.4, linetype = "dashed") +
geom_smooth(data = feYearlyMPV,
aes(x = year, y = fe_count),
method = "lm", se = TRUE, color = "steelblue") +
ggpubr::stat_cor(data = feYearlyMPV,
aes(x = year, y = fe_count),
color = "steelblue",
label.x.npc = 0.05, label.y.npc = 0.95) +
geom_point(data = feYearly,
aes(x = year, y = fe_count), color = "firebrick") +
geom_line(data = feYearly,
aes(x = year, y = fe_count),
color = "firebrick", alpha = 0.4, linetype = "dashed") +
geom_smooth(data = feYearly,
aes(x = year, y = fe_count),
method = "lm", se = TRUE, color = "firebrick") +
ggpubr::stat_cor(data = feYearly,
aes(x = year, y = fe_count),
color = "firebrick",
label.x.npc = 0.05, label.y.npc = 0.85) +
scale_x_continuous(breaks = seq(2000, 2025, by = 2)) +
labs(
x = "Year",
y = "Annual Deaths",
title = "Police-Related Deaths in the United States",
subtitle = "Blue: Mapping Police Violence | Red: Fatal Encounters"
) +
annotate(
"label",
x = 2005, y = 950,
label = paste0("MPV slope: ", round(slope_mpv, 2), " deaths/year"),
color = "steelblue",
hjust = 0, vjust = 1.5
) +
annotate(
"label",
x = 2000, y = 1800,
label = paste0("FE slope: ", round(slope_fe, 2), " deaths/year"),
color = "firebrick",
hjust = 0, vjust = 3.5
) +
theme_minimal()
Fatal Encounters and Mapping Police Violence both document substantial increases in recorded police-related deaths in the United States. The consistency in trend direction suggests that the increase is not solely an artifact of Fatal Encounters’ reporting practices. However, because the datasets share historical sources and use different inclusion criteria, improvements in data coverage cannot be ruled out as a contributing factor.
The rest of the analysis will only use Fatal Encounters data.
# Count fatal encounters by state and year
fe_state_counts <- fe_county |>
st_drop_geometry() |>
mutate(
STATEFP = substr(GEOID, 1, 2),
year = year(mdy(date))
) |>
filter(between(year, 2000, 2021)) |>
count(STATEFP, year) |>
pivot_wider(
names_from = year,
values_from = n,
names_prefix = "fe_count_",
values_fill = 0
)
# Join counts to state geometry and calculate percentage change
stateFE_change <- stateFE |>
left_join(fe_state_counts, by = "STATEFP") |>
mutate(
across(starts_with("fe_count_"), ~ replace_na(.x, 0L))
) |>
mutate(
fe_count_change = if_else(
fe_count_2000 > 0,
round(
100 * (fe_count_2021 - fe_count_2000) / fe_count_2000,
2
),
NA_real_
)
)
#Pivoting to plot data
stateFE_change_long <- stateFE_change |>
st_drop_geometry() |>
select(-fe_count) |>
pivot_longer(
cols = all_of(paste0("fe_count_", 2000:2021)),
names_to = "year",
names_prefix = "fe_count_",
values_to = "fe_count"
) |>
mutate(
year = as.integer(year),
fe_per_10k = if_else(
pop > 0,
10000 * fe_count / pop,
NA_real_
)
)State-wise Heatmap
ggplot(
stateFE_change_long,
aes(
x = year,
y = reorder(STATE, fe_per_10k, FUN = mean, na.rm = TRUE),
fill = fe_per_10k
)
) +
geom_tile(color = "white", linewidth = 0.1) +
scale_fill_viridis_c(
option = "magma",
direction = -1,
trans = "sqrt",
breaks = c(0, 0.02, 0.05, 0.10, 0.15, 0.20),
labels = scales::label_number(accuracy = 0.01),
name = "Deaths per 10k"
) +
scale_x_continuous(
breaks = seq(2000, 2021, 2),
expand = c(0, 0)
) +
scale_y_discrete(expand = c(0, 0)) +
labs(
title = "Annual recorded deaths per 10,000 residents, 2000–2021",
x = "Year",
y = NULL
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(face = "bold", size = 16),
plot.subtitle = element_text(size = 11, color = "gray40"),
axis.text.y = element_text(size = 8),
axis.text.x = element_text(size = 9),
panel.grid = element_blank(),
legend.position = "bottom",
legend.key.width = unit(2, "cm"),
plot.margin = margin(10, 15, 10, 10)
) +
guides(
fill = guide_colorbar(
title.position = "top",
barwidth = unit(10, "cm"),
barheight = unit(0.4, "cm")
)
)
I examined annual fatal encounter rates across states from 2000 to 2021 to identify geographic and temporal variation. The heatmap reveals substantial differences between states, with some exhibiting consistently elevated rates and others showing considerable year-to-year fluctuations. To narrow the analysis, I selected states with both above-average fatal encounter rates and above-average rates of increase for further county-level investigation.
#Creating filter to select states for analysis
stateTrends <- stateFE_change_long |>
filter(pop > 1500000) |> #Filtering out low-population states
group_by(STATE) |>
summarise(
slope = coef(lm(fe_per_10k ~ year))[2],
mean_rate = mean(fe_per_10k, na.rm = TRUE),
.groups = "drop"
) |>
mutate(
above_avg_rate = mean_rate > mean(mean_rate, na.rm = TRUE), #Filter for above average fatal encounter rate
above_avg_slope = slope > mean(slope, na.rm = TRUE), #Filter for above average increase
selected = above_avg_rate & above_avg_slope
) |>
arrange(desc(mean_rate))
selectedStates <- stateTrends |>
filter(selected)
ggplot(stateTrends, aes(x = mean_rate, y = slope)) +
geom_point(aes(colour = selected), size = 2.5, alpha = 0.8) +
geom_text_repel(
data = filter(stateTrends, selected),
aes(label = STATE),
size = 3,
max.overlaps = Inf
) +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_vline(xintercept = mean(stateTrends$mean_rate),
linetype = "dashed", alpha = 0.1) +
geom_hline(yintercept = mean(stateTrends$slope),
linetype = "dashed", alpha = 0.1) +
labs(
x = "Mean Annual Fatal Encounters per 10,000",
y = "Annual Change in Rate (Regression Slope)",
title = "Identifying States for Further Analysis"
) +
scale_color_manual(
values = c(`TRUE` = "firebrick", `FALSE` = "steelblue"),
labels = c(`TRUE` = "Selected States", `FALSE` = "Other States"),
name = NULL
) +
theme_minimal()
To identify states for further analysis, I calculated each state’s mean annual fatal encounter rate and estimated its rate of change using linear regression from 2000 to 2021. I restricted the comparison to states with populations exceeding 1.5 million and selected those exceeding both the average rate and average regression slope. New Mexico emerged as a particularly notable outlier, exhibiting the highest mean rate and steepest increase among the states considered. The 12 selected states were subsequently examined at the county level to investigate geographic and socioeconomic variation.
Summary of Selected States
selectedStateList <- selectedStates$STATE
selectedStatesFE <- stateFE |>
filter((STATE %in% selectedStateList)) |>
inner_join((selectedStates |>
select(STATE, slope, mean_rate)), by = "STATE")
ggplot(selectedStatesFE, aes(x=pop, y= reorder(STATE, pop))) + geom_col() +
theme_minimal()+
labs(
title = "Population by State",
x = "Population",
y = "State"
)
Demographic Comparison of Selected States
# Calculating population density for all counties
countyFE <- countyFE |>
mutate(
area_km2 = as.numeric(st_area(geometry)) / 1e6,
pop_density = pop / area_km2
)
selectedCountiesFE <- countyFE |>
mutate(
STATE = str_trim(sub(".*,", "", NAME))
) |>
filter(STATE %in% selectedStateList)
# Demographic composition
selectedStatesFE |>
st_drop_geometry() |>
pivot_longer(
cols = c(white_pct, black_pct, asian_pct, hispanic_pct, poverty_pct),
names_to = "demcol",
values_to = "percentage"
) |>
ggplot(aes(
x = percentage,
y = reorder(STATE, pop),
color = demcol
)) +
geom_point(position = position_dodge(width = 0.5), size = 2) +
theme_minimal() +
labs(
title = "Demographic Composition by State",
x = "Population (%)",
y = "State",
color = "Demographic Measure"
)
# Fatal encounters per 10,000 residents
selectedStatesFE |>
st_drop_geometry() |>
ggplot(aes(
x = fe_per_10k,
y = reorder(STATE, pop)
)) +
geom_col(fill = "steelblue", width = 0.7) +
theme_minimal() +
labs(
title = "Fatal Encounter Rates by State (Ordered by Population)",
x = "Fatal Encounters per 10,000 Residents",
y = "State"
)
New Mexico exhibits the highest rate (approximately 2.1 per 10,000 residents), followed by Oklahoma and Mississippi. Despite these differences, the states vary considerably in racial/ethnic composition and poverty levels, with no immediately apparent demographic pattern explaining the variation in fatal encounter rates.
Identifying Counties with Highest Fatality Rates
selectedCountiesFE |>
st_drop_geometry() |>
filter(pop >= 20000) |>
slice_max(fe_per_10k, n = 15) |>
ggplot(aes(
x = fe_per_10k,
y = reorder(NAME, fe_per_10k)
)) +
geom_col(fill = "steelblue") +
labs(
title = "Counties with the Highest Fatal Encounter Rates",
subtitle = "Selected states, counties with 20,000+ residents",
x = "Fatal Encounters per 10,000 Residents",
y = NULL
) +
theme_minimal()
Southern states’ counties constitute most number of fatalities.
Map of Selected States’ Counties
tmap_mode("view")
tm_shape(selectedCountiesFE) +
tm_polygons(
fill = "fe_per_10k",
fill.scale = tm_scale_continuous(
values = "brewer.oranges"
),
fill.legend = tm_legend("Fatal Encounters per 10,000"),
col = "white",
lwd = 0.3
) +
tm_layout(
title = "Fatal Encounter Rates in Selected States",
legend.outside = TRUE
)The map reveals considerable differences between counties, with several localized areas exhibiting substantially higher rates than surrounding counties. This suggests that statewide averages may obscure important geographic disparities, motivating further investigation of county-level socioeconomic and demographic characteristics.
County Level White vs Non-white Comparison
# Classifying victim race and calculate county-level proportions
fe_race_counts <- fe_county |>
st_drop_geometry() |>
mutate(
race_clean = str_to_lower(str_trim(race)),
race_group = case_when(
is.na(race_clean) |
race_clean == "" |
str_detect(race_clean, "unknown|unspecified|undetermined|unreported|not available|race not") ~ NA_character_, # Disclosure: I used ChatGPT to create a list of entries in the data that correspond to NA for race
str_detect(race_clean, "white") &
!str_detect(race_clean, "hispanic|latino") ~ "White",
TRUE ~ "Non-White"
)
) |>
group_by(GEOID) |>
summarise(
victims_known = sum(!is.na(race_group)),
victims_nonwhite = sum(race_group == "Non-White", na.rm = TRUE),
.groups = "drop"
)
# Joining victim proportions to Census demographics
county_race_comparison <- countyFE |>
st_drop_geometry() |>
left_join(fe_race_counts, by = "GEOID") |>
mutate(
population_nonwhite_pct = 100 - white_pct,
victims_nonwhite_pct = if_else(
victims_known > 0,
100 * victims_nonwhite / victims_known,
NA_real_
)
)
# Preparing national and selected-state comparisons
race_plot <- bind_rows(
county_race_comparison |>
mutate(Group = "All U.S. Counties"),
county_race_comparison |>
filter(GEOID %in% selectedCountiesFE$GEOID) |>
mutate(Group = "Selected States")
) |>
filter(
pop >= 20000,
victims_known >= 10,
!is.na(population_nonwhite_pct),
!is.na(victims_nonwhite_pct)
)
# Comparing victim demographics with resident demographics
ggplot(
race_plot,
aes(x = population_nonwhite_pct, y = victims_nonwhite_pct)
) +
geom_point(alpha = 0.3, color = "steelblue") +
geom_abline(
slope = 1, intercept = 0,
linetype = "dashed", color = "gray40"
) +
geom_smooth(
method = "lm", se = TRUE,
color = "firebrick"
) +
ggpubr::stat_cor(
label.x.npc = 0.05,
label.y.npc = 0.95
) +
facet_wrap(~Group, ncol = 2) +
coord_equal(xlim = c(0, 100), ylim = c(0, 100)) +
labs(
title = "Racial Composition of Fatal Encounters vs. County Population",
subtitle = "Counties with 20,000+ residents and at least 10 victims with known race",
x = "Non-White Share of County Population (%)",
y = "Non-White Share of Recorded Victims (%)"
) +
theme_minimal()
I compared the non-White share of recorded fatal encounter victims with the non-White share of county populations. Both the national sample (R = 0.80) and selected states (R = 0.78) exhibit strong positive correlations. However, the fitted regression lines generally lie above the equality line, suggesting that non-White individuals constitute a disproportionately large share of recorded victims relative to their population share in many counties. These comparisons are limited by missing racial information and differences in racial classification between datasets.
# Calculate racial composition disparity by county
county_disparity <- countyFE |>
left_join(
county_race_comparison |>
select(GEOID, victims_known, victims_nonwhite_pct,
population_nonwhite_pct),
by = "GEOID"
) |>
mutate(
disparity = victims_nonwhite_pct - population_nonwhite_pct
) |>
filter(
victims_known >= 15,
!is.na(disparity)
)
# Map disparities across counties
tm_shape(county_disparity) +
tm_polygons(
fill = "disparity",
fill.scale = tm_scale_continuous(
values = "-RdBu",
midpoint = 0
),
fill.legend = tm_legend(
"Disparity (percentage points)"
),
col = "white",
lwd = 0.2,
popup.vars = c(
"County" = "NAME",
"Non-White Population (%)" = "population_nonwhite_pct",
"Non-White Victims (%)" = "victims_nonwhite_pct",
"Disparity (pp)" = "disparity",
"Victims with Known Race" = "victims_known"
)
) +
tm_layout(
title = "Racial Composition Disparities in Fatal Encounters",
legend.outside = TRUE
)The map shows percentage-point difference between the non-White share of recorded fatal encounter victims and the non-White share of county populations. Among counties with sufficient data, disparities appear predominantly skewed toward greater non-White representation among victims, with positive differences particularly concentrated around major urban areas. However, the geographic coverage is limited by minimum population and victim-count thresholds, and these differences do not account for variation in exposure to police encounters.
County Characteristics and Fatal Encounter Rates
# Preparing national data
national_plot <- countyFE |>
st_drop_geometry() |>
filter(pop >= 20000) |>
pivot_longer(
cols = c(poverty_pct, black_pct, white_pct, pop_density),
names_to = "variable",
values_to = "value"
) |>
mutate(
value = if_else(variable == "pop_density", log10(value), value)
)
# Preparing selected states data
selected_plot <- selectedCountiesFE |>
st_drop_geometry() |>
filter(pop >= 20000) |>
pivot_longer(
cols = c(poverty_pct, black_pct, white_pct, pop_density),
names_to = "variable",
values_to = "value"
) |>
mutate(
value = if_else(variable == "pop_density", log10(value), value)
)
# Faceted scatterplots
ggplot(selected_plot, aes(x = value, y = fe_per_10k)) +
geom_point(alpha = 0.3, color = "steelblue") +
geom_smooth(
data = national_plot,
method = "lm", se = FALSE,
color = "firebrick"
) +
geom_smooth(
method = "lm", se = FALSE,
color = "steelblue"
) +
stat_cor(
data = national_plot,
color = "firebrick",
label.x.npc = 0.05,
label.y.npc = 0.95,
color = "black",
r.accuracy = 0.001,
p.accuracy = 0.0001,
show.legend = FALSE,
size = 3
) +
stat_cor(
color = "steelblue",
label.x.npc = 0.05,
label.y.npc = 0.85,
color = "black",
r.accuracy = 0.001,
p.accuracy = 0.0001,
show.legend = FALSE,
size = 3
) +
facet_wrap(
~variable,
scales = "free_x",
ncol = 2,
labeller = as_labeller(c(
poverty_pct = "Poverty (%)",
black_pct = "Black Population (%)",
white_pct = "White Population (%)",
pop_density = "Log10 Population Density (people/km²)"
))
) +
labs(
title = "County Characteristics and Fatal Encounter Rates",
subtitle = "Blue: Selected States | Red: All U.S. Counties (Population ≥ 20,000)",
x = NULL,
y = "Fatal Encounters per 10,000 Residents"
) +
theme_minimal()
Next, I compared county-level fatal encounter rates against poverty, racial composition, and population density, both nationally and within the selected states. Black population share shows a positive correlation with encounter rates in both samples, while White population share and density exhibit negative correlations. Notably, poverty has a substantially stronger positive association in selected states (R = 0.327) than nationally (R = 0.09). Overall, the relationships remain relatively weak to moderate, suggesting that these characteristics alone cannot adequately explain geographic variation in fatal encounters. Between both selected states and national level data, spatial deviation does not allow for a definite conclusion on a correlation.
Regression Analysis
# National county-level regression
model_national <- lm(
fe_per_10k ~ poverty_pct + black_pct + white_pct + pop_density,
data = countyFE |>
st_drop_geometry() |>
filter(pop >= 20000, pop_density > 0)
)
# Selected states county-level regression
model_selected <- lm(
fe_per_10k ~ poverty_pct + black_pct + white_pct + pop_density,
data = selectedCountiesFE |>
st_drop_geometry() |>
filter(pop >= 20000, pop_density > 0)
)
model_comparison <- bind_rows(
tidy(model_national) |>
mutate(Model = "National"),
tidy(model_selected) |>
mutate(Model = "Selected States")
) |>
select(Model, term, estimate, p.value) |>
pivot_wider(
names_from = term,
values_from = c(estimate, p.value),
names_glue = "{.value}_{term}"
) |>
left_join(
bind_rows(
glance(model_national) |> mutate(Model = "National"),
glance(model_selected) |> mutate(Model = "Selected States")
) |>
select(
Model,
R_squared = r.squared,
Adjusted_R2 = adj.r.squared,
F_statistic = statistic,
Model_p = p.value
),
by = "Model"
) |>
mutate(
N = c(nobs(model_national), nobs(model_selected))
) |>
relocate(Model, N, R_squared, Adjusted_R2, F_statistic, Model_p)
# Disclosure: as I was running out of time, I used ChatGPT to finalize the visual presentation of this table to make it visually clear and readable:
model_table <- model_comparison |>
pivot_longer(
cols = -Model,
names_to = "Statistic",
values_to = "Value"
) |>
mutate(
Statistic = recode(
Statistic,
N = "Observations",
R_squared = "R²",
Adjusted_R2 = "Adjusted R²",
F_statistic = "F-statistic",
Model_p = "Model p-value",
`estimate_(Intercept)` = "Intercept",
estimate_poverty_pct = "Poverty (%)",
estimate_black_pct = "Black Population (%)",
estimate_white_pct = "White Population (%)",
estimate_pop_density = "Population Density",
`p.value_(Intercept)` = "Intercept p-value",
p.value_poverty_pct = "Poverty p-value",
p.value_black_pct = "Black Population p-value",
p.value_white_pct = "White Population p-value",
p.value_pop_density = "Population Density p-value"
),
Value = if_else(
str_detect(Statistic, "p-value"),
format.pval(Value, digits = 3, eps = 0.001),
format(round(Value, 4), nsmall = 2, trim = TRUE)
)
) |>
pivot_wider(
names_from = Model,
values_from = Value
) |>
mutate(
Statistic = factor(
Statistic,
levels = c(
"Observations", "R²", "Adjusted R²",
"F-statistic", "Model p-value",
"Intercept", "Intercept p-value",
"Poverty (%)", "Poverty p-value",
"Black Population (%)", "Black Population p-value",
"White Population (%)", "White Population p-value",
"Population Density", "Population Density p-value"
)
)
) |>
arrange(Statistic)
model_table |>
kbl(
caption = "Comparison of County-Level Regression Models",
align = c("l", "r", "r"),
col.names = c("Statistic", "National", "Selected States")
) |>
kable_styling(
bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE,
position = "center"
) |>
pack_rows("Model Fit", 1, 5) |>
pack_rows("Regression Coefficients", 6, 15)| Statistic | National | Selected States |
|---|---|---|
| Model Fit | ||
| Observations | 1879.0000 | 431.0000 |
| R² | 0.0994 | 0.1778 |
| Adjusted R² | 0.0974 | 0.1700 |
| F-statistic | 51.6919 | 23.0244 |
| Model p-value | < 0.001 | < 0.001 |
| Regression Coefficients | ||
| Intercept | 0.6290 | 1.7063 |
| Intercept p-value | < 0.001 | < 0.001 |
| Poverty (%) | 0.0041 | 0.0329 |
| Poverty p-value | 9.83e-02 | < 0.001 |
| Black Population (%) | 0.0163 | -0.0055 |
| Black Population p-value | < 0.001 | 2.10e-02 |
| White Population (%) | 0.0008 | -0.0136 |
| White Population p-value | 3.40e-01 | < 0.001 |
| Population Density | -0.0001 | 0.0000 |
| Population Density p-value | 2.43e-03 | 9.90e-01 |
Finally, I constructed two multiple linear regression models to examine how county poverty rates, Black and White population shares, and population density relate to fatal encounters per 10,000 residents. Both models were statistically significant but explained relatively little variation (adjusted 9.7% nationally and 17% in selected states). Poverty exhibited a stronger positive association in selected states, while racial composition coefficients differed in direction and magnitude between the two samples. These differences highlight the importance of geographic context and suggest that relationships observed nationally may not hold uniformly across regions. However, substantial unexplained variation remains, and the models do not account for factors such as data quality and policing practices.
Conclusion
This analysis identifies a substantial increase in recorded fatal police encounters in the United States between 2000 and 2021. While a similar upward trend appears in Mapping Police Violence data, improvements in reporting and data collection may contribute to the observed increase. Hence, it is difficult to determine the precise magnitude of increase.
The analysis also reveals associations between fatal encounter rates and socioeconomic and racial/ethnic characteristics, although these relationships vary considerably across geographic contexts. The regression models explain relatively little variation, suggesting that national-level patterns alone are insufficient in explaining police-related fatalities.
Comparisons of victim demographics with county populations indicate a disproportionate representation of non-White individuals among recorded victims in many counties, especially around major cities. The correlations must also be cross-checked with other social, cultural, and historical contexts.
An important limitation is the use of fatal encounters as an indicator of police violence. Fatalities represent one of the most extreme outcomes of police interactions and most likely does not capture broader patterns of excessive force, discriminatory enforcement, or unequal treatment.
With additional time, I would have examined the geographic distribution of fatal encounters within urban and rural environments, particularly in states containing major metropolitan areas. I would have also tried regression models with adjusted variables for different places. The primary finding is that recorded fatal police encounters exhibit substantial geographic and racial disparities, but their underlying causes remain difficult to establish through nation-level analysis.