Dataset Profile

The dataset contains 48 variables.

## # A tibble: 5 × 48
##   CRASH_RECORD_ID                 CRASH_DATE_EST_I CRASH_DATE POSTED_SPEED_LIMIT
##   <chr>                           <chr>            <chr>                   <dbl>
## 1 23a79931ef555d54118f64dc9be2cf… <NA>             09/05/202…                 30
## 2 2675c13fd0f474d730a5b780968b3c… <NA>             09/22/202…                 50
## 3 5f54a59fcb087b12ae5b1acff96a3c… <NA>             07/29/202…                 30
## 4 7ebf015016f83d09b321afd671a836… <NA>             08/09/202…                 30
## 5 6c1659069e9c6285a650e70d6f9b57… <NA>             08/18/202…                 15
## # ℹ 44 more variables: TRAFFIC_CONTROL_DEVICE <chr>, DEVICE_CONDITION <chr>,
## #   WEATHER_CONDITION <chr>, LIGHTING_CONDITION <chr>, FIRST_CRASH_TYPE <chr>,
## #   TRAFFICWAY_TYPE <chr>, LANE_CNT <dbl>, ALIGNMENT <chr>,
## #   ROADWAY_SURFACE_COND <chr>, ROAD_DEFECT <chr>, REPORT_TYPE <chr>,
## #   CRASH_TYPE <chr>, INTERSECTION_RELATED_I <chr>, NOT_RIGHT_OF_WAY_I <chr>,
## #   HIT_AND_RUN_I <chr>, DAMAGE <chr>, DATE_POLICE_NOTIFIED <chr>,
## #   PRIM_CONTRIBUTORY_CAUSE <chr>, SEC_CONTRIBUTORY_CAUSE <chr>, …

Before we embark on any predictive analysis, it is good practice to get rid of non-informative features that are to have little or no effect to our modelling such as CRASH_RECORD_ID which is an identifier along with other sparse features in the dataset.

data <- data %>% select(-c(CRASH_RECORD_ID,STREET_NAME,LOCATION,PHOTOS_TAKEN_I,STATEMENTS_TAKEN_I,
                   DOORING_I,WORK_ZONE_I,WORK_ZONE_TYPE,WORK_ZONE_TYPE,WORKERS_PRESENT_I,
                   INTERSECTION_RELATED_I,NOT_RIGHT_OF_WAY_I,HIT_AND_RUN_I,CRASH_DATE_EST_I,LANE_CNT,
                   STREET_NAME,STREET_NO,LATITUDE, LONGITUDE,DATE_POLICE_NOTIFIED,
                   BEAT_OF_OCCURRENCE,NUM_UNITS))
data %>% head

## # A tibble: 6 × 28
##   CRASH_DATE          POSTED_SPEED_LIMIT TRAFFIC_CONTROL_DEVICE DEVICE_CONDITION
##   <chr>                            <dbl> <chr>                  <chr>           
## 1 09/05/2023 07:05:0…                 30 TRAFFIC SIGNAL         FUNCTIONING PRO…
## 2 09/22/2023 06:45:0…                 50 NO CONTROLS            NO CONTROLS     
## 3 07/29/2023 02:45:0…                 30 TRAFFIC SIGNAL         FUNCTIONING PRO…
## 4 08/09/2023 11:00:0…                 30 NO CONTROLS            NO CONTROLS     
## 5 08/18/2023 12:50:0…                 15 OTHER                  FUNCTIONING PRO…
## 6 11/26/2019 08:38:0…                 25 NO CONTROLS            NO CONTROLS     
## # ℹ 24 more variables: WEATHER_CONDITION <chr>, LIGHTING_CONDITION <chr>,
## #   FIRST_CRASH_TYPE <chr>, TRAFFICWAY_TYPE <chr>, ALIGNMENT <chr>,
## #   ROADWAY_SURFACE_COND <chr>, ROAD_DEFECT <chr>, REPORT_TYPE <chr>,
## #   CRASH_TYPE <chr>, DAMAGE <chr>, PRIM_CONTRIBUTORY_CAUSE <chr>,
## #   SEC_CONTRIBUTORY_CAUSE <chr>, STREET_DIRECTION <chr>,
## #   MOST_SEVERE_INJURY <chr>, INJURIES_TOTAL <dbl>, INJURIES_FATAL <dbl>,
## #   INJURIES_INCAPACITATING <dbl>, INJURIES_NON_INCAPACITATING <dbl>, …

Feature Engineering

This now leaves us with 28 variables. The next step we need to take is to ensure we format the variable into their appropriate data types. For instance, variable CRASH_MONTH is formatted as a double, and feeding this dataset in its current state would treat such a variable as a continuous numeric double rather than a discrete value. While at it, I will add the number of minutes to the variable CRASH_HOUR so it’s formatted as a continuous value from 0-24 expressing the number of hours from midnight an accident occurred. We will also add CRASH_YEARfrom the variable CRASH_DATE.

data<-data %>% mutate(CRASH_DATE=mdy_hms(CRASH_DATE),
                      CRASH_YEAR = as_factor(year(CRASH_DATE)),
                CRASH_HOUR=round(hour(CRASH_DATE)+minute(CRASH_DATE)/60,2),
                CRASH_YEAR=as_factor(year(CRASH_DATE))) %>% select(CRASH_DATE,everything())

data<-data %>% mutate(across(where(is.character), as.factor),
                INJURIES_TOTAL=as.integer(INJURIES_TOTAL),
                INJURIES_FATAL=as.integer(INJURIES_FATAL),
                INJURIES_INCAPACITATING=as.integer(INJURIES_INCAPACITATING),
                INJURIES_NON_INCAPACITATING=as.integer(INJURIES_NON_INCAPACITATING),
                INJURIES_REPORTED_NOT_EVIDENT=as.integer(INJURIES_REPORTED_NOT_EVIDENT),
                INJURIES_NO_INDICATION=as.integer(INJURIES_NO_INDICATION),
                INJURIES_UNKNOWN=as.integer(INJURIES_UNKNOWN),
                CRASH_DAY_OF_WEEK=as.factor(CRASH_DAY_OF_WEEK),
                CRASH_MONTH=as.factor(CRASH_MONTH))

Then we will remove CRASH_DATE since it’s now redundant. We could also use one of imputation methods like regression imputation to fill missing values. However, we will omit missing values for the sake of simplicity.

data<-data %>% 
  select(-CRASH_DATE) %>% 
  select(CRASH_DAY_OF_WEEK,CRASH_MONTH,CRASH_YEAR,CRASH_HOUR,everything())
data <- na.omit(data)
data %>% head

## # A tibble: 6 × 28
##   CRASH_DAY_OF_WEEK CRASH_MONTH CRASH_YEAR CRASH_HOUR POSTED_SPEED_LIMIT
##   <fct>             <fct>       <fct>           <dbl>              <dbl>
## 1 3                 9           2023            19.1                  30
## 2 6                 9           2023            18.8                  50
## 3 7                 7           2023            14.8                  30
## 4 4                 8           2023            23                    30
## 5 6                 8           2023            12.8                  15
## 6 3                 11          2019             8.63                 25
## # ℹ 23 more variables: TRAFFIC_CONTROL_DEVICE <fct>, DEVICE_CONDITION <fct>,
## #   WEATHER_CONDITION <fct>, LIGHTING_CONDITION <fct>, FIRST_CRASH_TYPE <fct>,
## #   TRAFFICWAY_TYPE <fct>, ALIGNMENT <fct>, ROADWAY_SURFACE_COND <fct>,
## #   ROAD_DEFECT <fct>, REPORT_TYPE <fct>, CRASH_TYPE <fct>, DAMAGE <fct>,
## #   PRIM_CONTRIBUTORY_CAUSE <fct>, SEC_CONTRIBUTORY_CAUSE <fct>,
## #   STREET_DIRECTION <fct>, MOST_SEVERE_INJURY <fct>, INJURIES_TOTAL <int>,
## #   INJURIES_FATAL <int>, INJURIES_INCAPACITATING <int>, …

Our dataset now has the following metadata: CRASH_DAY_OF_WEEK: day of crash.

CRASH_MONTH: month of crash

CRASH_YEAR: Year of crash

CRASH_HOUR: Number of hours since midnight when crash occurred

POSTED_SPEED_LIMIT: Speed limit of where crash occurred

TRAFFIC_CONTROL_DEVICE - The control device at the crash site

DEVICE_CONDITION - The condition of the device

WEATHER_CONDITION - weather condition

LIGHTING_CONDITION - Lighting condition

FIRST_CRASH_TYPE - first crash type

TRAFFICWAY_TYPE - the traffic type at the scene

ALIGNMENT - alignment of the vehicle

ROADWAY_SURFACE_COND - road condition

ROAD EFFECT - road defect

REPORT_TYPE - whether report type on scene

CRASH_TYPE - the crash type

DAMAGE - damage value

PRIM_CONTRIBUTORY_CAUSE - primary contributory cause

SEC_CONTRIBUTORY_CAUSE - secondary contributory cause

STREET_DIRECTION - street direction

MOST_SEVERE_INJURY - most severe injury

INJURIES_TOTAL - total number of injuries

INJURIES_FATAL - fatal injuries

INJURIES_INCAPACITATING - incapacitating injuries

INJURIES_NON_INCAPACITATING - non incapacitating injuries

INJURIES_REPORTED_NOT_EVIDENT - accidents without injuries reported

INJURIES_NO_INDICATION - no indication of injuries

INJURIES_UNKNOWN - unknown number of injuries

Explanatory Data Analysis

Since 2013, most crashes in Chicago have occurred in the 11-17 hrs window. This is expected because traffic is busy at this time when most people are getting off work. Conversely, hrs 0-7 and 20-24 see lesser crashes since most people are at home and there is less traffic.

data %>% 
  ggplot(aes(x=CRASH_HOUR)) +
  geom_histogram(fill = "skyblue", color = "black", bins = 20) +
  labs(title = "Chicago Traffic Crash Hours Distribution",
       x = "Crash Hour",
       y = "Frequency",
       caption="Source: data.gov") +
  theme_minimal()

Most causes of crashes could not be determined. The main primary contributory cause was failing to yield right of way.

data %>% group_by(PRIM_CONTRIBUTORY_CAUSE) %>%
  summarize(count=n(),fatalities=sum(INJURIES_FATAL)) %>%
  arrange(desc(count)) %>% head(5)

## # A tibble: 5 × 3
##   PRIM_CONTRIBUTORY_CAUSE        count fatalities
##   <fct>                          <int>      <int>
## 1 UNABLE TO DETERMINE           306564        257
## 2 FAILING TO YIELD RIGHT-OF-WAY  86219         60
## 3 FOLLOWING TOO CLOSELY          77288          4
## 4 NOT APPLICABLE                 41595         69
## 5 IMPROPER OVERTAKING/PASSING    38807          9

Physical condition of the driver is also a leading cause of crash fatalities.

data %>% group_by(PRIM_CONTRIBUTORY_CAUSE) %>%
  summarize(count=n(),fatalities=sum(INJURIES_FATAL)) %>%
  arrange(desc(fatalities)) %>% head(5)

## # A tibble: 5 × 3
##   PRIM_CONTRIBUTORY_CAUSE                 count fatalities
##   <fct>                                   <int>      <int>
## 1 UNABLE TO DETERMINE                    306564        257
## 2 PHYSICAL CONDITION OF DRIVER             4682         76
## 3 NOT APPLICABLE                          41595         69
## 4 FAILING TO REDUCE SPEED TO AVOID CRASH  33341         63
## 5 FAILING TO YIELD RIGHT-OF-WAY           86219         60

vals<-data %>% group_by(MOST_SEVERE_INJURY) %>% 
  summarize(proportion=n()/nrow(.)) 

plot_ly(labels = vals$MOST_SEVERE_INJURY, values = vals$proportion) %>% 
  add_pie(hole = 0.4)%>% 
  layout(title = "Most Severe Injury by Proportion:Hover for More Details",  showlegend = F,
                      xaxis = list(showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE),
                      yaxis = list(showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE))

data %>% group_by(DAMAGE) %>% 
  summarize(count=n()) %>% 
  ggplot(aes(x=DAMAGE,y=count)) +
  geom_col(width = 0.8, fill="grey")+
  scale_y_continuous(labels = scales::number_format()) +
  labs(title = "Traffic Crash Damage Frequency Bar",
       x = "Damage Value",
       y = "Frequency",
       caption="Source: data.gov") +
  theme_minimal() +                       
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Train-test Splitting

One of the main is splitting a dataset into a train and test dataset. The train is used to fit the model and test dataset is used to evaluate the performance of the model. We make use of the caret library to achieve this purpose. To avoid lengthy training times when dealing with huge amounts of data, it is common to sample a proportion of the dataset to build a model before training the model on the entire dataset. We will sample 10% of the data for our modelling.

set.seed(123)
sample_data<-sample_frac(data,size=.1,replace=F)
train_indices <- createDataPartition(sample_data$CRASH_TYPE, p = 0.8, list = FALSE)
train_data <- sample_data[train_indices, ]
test_data <- sample_data[-train_indices, ]

Then, we would need to fit the model and evaluate. Common metrics for regression include mean absolute error (MAE) and Root Mean Squared Error(RMSE).

lm_model <- linear_reg() %>%
  set_engine("lm") %>%
  set_mode("regression")

#fit the mlr model
lm_spec <- linear_reg() %>%
  set_mode("regression") %>%
  set_engine("lm")

lm_spec

## Linear Regression Model Specification (regression)
## 
## Computational engine: lm

mlr_mod <- lm_spec %>% 
fit(CRASH_HOUR ~ ., data = train_data)

test_aug <- augment(mlr_mod, new_data = test_data)
paste("mae",mean(abs(test_aug$.resid)))

## [1] "mae 4.20686181894546"

paste("rmse:",sqrt(mean((test_aug$.resid)**2)))

## [1] "rmse: 5.35141005390315"

The model has an MAE of 4.22 meaning that its prediction of the time the crash occurred is off by 4.22 hours on average. So, how did the algorithm estimate the time a crash occurred. To understand how the modelling process we will break down our illustrations into five components/objectives.

Objective 1: Describe probability as a foundation of statistical modeling, including inference and maximum likelihood estimation

When trying to predict a value, we try to predict what value a certain sample takes given a set of features. Understanding the distribution of the underlying distribution will help model the likelihood function which maximizes the probability that the predicted value is as close as possible to the observed value. We achieve this by deriving the likelihood function from the probability distribution assumed from the data then estimating the parameters that are most suitable for the given values. In our traffic crash data regression problem we are maximizing the coefficients independent variables which are the parameters in this case in order to give the best guess of when a crash occurred.

\[\hat{CRASH\_HOUR} = \beta_1CRASH\_MONTH12 +\beta_1CRASH\_MONTH11+...+Intercept\_term\]

We want to maximize the likelihood function by finding the best parameters beta_n such that the CRASH_HOUR_hat is close as possible to crash_hour. Finding the best parameters is known as maximum likelihood estimation. Note that for qualitative variables like CRASH_MONTH, we find a parameter for each level of the dummy variable.

Furthermore, we could also model predicting a binary variable. Let’s assume we need to know the variable ‘CRASH_TYPE’. It is a binary variable with values INJURY AND / OR TOW DUE TO CRASH and NO INJURY / DRIVE AWAY. To model the probability of this model we need to know the unconditional probability of a sample k from the population.

paste("Unconditional probability:",mean(data$CRASH_TYPE == "NO INJURY / DRIVE AWAY"))

## [1] "Unconditional probability: 0.737670957437766"

The unconditional probability of a sample k being classified as NO INJURY / DRIVE AWAY is

\[P(k)=0.738\]

In other Words we can guess the value of CRASH_TYPE to be NO INJURY / DRIVE AWAY and we would be correct 73.8% of the time. This probability is known as the naive model.

Given now a vector X of features such as CRASH_MONTH, PRIM_CONTRIBUTORY_CAUSE,DAMAGE we can reduce the Total Probability of Misclassification by fitting a logistic regression which maximizes the following likelihood function by finding the best parameters \[\begin{equation*} \log\left(\frac{p_i}{1-p_i}\right) = \beta_0 + \beta_1 X \end{equation*}\]

which is the log-odds of “success” derived from the odds of success and defined by:

\[ \frac{Pr(Y_i = NO\ INJURY / DRIVE AWAY)} {Pr(Y_i = INJURY\ AND\ /\ OR\ TOW\ DUE\ TO\ CRASH)} = \frac{p_i}{1-p_i}\]

We can further simplify this by encoding our levels of the response variable to dummy variable such that the probability of “sucesss” is 1 and “failure” is 0. \[ \frac{Pr(Y_i = 1)} {Pr(Y_i = 0)} = \frac{p_i}{1-p_i}\] Ultimately, by taking the inverse of the natural log which is the exponent we end up with this likelihood function. \[p(x) = \frac{1}{1+e^{-(\beta_0+\beta_1x_1+\beta_2x_2+....)}}\]
where \(x_1,x_2,....=X\) are the corresponding independent variables. This log-likelihood function is a sigmoid function that gets the probability \(p(x)\) as close to 1 as possible if it’s a “success” and \(1-p(x)\) to 0 if it’s a “failure”. This will now yield the conditional/posterior probability: \[Pr(k|X) = q\] where q is hopefully a number higher than that given by the naive model.

Objective 2: Determine and apply the appropriate generalized linear model for a specific data context

# Specify the logistic regression model
log_reg_spec <- logistic_reg() %>%
  set_engine("glm") %>%
  set_mode("classification")

# Fit the logistic regression model
log_reg_fit <- log_reg_spec %>%
  fit(CRASH_TYPE ~ ., data = train_data)

# Print the model summary
summary(log_reg_fit)

##              Length Class        Mode     
## lvl           2     -none-       character
## spec          8     logistic_reg list     
## fit          30     glm          list     
## preproc       1     -none-       list     
## elapsed       1     -none-       list     
## censor_probs  0     -none-       list

# Make predictions on the testing set
predictions <- predict(log_reg_fit, test_data) %>%
  bind_cols(test_data) %>%
  select(.pred_class, CRASH_TYPE)

paste("Logistic regression on predicting crashtype:",
      mean(predictions$.pred_class==predictions$CRASH_TYPE))

## [1] "Logistic regression on predicting crashtype: 0.899841420869014"

Objective 3: Conduct model selection for a set of candidate models

train_data <- train_data %>% select(-c(MOST_SEVERE_INJURY:INJURIES_UNKNOWN))
test_data<- test_data %>% select(-c(MOST_SEVERE_INJURY:INJURIES_UNKNOWN))
train_data %>% head(5)

## # A tibble: 5 × 20
##   CRASH_DAY_OF_WEEK CRASH_MONTH CRASH_YEAR CRASH_HOUR POSTED_SPEED_LIMIT
##   <fct>             <fct>       <fct>           <dbl>              <dbl>
## 1 7                 4           2020            17                    30
## 2 7                 4           2018            11.5                  30
## 3 3                 8           2018            16.7                  30
## 4 5                 4           2019            18.2                  30
## 5 1                 2           2023             7.77                 10
## # ℹ 15 more variables: TRAFFIC_CONTROL_DEVICE <fct>, DEVICE_CONDITION <fct>,
## #   WEATHER_CONDITION <fct>, LIGHTING_CONDITION <fct>, FIRST_CRASH_TYPE <fct>,
## #   TRAFFICWAY_TYPE <fct>, ALIGNMENT <fct>, ROADWAY_SURFACE_COND <fct>,
## #   ROAD_DEFECT <fct>, REPORT_TYPE <fct>, CRASH_TYPE <fct>, DAMAGE <fct>,
## #   PRIM_CONTRIBUTORY_CAUSE <fct>, SEC_CONTRIBUTORY_CAUSE <fct>,
## #   STREET_DIRECTION <fct>

# Specify the logistic regression model
log_reg_spec <- logistic_reg() %>%
  set_engine("glm") %>%
  set_mode("classification")

# Fit the logistic regression model
log_reg_fit <- log_reg_spec %>%
  fit(CRASH_TYPE ~ ., data = train_data)

# Print the model summary
summary(log_reg_fit)

##              Length Class        Mode     
## lvl           2     -none-       character
## spec          8     logistic_reg list     
## fit          30     glm          list     
## preproc       1     -none-       list     
## elapsed       1     -none-       list     
## censor_probs  0     -none-       list

# Make predictions on the testing set
predictions <- predict(log_reg_fit, test_data) %>%
  bind_cols(test_data) %>%
  select(.pred_class, CRASH_TYPE)

paste("Logistic regression on predicting crashtype of reduced model:",
      sum(predictions$.pred_class==test_data$CRASH_TYPE)/length(test_data$CRASH_TYPE))

## [1] "Logistic regression on predicting crashtype of reduced model: 0.845543926419283"

nn_model <- nnet(CRASH_TYPE ~ ., 
                 data = train_data, 
                  size = 2,  # number of units in hidden layer
                 linout = FALSE)  # set linear output to false since this is a classification problem

## # weights:  423
## initial  value 57666.590131 
## iter  10 value 32241.683740
## iter  20 value 25502.603773
## iter  30 value 24456.627129
## iter  40 value 24153.702556
## iter  50 value 24030.958886
## iter  60 value 23940.880556
## iter  70 value 23871.464482
## iter  80 value 23819.204058
## iter  90 value 23791.132447
## iter 100 value 23782.698298
## final  value 23782.698298 
## stopped after 100 iterations

# Make predictions on the test set
predicted_classes <- predict(nn_model, select(test_data,-CRASH_TYPE), type = "class")

paste("Accuracy with neural network model:",mean(predicted_classes==test_data$CRASH_TYPE))

## [1] "Accuracy with neural network model: 0.839137329527434"

log_reg_spec <- logistic_reg() %>%
  set_engine("glm") %>%
  set_mode("classification")

data_recipe <- recipe(CRASH_TYPE ~ ., data = train_data)

log_reg_wf <- workflow() %>%
  add_model(log_reg_spec) %>%
  add_recipe(data_recipe)

log_reg_res <- log_reg_wf %>%
  fit_resamples(resamples = vfold_cv(train_data, v = 5))

## → A | error:   factor CRASH_YEAR has new levels 2014

## 
There were issues with some computations   A: x1

There were issues with some computations   A: x1

# Get performance metrics
log_reg_metrics <- collect_metrics(log_reg_res)
paste("Cross validated accuracy",log_reg_metrics$mean[1])

## [1] "Cross validated accuracy 0.842557045802929"

Objective 4: Communicate the results of statistical models to a general audience

Communicating the results of your analysis to a non-technical audience is one of the hallmarks of effective statistical modelling. This will ensure that stakeholders involved can turn insights about the data into actionable items.

When communicating results to a general audience it is important to include visualizations and descriptions that would help interpret technical terms in an easy to understand manner.

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot. One such plot is the ROC curve which shows how our model is performing with respect to True Positive Rate otherwise known as the recall against the False Positive Rate i.e type I error rate.

\[TPR=Sensitivity = \frac{True\ Positives}{True\ Positives + False\ Negatives}\]

\[FPR = 1-Specificity = \frac{False\ Positives}{False\ Positives + True\ Negatives}\]

roc_curve <- roc(as.integer(test_data$CRASH_TYPE), as.integer(predictions$.pred_class))

plot(roc_curve, main = "ROC Curve for Different Threshold Values", col = "blue")

legend("bottomright", legend = paste("AUC =",
  round(auc(roc_curve), 2)), col = "blue", lwd = 2)

The solid line represents the hypothetical worst model at each threshold value. We aim for a model that extends as far as possible to the top-left corner.The top left is the best possible model giving 100% sensitivity and 0% FPR giving a unit area under the curve. It’s important to choose a suitable threshold value. We sometimes prefer to reduce one kind of error over the other. For example, we do not want a crash that resulted with injuries to be labelled as no injuries, so we try to minimize false negatives.

table(predictions$.pred_class, test_data$CRASH_TYPE)

##                                   
##                                    INJURY AND / OR TOW DUE TO CRASH
##   INJURY AND / OR TOW DUE TO CRASH                             2789
##   NO INJURY / DRIVE AWAY                                       1374
##                                   
##                                    NO INJURY / DRIVE AWAY
##   INJURY AND / OR TOW DUE TO CRASH                   1061
##   NO INJURY / DRIVE AWAY                            10541

Objective 5: Use programming software (i.e., R) to fit and assess statistical models

After trying various models, we can now train our choice model on the full data set rather than a sample using a programming software such as R or Python. Here is the full end-to-end procedure of implementing our logistic regression problem in R.

train_indices <- createDataPartition(data$CRASH_TYPE, p = 0.7, list = FALSE)
#remove variables that cause information leakage
train_data <- data[train_indices, ] %>% 
  select(-c(MOST_SEVERE_INJURY:INJURIES_UNKNOWN))
test_data <- data[-train_indices, ] %>% 
  select(-c(MOST_SEVERE_INJURY:INJURIES_UNKNOWN))
train_data %>% head

## # A tibble: 6 × 20
##   CRASH_DAY_OF_WEEK CRASH_MONTH CRASH_YEAR CRASH_HOUR POSTED_SPEED_LIMIT
##   <fct>             <fct>       <fct>           <dbl>              <dbl>
## 1 6                 9           2023            18.8                  50
## 2 7                 7           2023            14.8                  30
## 3 4                 8           2023            23                    30
## 4 6                 8           2023            12.8                  15
## 5 3                 11          2019             8.63                 25
## 6 2                 2           2023            17.5                  30
## # ℹ 15 more variables: TRAFFIC_CONTROL_DEVICE <fct>, DEVICE_CONDITION <fct>,
## #   WEATHER_CONDITION <fct>, LIGHTING_CONDITION <fct>, FIRST_CRASH_TYPE <fct>,
## #   TRAFFICWAY_TYPE <fct>, ALIGNMENT <fct>, ROADWAY_SURFACE_COND <fct>,
## #   ROAD_DEFECT <fct>, REPORT_TYPE <fct>, CRASH_TYPE <fct>, DAMAGE <fct>,
## #   PRIM_CONTRIBUTORY_CAUSE <fct>, SEC_CONTRIBUTORY_CAUSE <fct>,
## #   STREET_DIRECTION <fct>

# Specify the logistic regression model
log_reg_spec <- logistic_reg() %>%
  set_engine("glm") %>%
  set_mode("classification")

# Fit the logistic regression model
log_reg_fit <- log_reg_spec %>%
  fit(CRASH_TYPE ~ ., data = train_data)

# Print the model summary
summary(log_reg_fit)

##              Length Class        Mode     
## lvl           2     -none-       character
## spec          8     logistic_reg list     
## fit          30     glm          list     
## preproc       1     -none-       list     
## elapsed       1     -none-       list     
## censor_probs  0     -none-       list

# Make predictions on the testing set
predictions <- predict(log_reg_fit, test_data) %>%
  bind_cols(test_data) %>%
  select(.pred_class, CRASH_TYPE)

calculate_metrics <- function(true_labels, predicted_labels) {
  TP <- sum(true_labels == "NO INJURY / DRIVE AWAY" & 
              predicted_labels == "NO INJURY / DRIVE AWAY")
  FP <- sum(true_labels == "INJURY AND / OR TOW DUE TO CRASH" &
              predicted_labels == "NO INJURY / DRIVE AWAY")
  FN <- sum(true_labels == "NO INJURY / DRIVE AWAY" &
              predicted_labels == "INJURY AND / OR TOW DUE TO CRASH")
  TN <- sum(true_labels == "INJURY AND / OR TOW DUE TO CRASH" & 
              predicted_labels == "INJURY AND / OR TOW DUE TO CRASH")
  
  precision <- TP / (TP + FP)
  recall <- TP / (TP + FN)
  accuracy <- (TP + TN) / length(true_labels)
  
  return(list(precision = precision, recall = recall, accuracy = accuracy))
}

metrics <- calculate_metrics(test_data$CRASH_TYPE, predictions$.pred_class)
print(paste("Precision:", metrics$precision))

## [1] "Precision: 0.884771341122852"

print(paste("Recall:", metrics$recall))

## [1] "Recall: 0.908086696819091"

print(paste("Accuracy:", metrics$accuracy))

## [1] "Accuracy: 0.844957142796734"