Gun Violence Report from 2013-2018

Load the libraries

library(tidyr)
library(ggplot2)
library(skimr)
library(readr)
library(dplyr)
library(stringr)
library(lubridate)
library(gt)
library(gtExtras)
library(plotly)
library(maps)
library(gridExtra)
library(archive)
library(readr)
library(shiny)
library(rmarkdown)
library(vcd)
library(ggmap)
library(wordcloud2)
library(tidytext)
library(boot)
library(naniar)

Load the dataset

crime_stat <- read_csv(
  archive_read("DATA_01-2013_03-2018.tar.gz", file = 1), 
  col_types = cols())

Abstract

Here we have an extract of the gun violence related incidents data set from March 2013 to March 2018. It records each incident by a primary key called incident id. Most of the data has been scraped from website therefore there are records that are missing, when we perform the analysis we will be removing all incidents that have missing values. We will also split the some of the fields that are delimited by || and store in a key value format separated by ::.

The purpose of the analysis is to get a good understanding of the crime that occurs in each state and see which areas have more crime vs which areas have less. We also want to answer the following questions when it comes to crime related activity.

1. How have gun-related incidents changed over time? Is there a seasonal pattern?
2. Which states or cities have the highest number of incidents?
3. What is the age distribution of the participants involved in these incidents? How does age correlate with the role (victim, suspect, etc.)?
4. Analysis of notes and descriptions provided in the data set to understand the context or background of incidents.
5. Is there a correlation between the number of incidents and congressional districts?

Data Dictionary

Here we have the data dictionary, specifying the field name, data-type, brief description and if the field is required.

Field	Data Type	Description	Required?
incident_id	int	gunviolencearchive.org ID for incident	Yes
date	str	date of occurrence	Yes
state	str	State in which the incident occurred	Yes
city_or_county	str	City or county in which the incident occurred	Yes
address	str	address where incident took place	Yes
n_killed	int	number of people killed	Yes
n_injured	int	number of people injured	Yes
incident_url	str	link to gunviolencearchive.org webpage containing details of incident	Yes
source_url	str	link to online news story concerning incident	No
incident_url_fields_missing	bool	ignore, always False	Yes
congressional_district	int	Number of districts	No
gun_stolen	dict[int, str]	key: gun ID, value: ‘Unknown’ or ‘Stolen’	No
gun_type	dict[int, str]	key: gun ID, value: description of gun type	No
incident_characteristics	list[str]	list of incident characteristics	No
latitude	float	geolocation	No
location_description	str	description of location where incident took place	No
longitude	float	geolocation	No
n_guns_involved	int	number of guns involved	No
notes	str	additional notes about the incident	No
participant_age	dict[int, int]	key: participant ID	No
participant_age_group	dict[int, str]	key: participant ID, value: description of age group, e.g. ‘Adult 18+’	No
participant_gender	dict[int, str]	key: participant ID, value: ‘Male’ or ‘Female’	No
participant_name	dict[int, str]	key: participant ID	No
participant_relationship	dict[int, str]	key: participant ID, value: relationship of participant to others	No
participant_status	dict[int, str]	key: participant ID, value: ‘Arrested’, ‘Killed’, ‘Injured’, or ‘Unharmed’	No
participant_type	dict[int, str]	key: participant ID, value: ‘Victim’ or ‘Subject-Suspect’	No
sources	list[str]	links to online news stories concerning incident	No
state_house_district	int	Number of state house districts	No
state_senate_district	int	Number of state senate districts	No

Missingness of data

gg_miss_var(crime_stat, show_pct = T)

According to the graph we do have some missing data present in our dataset. The following variables seem to have the most missing data:

1. participant_relatinship
2. location_description
3. participant_name
4. gun_stolen
5. n_gun_involved

Even though we have some missing records in our dataset. It will not really impact our analysis as we may not even use some of the data. Therefore lets move onto the data cleansing portion!

Data Cleaning

Some of the columns are concatenated by || and are stored in a key-value pair format. Lets create a custom function that runs on the concatenated columns in order to perform some meaningful analysis.

Generate sub tables from the current data set.

columns_to_process = c("gun_stolen", "gun_type", "participant_age",
                       "participant_age_group", "participant_gender",
                       "participant_name", "participant_relationship",
                       "participant_status", "participant_type")

split_tables_func <- function(df, column_name_list) {
  # Initialize a list to store the data frames
  list_of_dfs <- list()

  # Check if all column names in the list are actually present in the data frame
  valid_columns <- column_name_list[column_name_list %in% names(df)]

  for (i in valid_columns) {
    # Process each column and store in the list
    list_of_dfs[[i]] <- df |>
      select(incident_id, all_of(i)) |>
      separate_rows(all_of(i), sep = "\\|\\|") |>
      separate(
        all_of(i),
        into = c("person_id",
                 paste0(i, "_status")),
        sep = "::",
        convert = TRUE) # Ensure person_id is an integer
  }

  # Return the list of data frames
  return(list_of_dfs)
}

Now that we have the function defined and the list of columns to split defined we will go ahead and create the different tables to perform specific analysis. The idea is to split the columns into different tables using a primary key “incident_id” and a surrogate key “person_id” to correctly merge data-sets.

# run the function to split all columns into sub tables
processed_dfs <- split_tables_func(crime_stat,columns_to_process)

# Taking a look at the output of one of the generated tables.
head(processed_dfs$participant_age) |>
  gt() |>
  gt_theme_guardian() |>
  tab_header(title = "First 6 records from Participant Age Table")

First 6 records from Participant Age Table
incident_id	person_id	participant_age_status
461105	0	20
460726	0	20
478855	0	25
478855	1	31
478855	2	33
478855	3	34

# Once the list of data frames is generated we can split the list
# by using the list2env function
list2env(processed_dfs, envir = .GlobalEnv)

## <environment: R_GlobalEnv>

Now that our data pre-processing is completed and the data-set looks ready for analysis lets try to answer a few questions to get better insights of the data.

Exploratory Data Analysis

1. How have gun-related incidents changed over time? Is there a seasonal pattern?

Plot 1: Annual Trend of Incidents - Plot the number of incidents over years.

We want to see the seasonality of the gun violence occurring over the years. This can give us a good idea of when gun violence was the highest vs lowest from 2013-2018

# Extract the year using the lubridate package. 
monthly_trend <- crime_stat |>
  group_by(date) %>%
  summarise(incident_count = n()) %>%
  ungroup()

monthly_trend <- monthly_trend |>
  mutate(year = year(date))

# Get the list of all unique years in the data-set.
years <- unique(monthly_trend$year)

# Create the list to make the drop down buttons. Each button is a year filter.
year_buttons <- lapply(years, function(y){
  list(
    method = "restyle",
    args = list("transforms[0].value", y),
    label = as.character(y)
  )
})

fig <- plot_ly(monthly_trend,
               x = ~date,
               y = ~incident_count,
               type = 'scatter',
               mode = 'lines+markers',
               line = list(color='steelblue'), marker = list(color='red'),
               transforms = list(
                 list(
                   type = 'filter',
                   target = ~year,
                   operation = '=',
                   value = years[1]
                 )
               )) |>
                 layout(title='Trend of Incidents',
                        xaxis= list(title='Date'),
                        yaxis= list(title='Number of Indcidents'),
                        updatemenus = list(
                          list(
                            y = 0.8, 
                            buttons = year_buttons
                            )
                          )
                        )
fig

When we look at the count of incidents over the year, we notice some patterns in the data. Once obvious point is that crime is consistently lower during the holidays. If you filter the chart by selecting the year drop-down you can notice that crime is always the lowest during and around Christmas.

Another note able trend we can notice is crime becomes an all time high during the month of August 2016. Based on wikipedia this was around the time of the George Floyd incident, 2016 Louisiana Floods that miss-placed people from 146,000 homes and also the California fires that displaced over 82,000 residents. Another massive uptick in crime is seen on Jan 1st of 2017 right about the time President Trump took office. This date has the highest incident reports. This trends need to be further analyzed and look at by state to get a good understanding of what incidents exactly caused the huge rise in crime at the end of 2016.

2. Which states or cities have the highest number of incidents?

summary statistics

We would first like to see the summary statistics of the number of incidents, people killed and people injured by state.

summary_by_state <- crime_stat |>
  group_by(state) |>
  summarise(
      incidents_count = n(),
      avg_killed = round(mean(n_killed, na.rm = T), 2),
      std_killed = round(sd(n_killed, na.rm = T), 2),
      avg_injured = round(mean(n_injured, na.rm = T), 2),
      std_injured = round(sd(n_injured, na.rm = T), 2)
            ) |>
  arrange(desc(incidents_count)) |>
  head()|>
  gt()|>
  cols_label(
    incidents_count = "Total Count",
    avg_killed = "Mean Deaths",
    std_killed = "SD Deaths",
    avg_injured = "Mean Injured",
    std_injured = "SD Injured") |>
  gt_theme_guardian() |>
  tab_header(title = "Count, Mean and SD of killed & Injured by State") |>
  tab_caption(caption = md("Table1: Summary Statistics of the states in USA"))

summary_by_state

Table1: Summary Statistics of the states in USA

Count, Mean and SD of killed & Injured by State
state	Total Count	Mean Deaths	SD Deaths	Mean Injured	SD Injured
Illinois	17556	0.19	0.44	0.77	0.76
California	16306	0.34	0.57	0.47	0.79
Florida	15029	0.26	0.68	0.47	0.87
Texas	13577	0.37	0.64	0.45	0.75
Ohio	10244	0.24	0.50	0.56	0.72
New York	9712	0.18	0.43	0.52	0.73

Looking at the summary statistic something interesting stands out to me if you look at the number of incidents, Illinois has the highest count but notice the means and standard deviations of killed vs injured. The counts of people being killed is significantly lower as compared to the next largest which is California. They do however have the largest injuries out of all the other incidents. Taking a look at the breakdown by city_or_county field (Not shown) maybe it could just be that Chicago has better policing laws and there is more petty theft in the city rather than the other parts of the state. Similarly California has a similar trend of having more people injured in the city rather than the outskirts where there is a higher incidents of death by gun violence. Lets build a bar graph to represent this and further drill down if that is really the case.

Plot 2: Gun Violence per state and city.

# crime_stat <- crime_stat |>
#   arrange(state)
# 
# ui <- fluidPage(
#   titlePanel("Top 20 Cities by State & Incident Counts"),
#   selectInput("stateInput", "Choose a State:",
#               choices = unique(crime_stat$state)),
#   plotlyOutput("incidentPlot")
# )
# 
# server <- function(input, output) {
#   # Reactive expression to filter the data-set based on
#   # selected state and calculate top 20 cities
#   filtered_data <- reactive({
#     crime_stat |>
#       filter(state == input$stateInput) |>
#       group_by(city_or_county) |>
#       summarise(incidents = n(), .groups = 'drop') |>
#       slice_max(n = 20, order_by = incidents) |>
#       arrange(desc(incidents))
#   })
# 
#   # Render the plot using Plotly
#   output$incidentPlot <- renderPlotly({
#     fig <- plot_ly(data = filtered_data(),
#                    y = ~city_or_county,
#                    x = ~incidents,
#                    type = 'bar',
#                    orientation = 'h',
#                    marker = list(color = 'steelblue')) |>
#       layout(title = paste("Top 20 Cities in", input$stateInput, "by Incident Counts"),
#              xaxis = list(title = "Number of Incidents"),
#              yaxis = list(title = "City or County", automargin = TRUE)) |>
#       layout(xaxis = list(categoryorder = "total descending"))
# 
#     fig
#   })
# }
# 
# # Run the application
# shinyApp(ui = ui, server = server)

(adding Screenshots of the output as I was not able to get shiny to work on rpubs)

Graph of Alabama

Graph of California

Graph of Michigan

When we filter down the graphs based on the specific state it seems evident that there are a significant number of incidents in Urban areas rather than in the small towns and rural areas.

Lets also look at the geographical heat-map of areas that had the most number of incidents from 2013 to 2018 based on the urban and rural areas.

Plot 3: Heat Map of crime in Rural vs Urban Areas.

# Lets merge a dataset from https://worldpopulationreview.com/states/state-abbreviations
# to get the matching state names.

states <- read_csv("table-data.csv")

## Rows: 51 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): state, abbrev, code
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

location_smmry <- crime_stat |>
  mutate(year = year(crime_stat$date))|>
  group_by(state, year) |>
  summarise(total_incidents = n(),
            total_killed = sum(n_killed),
            total_injured = sum(n_injured)) |>
  ungroup() |>
  left_join(states, by = 'state') |> #merge the tables
  mutate(hover = paste0("Total Incidents in ", state,":", total_incidents,
                        "\n", "Total Killed: ", total_killed,
                        "\n", "Total Injured: ", total_injured))

## `summarise()` has grouped output by 'state'. You can override using the
## `.groups` argument.

#Lets plot the table using plotly
heatmap <- plot_geo(location_smmry,
                    locationmode = 'USA-states',
                    frame = ~year) |>
  add_trace(locations = ~code,
            z = ~total_incidents,
            zmin = 0,
            zmax = max(location_smmry$total_incidents),
            color = ~total_incidents,
            colorscale = 'Reds',
            text = ~hover,
            hoverinfo = 'text') |>
  layout(geo = list(
    scope = 'usa'
  ),
  title = list(text = 'Choropleth Heat map of incidents from 2013 to 2018',
  y = 0.2))

heatmap

based on this graph we do see a significant increase in the crime rates by state. Lets further drill down and see the exact location of where these incidents are taking place. Is it mostly in Urban areas or rural areas?

incident_loc <- crime_stat |>
  mutate(year = year(date)) |>
  filter(n_killed > 0) |>
  select(year, longitude, latitude, n_killed, n_injured) |>
  group_by(year)
  


geo_properties <- list(
  scope = 'usa',
  projection = list(type = 'albers usa'),
  showland = T,
  showsubunits = F,
  landcolor = toRGB('gray10'),
  showlakes = T,
  lakecolor = toRGB('white')
)

urban_rural <- plot_geo(incident_loc,
                        lat = ~latitude,
                        lon = ~longitude,
                        frame = ~year,
                        marker = list(
                          size = 2,
                          color = "#FFFFCC",
                          opacity = 0.25
                        )) |>
  add_markers(hoverinfo = "none")|>
  layout(geo = geo_properties,
         title = list(
           text = 'Choropleth Heat map of incidents in Rural vs Urban areas.',
           y = 0.2))

urban_rural

## Warning: Ignoring 3004 observations

Based on this heatmap it is evident that there definitly is more activity and incidents in the Urban areas rather than the rural areas. The high desnity clusters are all major cities such as Chicago, Atlanta, Detroit and Los Angeles.

3. What is the age distribution of the participants involved in these incidents? How does age correlate with the role (victim, suspect, etc.)?

summary statistics

In order to properly answer this question lets take a look at some summary statistics of gun related incidents by age groups.

#select the specified columns from the table
age_grp_smmry <- crime_stat |>
  select(incident_id, n_killed, n_injured)

age_grp_smmry <- merge(age_grp_smmry, participant_age_group, 
                       by=c("incident_id"), all = TRUE)
  
age_grp_smmry <- merge(age_grp_smmry, participant_type,  
                       by=c('incident_id', 'person_id'), all = T)

age_grp_smmry <- merge(age_grp_smmry, participant_status,  
                       by=c('incident_id', 'person_id'), all = T)

age_grp_smmry <- merge(age_grp_smmry, participant_gender,  
                       by=c('incident_id', 'person_id'), all = T)

age_grp_smmry <- merge(age_grp_smmry, participant_age,  
                       by=c('incident_id', 'person_id'), all = T)

summary_stats <- age_grp_smmry |>
  filter(participant_type_status != 'Victim' & 
           participant_age_group_status != "NA" & 
           participant_status_status == 'Killed' & 
           participant_gender_status !="NA") |>
  group_by(participant_gender_status,
           participant_age_group_status) |>
  summarise(
    total_count = n(),
    mean_killed = round(mean(n_killed, na.rm = T), 2),
    sd_killed = round(sd(n_killed, na.rm = T), 2)
  ) |>
  arrange(desc(total_count)) |>
  gt() |>
  cols_label(participant_age_group_status = "Age Group") |>
  gt_theme_guardian() |>
  tab_header(title = "Count, Mean and SD of people killed in USA by Age Group") |>
  tab_caption(caption = md("Table2: Summary Statistics of the age group
                           distribution by male and female."))

## `summarise()` has grouped output by 'participant_gender_status'. You can
## override using the `.groups` argument.

summary_stats

Table2: Summary Statistics of the age group distribution by male and female.

Count, Mean and SD of people killed in USA by Age Group
Age Group	total_count	mean_killed	sd_killed
Male
Adult 18+	7953	1.38	0.93
Teen 12-17	171	1.24	0.62
Child 0-11	3	1.00	0.00
Female
Adult 18+	393	1.57	1.05
Teen 12-17	9	1.44	0.53
Child 0-11	2	2.00	1.41

Based on the table above we notice that the on average Males in the Adult 18+ age group are killed due to gun related violence. Followed by teenagers and lastly the children ages 0-11. Per incident we see that there is an average of at least 1 person being killed.

Females in the data set are quite less in comparison to the men but we do see on average per incident they have slightly more deaths. This figure may just be high due to the fact that we have a lower count of women in the data-set.

4. Analysis of notes and descriptions provided in the data set to understand the context or background of incidents.

We want to check the patterns and specific situations that lead to incidents. We would also like to see the commonality what trends are mostly reported when an incident occurs.

Lets look at the variables notes and the incident_characteristics column and create a word map.

wordFreqs <- crime_stat |>
  unnest_tokens(word, incident_characteristics) |>
  filter(!grepl('[0-9]', word)) |>
  count(word)

wordFreqs2 <- crime_stat |>
  unnest_tokens(word, notes) |>
  filter(!grepl('[0-9]', word)) |>
  count(word)

wordFreqs <- merge(wordFreqs, wordFreqs2, by = "word")

# Remove the stop words
wordFreqs <- wordFreqs |>
  anti_join(stop_words)

## Joining with `by = join_by(word)`

wordcloud2(wordFreqs, size =1.6,
           color = 'random-dark')

When observing these word clouds, there does seem to be somewhat of a pattern. One obvious point we can note is most of the gun violence incidents are done by accident. Typically it involves car, drugs and robbery.

Lets look at the next analysis!

5.Explore if there’s a correlation between the number of incidents and congressional districts.

Plot5: Scatter plot comparing incidents vs congressional districts.

Here we would like to see the correlation between the number of incidents vs the congressional districts. Which ruling party in a particular state is able to manage gun violence based on their particular states gun laws.

# Load the statewide senate data from #https://dataverse.harvard.edu/file.xhtml?fileId=7609736&version=7.0&toolType=PREVIEW

# We want to use this data to mark the states that are belonging to different
# parties to see if crime is more due to the congressional district being
# Democratic, Republic or Other.

state_poli <- read.csv("1976-2020-senate.csv", header = TRUE) |>
  filter(year == 2018) |>
  select(state, party_simplified, candidatevotes) |>
  group_by(state) |>
  top_n(1, candidatevotes)

state_poli$state <- tolower(state_poli$state)

district_summary <- crime_stat |>
  group_by(congressional_district, state) |>
  summarise(incidents = n())

## `summarise()` has grouped output by 'congressional_district'. You can override
## using the `.groups` argument.

district_summary$state <- tolower(district_summary$state)

# Merge the two tables to get the matching political party
combined <- left_join(district_summary, state_poli, by = 'state')

# Replace NA values in the party column with "other"
combined <- combined %>%
  mutate(party_simplified = if_else(is.na(party_simplified), "OTHER", party_simplified))

# Group the dataset
combined <- combined |>
  select(congressional_district, state, party_simplified, incidents)|>
  group_by(state, party_simplified) |>
  summarise(
    total_districts = sum(congressional_district, na.rm = TRUE),
    total_incidents = sum(incidents, na.rm = FALSE)
  )

## `summarise()` has grouped output by 'state'. You can override using the
## `.groups` argument.

# Ensure your 'party_simplified' column is factored with the correct order
combined$party_simplified <- factor(combined$party_simplified, levels = c("DEMOCRAT", "REPUBLICAN", "OTHER"))

Now that we have merged and combined our data sets, lets plot in a scatter graph.

# Create a Plotly scatter plot
fig <- plot_ly(data = combined,
               x = ~total_incidents, 
               y = ~total_districts,
               color = ~party_simplified,
               hoverinfo = "text",
               text = paste("State: ", combined$state,
                            "<br>",
                            "Districts: ", combined$total_districts,
                            "<br>",
                            "Incidents: ", combined$total_incidents)) |>
  add_markers(
    marker = list(
      colors = ~party_simplified)) |>
  layout(
    title = "Scatter plot comparing incidents vs total_districts",
    xaxis = list(title="Total Count of Incidents"),
    yaxis = list(title="Total Count of Districts"),
    hoverlabel = list(
      font = list(size = 10)))

# Print the plot
fig

Based on the scatter plot it seems evident that democratic states on average have more incidents of gun violence and crime. The highest being California. The lowest being Hawaii (Zoom into the graph.) It almost seems like there is somewhat of an inverse correlation between the gun law strictness and incidents reported. Since most blue states have stricter gun laws there seem to be more violent crimes.

Monte Carlo

Geographic Hotspots Analysis: Use Monte Carlo to help identify which areas (states, cities) might become hotspots for future incidents based on historical patterns.

Null Hypotheses: (H0): There is no difference in the distribution of incidents among the different geographic areas Alternate Hypothesis: (H1): Some areas have higher frequency of incidents than others.

# lets prepare the data for analysis

incident_summary <- crime_stat |>
  group_by(state) |>
  summarise(total_incidents = n()) |>
  # Finding the likelihood of incidents.
  mutate(probs = total_incidents / sum(total_incidents))

# Now that we have the probabilities by each state lets run the simulation
set.seed(123) #setting seed for reproducibility
n_simulations <- 1000
# Initialize the dataframe
simulated_data <- data.frame(state = character(0), count = integer(0))

for(i in 1:n_simulations) {
  simulated_states <- sample(incident_summary$state,
                             size = 100, # number of incidents per sim
                             replace = TRUE,
                             prob = incident_summary$probs)
  simulation_count <- as.data.frame(table(simulated_states))
  colnames(simulation_count) <- c("state", "count")
  simulation_count$simulation <- i
  simulated_data <- rbind(simulated_data, simulation_count)
}

# Aggregate the simulation results.
simulation_results <- simulated_data |>
  group_by(state) |>
  summarise(average_count = mean(count))

Now that we ran the simulation lets plot the results in graph to see potential hot-spots for future crime incidents based on the historical data.

simulation_results |>
  arrange(desc(average_count))|>
  plot_ly(
    x = ~state,
    y = ~average_count
  ) |>
  add_bars() |>
  layout(xaxis = list(categoryorder = "total descending"))

# ggplot(simulation_results, aes(x = reorder(state, average_count), y = average_count)) +
#   geom_bar(stat = "identity") +
#   coord_flip()+
#   labs(title = "Average Number of Incidents by State in Simulated Data",
#        x = "State",
#        y = "Average Number of Incidents")+
#   theme_bw()
# 
# options(repr.plot.width = 10, repr.plot.height = 25)

Based on the historical pattern it seems that Illinois is on the rise in terms of being a hot spot for crime followed by California and then Florida. Nothing seems unexpected as per the state but further analysis and drill-down might be needed by city and county to see the whole picture.

Bootstrapping.

Influence of Gun Availability: Explore how the n_guns_involved and gun_type variables influence the outcomes of incidents.

# Merge the split tables into one big table for analysis.
bts_data <- crime_stat |>
  select(incident_id, n_killed, n_guns_involved) |>
  filter(n_killed > 0) |>
  group_by(incident_id) |>
  left_join(gun_type, by = c("incident_id"))

# Define a statistic function to explore the influence of gun availability
effect_of_guns <- function(data, indices) {
  sampled_data <- data[indices, ] #Re-sampling with replacement
  lm_model <- lm(n_killed ~ n_guns_involved, data = sampled_data)
  return(coef(lm_model)) #Return the coefficients in a linear model.
}

#Perform the bootstrap
bootstrap_results <- boot(bts_data, statistic = effect_of_guns, R = 1000)

#View the results
print(bootstrap_results)

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = bts_data, statistic = effect_of_guns, R = 1000)
## 
## 
## Bootstrap Statistics :
##        original        bias     std. error
## t1* 1.145691356 -1.400447e-04 0.0037198889
## t2* 0.005177562  1.425358e-05 0.0007491153

Based on the bootstrap results our model looks like the following:

n_killed = Beta0 + Beta1 * n_guns_involved + beta2 * gun_type + e

The original estimates t1 and t1 are 1.145 and 0.005 respectively. The first value represents the number of people killed when no guns are involved. This is our baseline. The second value represents the change in expected number of people killed for each additional gun involved. The positive value does suggest that number of guns involved does increase the incidents and fatalities.

The small biases of 8.16e-05 and 1.31e-05 suggest that the estimates are close to the original estimate which is reliable for our analysis.

Lastly, the standard error being low does suggest that the intercept estimate is precise in accuracy.

plot(bootstrap_results)

Our qq-plot is perfectly in line suggesting that the data is normally distributed. So it the histogram.

#Calculate the confidence intervals for the coefficients
conf_intervals <- boot.ci(bootstrap_results, type = "norm")

conf_intervals

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = bootstrap_results, type = "norm")
## 
## Intervals : 
## Level      Normal        
## 95%   ( 1.139,  1.153 )  
## Calculations and Intervals on Original Scale

#graph the coefficients. 
boot_coefs <- as.data.frame(bootstrap_results$t)
colnames(boot_coefs) <- c("Intercept", "Slope")

ggplot(boot_coefs, aes(x = Slope)) +
  geom_histogram(binwidth = 0.0001, fill = "blue", color = "black") +
  geom_vline(xintercept = 0.005177562, color = "red", linetype = "dashed", size = 1) +
  labs(title = "Bootstrap Distribution of Slope (Effect of Number of Guns)",
       x = "Slope Coefficient",
       y = "Frequency")

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

According to the Confidence interval we are 95% confident that the true average number of fatalities when no guns are involved lies between 1.138 and 1.53. The interval is very narrow, suggesting a high level of precision in the estimate of the baseline number of fatalities.