library(knitr)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.4.4 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggplot2)
library(readxl)
wchd20_22 <- read_excel("C:/Users/latimerl/Downloads/wchd20_22.xlsx")
## Warning: Expecting numeric in A9182 / R9182C1: got 'Crash Instance'
## Warning: Expecting numeric in E9182 / R9182C5: got 'Number of Traffic Lanes'
## Warning: Expecting numeric in R9182 / R9182C18: got 'Total Units Reported'
## Warning: Expecting numeric in S9182 / R9182C19: got 'Total Motor Vehicles'
## Warning: Expecting numeric in T9182 / R9182C20: got 'Crash Day'
## Warning: Expecting numeric in AE9182 / R9182C31: got 'Crash Year'
## Warning: Expecting numeric in AQ9182 / R9182C43: got 'Crash Longitude'
## Warning: Expecting numeric in AR9182 / R9182C44: got 'Crash Latitude'
## Warning: Expecting numeric in BJ9182 / R9182C62: got 'Total Non-Motor Vehicles
## (2016+)'
## Warning: Expecting numeric in E22495 / R22495C5: got 'Crash Occurred on Other
## Than a Roadway'
filtered20_22 <- wchd20_22 |> filter(!row_number() %in% c(1981))
filterfinal <- filtered20_22 |> filter(!is.na(`Crash Instance`))
fall <- c('September', 'October','November')
winter <- c('December', 'January', 'February')
spring <- c('March', 'April', 'May')
summer <- c('June', 'July', 'August')
weekday <- c('Monday', 'Tuesday', 'Wednesday', 'Thursday')
weekend <- c('Friday', 'Saturday', 'Sunday')
wchd_new <- filterfinal %>%
mutate(quarter = case_when(`Crash Month`== 'January' | `Crash Month` == 'February' | `Crash Month` == 'March' ~ 'quarter1',
`Crash Month` == 'April' | `Crash Month` == 'May' | `Crash Month`== 'June' ~ 'quarter2',
`Crash Month`== 'July' | `Crash Month` == 'August' | `Crash Month`== 'September' ~ 'quarter3',
`Crash Month`== 'October' | `Crash Month` == 'November' | `Crash Month` == 'December' ~ 'quarter4')) %>%
mutate(weekend = case_when(`Day of Week` == 'Saturday' | `Day of Week` == 'Sunday' ~ 'weekend',
`Day of Week` == 'Monday' | `Day of Week` == 'Tuesday' | `Day of Week` == 'Wednesday' | `Day of Week` == 'Thursday' | `Day of Week` == 'Friday' ~ 'weekday')) %>%
mutate(driversatrisk_age = case_when(`Crash: Young Driver` == 'Driver Age 18-20' | `Crash: Young Driver` == 'Driver Age 21-24' | `Crash: Young Driver` == 'Driver Age 16' | `Crash: Young Driver`== 'Driver Age 17' | `Crash: Young Driver` == 'Driver Age 15' | `Crash: Young Driver`== 'Driver Age 18-21' | `Crash: Young Driver`== 'Driver Age 21-25' | `Crash: Senior Driver`== 'Driver Age 65-74' | `Crash: Senior Driver`== 'Driver Age 75-84' | `Crash: Senior Driver` == 'Driver Age 60-64' | `Crash: Senior Driver` == 'Driver Age 85 or Over' ~ 'atrisk_age',
`Crash: Young Driver`== 'No Driver Age 15-24' | `Crash: Senior Driver` == 'No Driver Under Age 60' ~ 'noatrisk_age')) %>%
mutate(unsafe_conditions = case_when(`Road Conditions` == 'Snow' | `Road Conditions` == 'Ice' | `Road Conditions` == 'Wet' | `Road Conditions` == 'Slush' | `Road Conditions` == 'Mud, Dirt, Gravel' | `Road Conditions` == 'Oily' | `Road Conditions` == 'Water (standing/moving)' | `Road Conditions` == 'Debris' | `Road Conditions` == 'Sand' | `Construction Type` == 'Utility' | `Construction Type` == 'Construction/Maintenance' | `Construction Type` == 'Construction/Maintenance' | `Lighting Conditions` == 'Dawn' | `Lighting Conditions` == 'Dark - Unlighted' | `Lighting Conditions` == 'Dusk' ~ 'unsafe_conditions',
TRUE ~ 'safe_conditions')) %>%
mutate(ped_bike = case_when(`Crash: Pedestrian` == 'Pedestrian Involved' | `Crash: Bicyclist` == 'Bicyclist Involved' ~ 'yes',
TRUE ~ 'no')) %>%
mutate(part_of_day = case_when(`Time of Day` == '5:00 AM - 5:59 AM' | `Time of Day` == '6:00 AM - 6:59 AM' | `Time of Day` == '7:00 AM - 7:59 AM' | `Time of Day` == '8:00 AM - 8:59 AM' | `Time of Day` == '9:00 AM - 9:59 AM' | `Time of Day` == '10:00 AM - 10:59 AM' | `Time of Day` == '11:00 AM - 11:59 AM' ~ 'Morning',
`Time of Day` == '12:00 Noon - 12:59 PM' | `Time of Day` == '1:00 PM - 1:59 PM' | `Time of Day` == '2:00 PM - 2:59 PM' | `Time of Day` == '3:00 PM - 3:59 PM' | `Time of Day` == '4:00 PM - 4:59 PM' ~ 'Afternoon',
`Time of Day` == '5:00 PM - 5:59 PM' | `Time of Day` == '6:00 PM - 6:59 PM' | `Time of Day` == '7:00 PM - 7:59 PM' | `Time of Day` == '8:00 PM - 8:59 PM' ~ 'Evening',
`Time of Day` == '9:00 PM - 9:59 PM' | `Time of Day` == '10:00 PM - 10:59 PM' | `Time of Day` == '11:00 PM - 11:59 PM' | `Time of Day` == '12:00 Midnight - 12:59 AM' | `Time of Day` == '1:00 AM - 1:59 AM' | `Time of Day` == '2:00 AM - 2:59 AM' | `Time of Day` == '3:00 AM - 3:59 AM' | `Time of Day` == '4:00 AM - 4:59 AM' ~ 'Night')) %>%
mutate(season = case_when(`Crash Month` %in% summer ~'Summer',
`Crash Month` %in% fall ~ 'Fall',
`Crash Month` %in% winter ~ 'Winter',
`Crash Month` %in% spring ~ 'Spring')) %>%
mutate(weekends_fri = case_when(`Day of Week` %in% weekday ~ 'Weekday',
`Day of Week` %in% weekend ~ 'Weekend'))
wchd_new2 <- wchd_new %>%
mutate(date_time = paste(`Day of Week`, `Time of Day`)) %>%
mutate(weekends_fri_night = case_when(`Day of Week` == "Saturday" | `Day of Week` == "Sunday" | date_time == "Friday 7:00 PM - 7:59 PM" | date_time == "Friday 5:00 PM - 5:59 PM" | date_time == "Friday 6:00 PM - 6:59 PM" | date_time == "Friday 8:00 PM - 8:59 PM" | date_time == "Friday 9:00 PM - 9:59 PM" | date_time == "Friday 10:00 PM - 10:59 PM" | date_time == "Friday 11:00 PM - 11:59 PM" ~ "Weekend",
TRUE ~ "Weekday"))
wchd_new2_ID <- wchd_new2 %>%
mutate(ID = row_number())
wchd_month <- wchd_new2_ID[c("ID", "Crash Month")]
wchd_year <- wchd_new2_ID[c("ID", "Crash Year")]
wchd_road_conditions <- wchd_new2_ID[c("ID", "Road Conditions")]
wchd_crash_month <- wchd_new2_ID[c("ID", "Crash Month")]
wchd_lighting_conditions <- wchd_new2_ID[c("ID", "Lighting Conditions")]
wchd_crash_day <- wchd_new2_ID[c("ID", "Crash Day")]
wchd_drinking <- wchd_new2_ID[c("ID", "Crash: Drinking")]
wchd_drug_use <- wchd_new2_ID[c("ID", "Crash: Drug Use")]
wchd_pedestrian <- wchd_new2_ID[c("ID", "Crash: Pedestrian")]
wchd_fatal <- wchd_new2_ID[c("ID", "Crash: Fatal Crash")]
wchd_time_of_day <- wchd_new2_ID[c("ID", "Time of Day")]
wchd_crash_longitude <- wchd_new2_ID[c("ID", "Crash Longitude")]
wchd_crash_latitude <- wchd_new2_ID[c("ID", "Crash Latitude")]
wchd_young_driver <- wchd_new2_ID[c("ID", "Crash: Young Driver")]
wchd_senior_driver <- wchd_new2_ID[c("ID", "Crash: Senior Driver")]
wchd_crash_latitude <- wchd_new2_ID[c("ID", "Crash Latitude")]
wchd_day_of_week <- wchd_new2_ID[c("ID", "Day of Week")]
wchd_driversatrisk_age <- wchd_new2_ID[c("ID", "driversatrisk_age")]
wchd_unsafe_conditions <- wchd_new2_ID[c("ID", "unsafe_conditions")]
wchd_ped_bike <- wchd_new2_ID[c("ID", "ped_bike")]
wchd_part_of_day <- wchd_new2_ID[c("ID", "part_of_day")]
wchd_season <- wchd_new2_ID[c("ID", "season")]
wchd_weekends_fri_night <- wchd_new2_ID[c("ID", "weekends_fri_night")]
wchd_weather_conditions <- wchd_new2_ID[c("ID", "Weather Conditions (2016+)")]
wchd_distracted <- wchd_new2_ID[c("ID", "Crash: Driver Distracted (2016+)")]
wchd_merged <- left_join(wchd_month, wchd_year,by = "ID") %>%
left_join(wchd_road_conditions, by = "ID") %>%
left_join(wchd_lighting_conditions, by = "ID") %>%
left_join(wchd_crash_day, by = "ID") %>%
left_join(wchd_drinking, by = "ID") %>%
left_join(wchd_pedestrian, by = "ID") %>%
left_join(wchd_fatal, by = "ID") %>%
left_join(wchd_time_of_day, by = "ID") %>%
left_join(wchd_crash_longitude, by = "ID") %>%
left_join(wchd_crash_latitude, by = "ID") %>%
left_join(wchd_young_driver, by = "ID") %>%
left_join(wchd_senior_driver, by = "ID") %>%
left_join(wchd_day_of_week, by = "ID") %>%
left_join(wchd_unsafe_conditions, by = "ID") %>%
left_join(wchd_driversatrisk_age, by = "ID") %>%
left_join(wchd_weather_conditions, by = "ID") %>%
left_join(wchd_ped_bike, by = "ID") %>%
left_join(wchd_part_of_day, by = "ID") %>%
left_join(wchd_season, by = "ID") %>%
left_join(wchd_weekends_fri_night, by = "ID") %>%
left_join(wchd_distracted, by = "ID")
game_day <- wchd_merged %>%
mutate(date = paste(`Crash Month`, `Crash Day`, `Crash Year`)) %>%
mutate(gameday = case_when(date == "September 4 2021" |date == "September 11 2021" | date == "September 18 2021" | date == "September 25 2021" | date == "October 23 2021" | date == "November 6 2021" | date == "November 27 2021" | date == "September 3 2022" | date == "September 10 2022" | date == "September 17 2022" | date == "September 24 2022" | date == "October 15 2022" | date == "October 29 2022" | date == "November 12 2022" | date == "November 19 2022" ~ "yes",
TRUE ~ "no")) %>%
filter(`Crash Year` != "2020")
wchd_merged <- wchd_merged |>
mutate(fatality_bin = ifelse(`Crash: Fatal Crash` == "Fatal Crash (may also be injuries)", 1, 0),
drinking_bin = ifelse(`Crash: Drinking` == "Drinking Involved", 1, 0),
weekends_bin = ifelse(weekends_fri_night == "Weekend", "1", "0"),
unsafe_bin = ifelse(unsafe_conditions == "unsafe_conditions", "1", "0"),
fall_bin = ifelse(season == "Fall", 1, 0),
night_bin = ifelse(part_of_day == "Night", 1, 0),
ped_bike_bin = ifelse(ped_bike == "yes", 1, 0),
ped_bin = ifelse(`Crash: Pedestrian` == "Pedestrian Involved", 1, 0),
distracted_bin = ifelse(`Crash: Driver Distracted (2016+)` == "Driver Distracted", 1, 0),
agerisk_bin = ifelse(driversatrisk_age == "atrisk_age", 1, 0),
morning_night_bin = ifelse(part_of_day == "Morning"|part_of_day == "Night", 1, 0))
wchd_merged |>
count(fatality_bin, `Crash: Fatal Crash`)
## # A tibble: 2 × 3
## fatality_bin `Crash: Fatal Crash` n
## <dbl> <chr> <int>
## 1 0 No Fatalities (may be injuries) 25014
## 2 1 Fatal Crash (may also be injuries) 72
wchd_merged |>
count(drinking_bin, `Crash: Drinking`)
## # A tibble: 2 × 3
## drinking_bin `Crash: Drinking` n
## <dbl> <chr> <int>
## 1 0 No Drinking Involved 24230
## 2 1 Drinking Involved 856
wchd_merged |>
count(weekends_bin, weekends_fri_night)
## # A tibble: 2 × 3
## weekends_bin weekends_fri_night n
## <chr> <chr> <int>
## 1 0 Weekday 17843
## 2 1 Weekend 7243
wchd_merged |>
count(unsafe_bin, unsafe_conditions)
## # A tibble: 2 × 3
## unsafe_bin unsafe_conditions n
## <chr> <chr> <int>
## 1 0 safe_conditions 14774
## 2 1 unsafe_conditions 10312
wchd_merged |>
count(fall_bin)
## # A tibble: 2 × 2
## fall_bin n
## <dbl> <int>
## 1 0 17555
## 2 1 7531
wchd_merged |>
count(ped_bike_bin)
## # A tibble: 2 × 2
## ped_bike_bin n
## <dbl> <int>
## 1 0 24665
## 2 1 421
wchd_merged |>
count(`Crash: Driver Distracted (2016+)`, distracted_bin)
## # A tibble: 2 × 3
## `Crash: Driver Distracted (2016+)` distracted_bin n
## <chr> <dbl> <int>
## 1 Driver Distracted 1 1653
## 2 No Driver Distraction 0 23433
wchd_merged |>
count(agerisk_bin, driversatrisk_age)
## # A tibble: 2 × 3
## agerisk_bin driversatrisk_age n
## <dbl> <chr> <int>
## 1 0 noatrisk_age 11670
## 2 1 atrisk_age 13416
game_day <- game_day %>%
mutate(gameday_bin = ifelse(gameday == "yes", 1, 0),
drinking_bin = ifelse(`Crash: Drinking` == "Drinking Involved", 1, 0),
fatal_bin = ifelse(`Crash: Fatal Crash` == "Fatal Crash (may also be injuries)", 1, 0),
distracted_bin = ifelse(`Crash: Driver Distracted (2016+)` == "Driver Distracted", 1, 0),
ped_bike_bin = ifelse(ped_bike == "yes", 1, 0))
log_reg <- function(x) {
{{x}} |>
broom::tidy() |>
mutate(or = exp(estimate))
}
model1 <- glm(fatality_bin ~ weekends_bin, data = wchd_merged, family = "binomial")
log_reg(model1)
## # A tibble: 2 × 6
## term estimate std.error statistic p.value or
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -6.15 0.162 -37.9 0 0.00213
## 2 weekends_bin1 0.793 0.236 3.35 0.000799 2.21
#The odds of a fatal crash are 2.2 times higher during the weekends than during weekdays.
model2 <- glm(fatality_bin ~ weekends_bin + drinking_bin, data = wchd_merged, family = "binomial")
log_reg(model2)
## # A tibble: 3 × 6
## term estimate std.error statistic p.value or
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -6.37 0.172 -37.0 3.05e-300 0.00171
## 2 weekends_bin1 0.504 0.244 2.07 3.87e- 2 1.66
## 3 drinking_bin 2.49 0.262 9.47 2.68e- 21 12.0
#The odds of a fatal crash are 1.65 times higher during the weekends than during weekdays when controlling for drinking.
model3 <- glm(fatality_bin ~ weekends_bin + drinking_bin + weekends_bin*drinking_bin, data = wchd_merged, family = "binomial")
log_reg(model3)
## # A tibble: 4 × 6
## term estimate std.error statistic p.value or
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -6.36 0.183 -34.8 1.07e-265 0.00172
## 2 weekends_bin1 0.487 0.294 1.66 9.72e- 2 1.63
## 3 drinking_bin 2.45 0.401 6.12 9.50e- 10 11.6
## 4 weekends_bin1:drinking_bin 0.0577 0.532 0.109 9.14e- 1 1.06
#Since the OR of the interaction between weekends and drinking is ~1, there is not a significant interaction between the two variables
model4 <- glm(fatality_bin ~ unsafe_bin, data = wchd_merged, family = "binomial")
log_reg(model4)
## # A tibble: 2 × 6
## term estimate std.error statistic p.value or
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -6.34 0.196 -32.3 6.29e-229 0.00176
## 2 unsafe_bin1 0.933 0.246 3.80 1.47e- 4 2.54
#The odds of a fatal crash are 2.54 times higher in unsafe conditions than safe conditions.
model5 <- glm(fatality_bin ~ drinking_bin + unsafe_bin + weekends_bin, data = wchd_merged, family = "binomial")
log_reg(model5)
## # A tibble: 4 × 6
## term estimate std.error statistic p.value or
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -6.74 0.224 -30.1 2.89e-199 0.00118
## 2 drinking_bin 2.41 0.263 9.16 5.05e- 20 11.2
## 3 unsafe_bin1 0.771 0.248 3.10 1.92e- 3 2.16
## 4 weekends_bin1 0.444 0.245 1.82 6.94e- 2 1.56
#The odds of a fatal crash are 11.16 times higher when drinking is involved compared to when drinking is not involved when controlling for conditions and the weekend vs weekday.
model6 <- glm(drinking_bin ~ night_bin + weekends_bin, data = wchd_merged, family = "binomial")
log_reg(model6)
## # A tibble: 3 × 6
## term estimate std.error statistic p.value or
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -4.26 0.0601 -70.9 0 0.0141
## 2 night_bin 1.96 0.0728 27.0 2.55e-160 7.13
## 3 weekends_bin1 0.721 0.0727 9.91 3.61e- 23 2.06
#The odds of a crash involving drinking are 7.1 times higher at night compared to other times of day when controlling for weekend vs weekday.
model7 <- glm(fatality_bin ~ ped_bike_bin, data = wchd_merged, family = "binomial")
log_reg(model7)
## # A tibble: 2 × 6
## term estimate std.error statistic p.value or
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -6.05 0.131 -46.0 0 0.00236
## 2 ped_bike_bin 2.68 0.302 8.88 6.81e-19 14.6
#The odds of a crash being fatal are 14.59 times higher if a pedestrian or biker is involved compared to no pedestrian or biker involved. --> is this caused by one more than the other?
model8 <- glm(fatality_bin ~ ped_bin + unsafe_bin + drinking_bin, data = wchd_merged, family = "binomial")
log_reg(model8)
## # A tibble: 4 × 6
## term estimate std.error statistic p.value or
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -6.72 0.210 -31.9 1.23e-223 0.00121
## 2 ped_bin 2.81 0.367 7.65 1.98e- 14 16.6
## 3 unsafe_bin1 0.795 0.250 3.18 1.48e- 3 2.21
## 4 drinking_bin 2.44 0.263 9.31 1.29e- 20 11.5
#The odds of a crash being fatal are 16.58 times higher if a pedestrian is involved compared to no pedestrian being involved when controlling for driving conditions and drinking.
model9 <- glm(drinking_bin ~ gameday_bin, data = game_day, family = "binomial")
log_reg(model9)
## # A tibble: 2 × 6
## term estimate std.error statistic p.value or
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -3.40 0.0430 -79.2 0 0.0333
## 2 gameday_bin 0.375 0.205 1.83 0.0677 1.46
#The odds of a crash involving drinking are 1.45 times higher on a game day
model10 <- glm(fatal_bin ~ gameday_bin, data = game_day, family = "binomial")
log_reg(model10)
## # A tibble: 2 × 6
## term estimate std.error statistic p.value or
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -5.95 0.149 -39.9 0 0.00260
## 2 gameday_bin 1.01 0.524 1.94 0.0526 2.76
#The odds of a fatal crash are 2.76 times higher on a game day
model11 <- glm(distracted_bin ~ gameday_bin, data = game_day, family = "binomial")
log_reg(model11)
## # A tibble: 2 × 6
## term estimate std.error statistic p.value or
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -2.64 0.0305 -86.7 0 0.0712
## 2 gameday_bin -0.233 0.190 -1.23 0.220 0.792
#The odds of a crash with a distracted driver are .79 times less on a game day
model12 <- glm(ped_bike_bin ~ gameday_bin, data = game_day, family = "binomial")
log_reg(model12)
## # A tibble: 2 × 6
## term estimate std.error statistic p.value or
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -4.10 0.0598 -68.5 0 0.0166
## 2 gameday_bin 0.183 0.310 0.588 0.556 1.20
#The odds of a crash involving a pedestrian or biker are 1.2 times higher on a game day