| Data Set Variables |
|---|
| Weather |
| Weather Frequency by Damage |
| Drivers License plate |
| Speed Limit |
| Vehicle Make |
| Latitude |
| Longitude |
| Year (of crash) |
| Damage |
| Damage Frequency |
(Total 39 variables; I’ve included the ones I use in this analysis)
For this analysis, I selected a data set from Data.gov that reports crash data in Montgomery County, Maryland. The data set, collected through the Automated Crash Reporting System (ACRS) managed by the Maryland State Police, includes detailed records of vehicle crashes categorized into three types: fatal, property damage, and injury-related crashes.
The data is compiled from reports filled out by law enforcement officers at the scene of a crash. These reports include a variety of factors such as driver and vehicle information, weather and road conditions, and contributing causes.
In this research, I aim to explore questions such as: “Does weather play a role in the frequency or type of crashes?” and “Is there a relationship between the make of a vehicle—like Toyota, Mercedes, or Tesla—and how often it’s involved in a crash?” These questions are not only data-driven but also reflect real-world concerns about safety and accountability on the road.
library(tidyverse)
library(viridis)
library(ggplot2)
library(DataExplorer)Warning: package 'DataExplorer' was built under R version 4.5.2library(leaflet)
library(sf)Warning: package 'sf' was built under R version 4.5.2library(plotly)Warning: package 'plotly' was built under R version 4.5.2crash <- read_csv("/Users/maishasubin/Desktop/DATA110/Crash_Reporting_Data.csv")# Removing columns not used in this analysis
crash_col <- crash[, -c(1, 2, 3, 7, 8, 9, 10, 11, 18, 19, 25, 26, 27, 28, 30, 32, 33, 34, 36 )] # Create separate columns for date and time
crash_time <- crash_col %>%
separate('Crash Date/Time', into = c("year", "time"), sep = " ")
# Extract year as numeric from date column
crash_time$year <- as.numeric(str_sub(crash_time$year, start = -4L))
# Sort the years in ascending order
sorted_years <- sort(crash_time$year)
unique(sorted_years) [1] 2015 2016 2017 2018 2019 2020 2021 2022 2023 2024 2025# Renaming variables
crash_rename <- crash_time |>
rename(damage = `ACRS Report Type`
, route_type = `Route Type`
, collision_type = `Collision Type`
, surface_condition = `Surface Condition`
, weather = `Weather`
, light = `Light`
, traffic_control = `Traffic Control`
, substance_abuse = `Driver Substance Abuse`
, driver_fault = `Driver At Fault`
, injury_severity = `Injury Severity`
, circumstance = `Circumstance`, distraction = `Driver Distracted By`
, license_state =`Drivers License State`
, vehicle_movement = `Vehicle Movement`
, speed_limit = `Speed Limit`
, vehicle_make = `Vehicle Make`
, latitude = `Latitude`
, longitude = `Longitude`
, location = `Location`)
# Removing missing values: NA's and manually string "N/A"
crash_nona <- crash_rename %>%
filter(!is.na(damage) & damage != "N/A", !is.na(year) & year != "N/A"
, !is.na(time) & time != "N/A"
, !is.na(route_type) & route_type != "N/A"
, !is.na(collision_type) & collision_type != "N/A"
, !is.na(surface_condition) & surface_condition != "N/A"
, !is.na(weather) & weather != "N/A"
, !is.na(light) & light!= "N/A"
, !is.na(traffic_control) & traffic_control != "N/A"
, !is.na(substance_abuse) & substance_abuse != "N/A"
, !is.na(driver_fault) & driver_fault != "N/A"
, !is.na(injury_severity) & injury_severity != "N/A"
, !is.na(circumstance) & circumstance != "N/A"
, !is.na(license_state) & license_state != "N/A"
, !is.na(vehicle_movement) & vehicle_movement != "N/A"
, !is.na(speed_limit) & speed_limit != "N/A"
, !is.na(distraction) & distraction != "N/A"
, !is.na(vehicle_make) & vehicle_make != "N/A"
, !is.na(latitude) & latitude != "N/A"
, !is.na(longitude) & longitude != "N/A"
, !is.na(location) & location != "N/A")# Removing additional columns
crash_nona <- crash_nona %>%
select(-starts_with("weather_freq_by_damage"))
# Counting how often each weather type appears within each damage type
weather_freq_by_damage <- crash_nona %>%
count(damage, weather, name = "weather_freq_by_damage") # Combining 2 variable's
# Joining those frequencies back into my original data
crash_nona <- crash_nona %>%
left_join(weather_freq_by_damage, by = c("damage", "weather"))# Creating a frequency table for the 'damage' column
damage_freq <- table(crash_nona$damage)
# Mapping frequencies back to the original data frame
crash_nona$damage_freq <- damage_freq[crash_nona$damage]# Creating frequency table for 'traffic_control'
traffic_control_freq <- table(crash_nona$traffic_control)
# Mapping frequency values back to the data frame
crash_nona$traffic_control_freq <- traffic_control_freq[crash_nona$traffic_control]# To analyze correlation of selected variables
crash_corr <- crash_nona [, c("speed_limit", "year", "latitude", "longitude", "weather_freq_by_damage", "damage_freq", "traffic_control_freq")]
plot_correlation(crash_corr)Warning in dummify(data, maxcat = maxcat): Ignored all discrete features since
`maxcat` set to 20 categories!
# Running linear regression models
lm_model <- lm(damage_freq ~ speed_limit + year
+ latitude + longitude + weather_freq_by_damage
+ traffic_control_freq, data = crash_nona)
# Summary of model
summary(lm_model)
Call:
lm(formula = damage_freq ~ speed_limit + year + latitude + longitude +
weather_freq_by_damage + traffic_control_freq, data = crash_nona)
Residuals:
Min 1Q Median 3Q Max
-12412 -2593 1080 2029 3784
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.327e+05 1.640e+04 -8.091 6.18e-16 ***
speed_limit -1.631e+01 1.943e+00 -8.395 < 2e-16 ***
year 8.295e+01 5.400e+00 15.360 < 2e-16 ***
latitude -4.620e+02 2.565e+02 -1.801 0.0717 .
longitude 5.605e+01 1.941e+02 0.289 0.7728
weather_freq_by_damage 2.770e-01 6.524e-03 42.460 < 2e-16 ***
traffic_control_freq -3.702e-03 5.221e-03 -0.709 0.4783
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2381 on 24855 degrees of freedom
Multiple R-squared: 0.07452, Adjusted R-squared: 0.0743
F-statistic: 333.6 on 6 and 24855 DF, p-value: < 2.2e-16Equation for my linear regression model:
damage_freq = − 132700 − (16.31 * speed_limit) + (82.95 * year) − (462 * latitude) + (56.05 * longitude)− 8.550 + (0.277 * weather_freq_by_damage) − (0.0037 * traffic_control_freq)
Analyzing model based on p-values:
speed_limit, year and weather_freq_by_damage type are statistically significant predictors of crash damage frequency with p values of < 2e-16, < 2e-16 and < 2e-16 respectively.
Variables like latitude, longitude, weather_freq_by_damage, and traffic_control_freq are not statistically significant based on their high p-values (around or greater than 0.05).
Adjusted R^2 values:
Adjusted R^2 = 0.0743 suggests this model explains about 7.4% of the variability in the response variable damage_freq, even after adjusting for the number of predictors. While the model is statistically significant (F-statistic p < 2.2e-16), it’s not strong, the predictive power is low.
# Grouping data by year and damage
crash_yearly <- crash_nona %>%
group_by(year, damage) %>%
summarise(damage_count = n(), .groups = "drop")
# Converting year to factor for discrete x-axis
crash_yearly$year <- as.factor(crash_yearly$year)
# Facet plot
ggplot(crash_yearly, aes(x = year, y = damage_count, color = damage, group = damage)) +
geom_point() +
geom_line() +
facet_wrap(~damage) +
labs(title = "Crash Count by Year and Damage Type",
x = "Year", y = "Number of Crashes",
caption = "Source: Data.gov") +
scale_color_viridis_d("Damage Type") +
theme_bw()
library(leaflet.extras) Warning: package 'leaflet.extras' was built under R version 4.5.2# Filtering for 2024 and scale severity as 2024 showed a unusal upward trend in damage
crash_heat_2024 <- crash_nona %>%
filter(year == 2024) %>%
mutate(intensity = case_when(
damage == "Fatal Crash" ~ 15,
damage == "Injury Crash" ~ 3,
damage == "Property Damage Crash" ~ 2, # For intensity of colors on heatmap
TRUE ~ 0
)) %>%
mutate(intensity = intensity * 20) # Increasing intensity scale
# GIS plot
leaflet(crash_heat_2024) %>%
addProviderTiles("CartoDB.DarkMatter") %>% # Selecting theme for map
# Heatmap layer
addHeatmap(
lng = ~longitude, lat = ~latitude,
intensity = ~intensity, # Controls how strong each point affects the heatmap
blur = 15, # Blur effect of heatmap
max = 10, # Adjust max intensity
radius = 20, # Radius of heatmap
gradient = c( "0" = "#440154", "0.2" = "#3b528b",
"0.6" = "#5dc962", "0.8" = "#fde725", "1" = "#ff5a00")
) %>%
# Adding interactive markers with popups
addCircleMarkers(
lng = ~longitude, lat = ~latitude,
radius = 5, # Smaller radius to avoid overcrowding
color = ~ifelse(damage == "Fatal Crash", "purple",
ifelse(damage == "Injury Crash", "orangered", "yellow")),
fillOpacity = 0.5,
stroke = FALSE,
popup = ~paste(
"Year: ", year, "<br>",
"Damage Type: ", damage, "<br>",
"Vehicle Make: ", vehicle_make, "<br>",
"Weather: ", weather
),
group = "markers" # Creating a separate group for the markers so heatmap and points can be viewed separately
) %>%
# Adding layer control to toggle between heatmap and markers
addLayersControl( # Adds a checkbox menu to turn the circle marker layer on/off.
overlayGroups = c("markers"),
options = layersControlOptions(collapsed = FALSE)
)%>%
# For legend in the bottom right
addLegend(
position = "bottomright",
colors = c("purple", "orangered", "yellow"),
labels = c("Fatal Crash", "Injury Crash", "Property Damage Crash"),
title = "Damage Type",
opacity = 1
)# Group and count crashes per speed limit
speed_summary <- crash_nona %>%
group_by(speed_limit, damage) %>%
summarise(crash_count = n()) %>%
arrange(speed_limit, damage)`summarise()` has grouped output by 'speed_limit'. You can override using the
`.groups` argument.ggplot(speed_summary, aes(x = factor(speed_limit), y = crash_count, fill = damage)) +
geom_bar(stat = "identity") +
labs(
title = "Crash Frequency by Speed Limit and Damage Type",
x = "Speed Limit (mph)",
y = "Number of Crashes",
caption = "Source: Data.gov"
) +
theme_minimal() +
scale_fill_manual(values = c("red", "maroon", "pink")) + # Adjust colors as needed
theme(legend.title = element_blank())
# Another GIS plot for speed limit
# Using color palettes to represent speed
pal_speed <- colorNumeric(palette = "viridis", domain = crash_nona$speed_limit) # Defining a numeric color scale (pal_speed) using the "viridis" palette.
leaflet(crash_nona) %>%
addProviderTiles("CartoDB.Positron") %>%
addCircleMarkers(
lng = ~longitude, lat = ~latitude,
color = ~pal_speed(speed_limit),
fillOpacity = 0.6,
radius = 3,
popup = ~paste("Speed Limit: ", speed_limit, "<br>Damage: ", damage)
) %>%
addLegend(
"bottomright", pal = pal_speed, values = ~speed_limit,
title = "Speed Limit (mph)"
) # Adds a legend in the bottom right showing the color gradient for speed limits.# Time series line graph
# Defining selected vehicle makes
selected_makes <- c("TOYOTA", "FORD", "HYUNDAI", "NISSAN", "HONDA",
"MERCEDES", "PORSCHE", "VOLKSWAGEN", "TESLA")
# Filtering and summarizing by year and make
summary_data <- crash_nona %>%
filter(vehicle_make %in% selected_makes) %>%
group_by(year, vehicle_make) %>%
summarise(avg_damage = mean(damage_freq, na.rm = TRUE), .groups = "drop") # Averaging
# Final line plot
p <- ggplot(summary_data, aes(x = factor(year), y = avg_damage, color = vehicle_make, group = vehicle_make)) +
geom_line(linewidth = 1) +
scale_color_viridis_d(option = "magma") +
labs(
title = "Average Damage per Crash by Vehicle Make (2015–2025)",
x = "Year",
y = "Average Damage Frequency",
caption = "Data.gov" # Caption not showing up because of interactivity
) +
theme(panel.background = element_blank(), # Removes background color
plot.background = element_blank(), # Removes plot area background
panel.grid = element_blank() # Removes gridlines
)
# Converts ggplot to interactive plotly plot
ggplotly(p)Understanding the spatial distribution of motor vehicle crashes and the factors associated with them is essential for improving road safety and guiding traffic policies. The data set I chose in this research contains crash records from Montgomery County, Maryland, with variables such as speed limit, crash year, weather condition, damage type, and geographic coordinates. Using R and the leaflet package, I created interactive GIS visualizations to explore how speed and damage type are spatially distributed across the county.
One visualization used a heatmap layered over a base map to show crash concentrations. Locations with a high density of crashes appear as intense areas, mostly clustered around major intersections and highways. Another visualization used circle markers colored by speed limit, allowing for quick assessment of where high-speed zones coincide with crashes. Interactive popups on these markers provided crash-specific information like year, weather, and damage severity.
According to the National Highway Traffic Safety Administration (NHTSA), speeding contributed to 29% of all traffic fatalities in 2021, and crash likelihood increases significantly with speed, especially in adverse weather or uncontrolled intersections (NHTSA, 2022). This supports the patterns seen in the GIS plot: higher-speed zones not only had more crashes but also a higher share of injury-related incidents.