The Application of Logistic Regression to Examine the Predictors of Car Crashes Caused by Alcohol

Keywords: Logistic regression, Akaike Information Criterion (AIC), ROC curve, Area Under Curve (AUC), Sensitivity, Specificity, Misclassification Rate

GitHub Repository: MUSA5000-Logistic-Regression | Website

Introduction

Motor vehicle crashes remain one of the leading causes of injury and death in urban environments across the United States. Among these incidents, crashes involving alcohol-impaired drivers are particularly concerning because they account for a disproportionately high number of fatalities and severe injuries. According to national statistics, nearly one-third of all traffic-related deaths involve a driver under the influence of alcohol, underscoring the continued public health risks associated with impaired driving. In Philadelphia, where traffic density, roadway design, and neighborhood characteristics vary widely, it is especially important to understand the factors that contribute to alcohol-related crashes. By identifying these factors, city officials and public health professionals can design more targeted interventions to reduce the prevalence and severity of such incidents. This analysis uses crash data from Philadelphia covering the years 2008 to 2012 to investigate how various characteristics of crashes and neighborhoods are associated with the presence of a drinking driver.

This study explores the factors associated with DUI in traffic crashes in Philadelphia by examining a set of binary and continuous predictors. These include indicators of crash severity (Crashes resulted in fatality or major injury), risky driving behaviors (Crashes involved an overturned vehicle, speeding car and aggressive driving or Driver was using cell phone), driver demographics (Crashes involved at least one driver who was 16 or 17 years old and Crashes involved at least one driver who was at least 65 years old), and neighborhood-level socioeconomic characteristics (percentage with bachelor’s degree or more, Median household income). Among them, crashes resulting in fatalities or major injuries are more likely to involve alcohol, as intoxicated drivers have slower reaction times and poorer judgment. Behaviors such as speeding, aggressive driving, and overturning are more common in high-risk situations often associated with alcohol impairment. Distractions like cellphone use may not be directly linked to alcohol use, but they reflect a broader disregard for road safety that might overlap with impaired driving tendencies. Age also plays a role as younger drivers (ages 16–17) are less experienced and may be more prone to risk-taking, while older drivers (65+) typically show lower alcohol involvement because of health concerns or lifestyle differences. Beyond driver behavior, neighborhood-level characteristics such as education and income may also influence the likelihood of alcohol-involved crashes. Areas with higher percentages of college-educated residents may have greater awareness of drunk driving risks. Wealthier neighborhoods might also exhibit lower rates of alcohol-related crashes as they have access to greater availability of resources and support systems.

In this report, we use R to perform logistic regression analysis, examining the statistical relationships between the selected predictors and the likelihood that a crash involved a drinking driver across Philadelphia.

Method

Issues with OLS Regression

The OLS regression model can be expressed as following formula:

\[ y = \beta_0 + \beta_1x_1 + \beta_2x_2 + ... + \beta_kx_k + \epsilon \]

where:

• \(y\) is the dependent variable
• \(x_1, x_2, ..., x_k\) are the independent variables
• \(\beta_0, \beta_1, \beta_2, ..., \beta_k\) are the coefficients
• \(\epsilon\) is the error term

Using OLS regression to model binary outcomes can lead to several issues:

• First, in the OLS regression, \(\beta_1\) is interpreted as the amount of the dependent variable \(y\) changes when \(x_1\) increases by one unit, holding all other variables constant. However, when the dependent variable is binary, a one unit increase in \(x_1\) results in \(\beta_1\) increase in \(y\) no longer makes sense, as \(y\) can change only from 0 to 1 or from 1 to 0. This leads to some interpretation issues.
• Second, OLS regression assumes a linear relationship between the predictors and the dependent variable. With a binary dependent variable, this assumption no longer holds, as binary outcomes are inherently no linear.
• Third, OLS also assumes normally distributed residuals. However, with binary outcome, the residuals are not normally distributed. Instead, they follow a binomial distribution.
• Finally, OLS regression has another assumption of homoscedasticity, which means that the variance of the residuals is constant across all levels of the independent variables. However, with binary outcomes, this assumption is violated, as the variance of the residuals is not constant. It varies with the predicted probability, leading to heteroscedasticity.

All those violations of the OLS regression assumptions affect the validity of the model and the interpretation of the coefficients. In this way, we need to use a different approach to model binary outcomes.

Logistic Regression Concepts

Logistic regression provide a way for us to get around above-mentioned issues. Instead of predicting the the dependent variable \(y\) directly, logistic regression predicts the probability of the dependent variable \(y\) being equal to 1. This can be explained through the concept of odds. The odds of an event is a ratio that compares the likelihood of the event happening to the likelihood of it not happening. The odds of an event can be expressed as:

\[ \text{Odds} = \frac{\text{Probability of Event Happening}}{\text{Probability of Event Not Happening}} = \frac{p}{1 - p} \] Where \(p\) is the probability of the event happening. The odds can take any value from 0 to infinity, where 0 means the event will never happen and infinity means the event will always happen. For example, if there is a 60% of an event happening (probability\(p\)=0.6), the odds of the event would be:

\[ \text{Odds} = \frac{0.6}{1 - 0.6} = \frac{0.6}{0.4} = 1.5 \] This mean the vent is 1.5 times more likely to happen than not happening. In logistic regression, instead of predicting the probability of the event directly, we predict the log odds (also known as logit) of the outcome. The logit model with multiple predictors is expressed as:

\[ \log \left( \frac{p}{1 - p} \right) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_k x_k \] In our specific case, it can be written as: \[ \log \left( \frac{P(\text{DRINKING_D} = 1)}{1 - P(\text{DRINKING_D} = 1)} \right) = \beta_0 + \beta_1 \cdot \text{FATAL_OR_M} + \beta_2 \cdot \text{OVERTURNED} + \dots + \beta_9 \cdot \text{MEDHHINC} \]

Where:

• \(P(\text{DRINKING_D} = 1)\) is the probability of a crash involving a drinking driver.
• \(\log \left( \frac{P(\text{DRINKING_D} = 1)}{1 - P(\text{DRINKING_D} = 1)} \right)\) is the log-odds of drinking involvement.
• \(\beta_0\) is the intercept of the model, representing the log-odds of drinking and drive involved when all predictors are zero
• \(\beta_1, \beta_2, \dots, \beta_k\) are the coefficients of each predictor variable, indicating the change in log-odds for a one unit increase in each predictors, holdings other constant.

The predictor variable in the model is defined as follows:

• FATAL_OR_M: whether the crash involved a fatality
• OVERTURNED: whether the crash involved an overturned vehicle
• CELL_PHONE: whether the driver was using a cell phone
• SPEEDING: whether the crash involved speeding
• AGGRESSIVE: whether the crash involved aggressive driving
• DRIVER1617: whether the crash involved at least one driver aged 16 or 17
• DRIVER65PLUS: whether the crash involved at least one driver aged 65 or older
• PCTBACHMOR: percentage of the population with a bachelor’s degree or more
• MEDHHINC: median household income in the area

Logistic Function: The logistic function is the inverse of the logit function, as shown in following graph. It takes the linear predictor (log-odds) and transforms it into a probability value between 0 and 1 [0,1]. This transformation is ideal for modeling binary outcomes as it transforms the linear combination of predictors into a probability value.

knitr::include_graphics("logisticvslogit.png")

The logistic function is defined as:

\[ P(Y = 1) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_k x_k)}} \] In this specific case, it can be written as:

\[ P(\text{DRINKING_D} = 1) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 \cdot \text{FATAL_OR_M} + \dots + \beta_9 \cdot \text{MEDHHINC})}} \] As shown in the graph above, the logistic function has an S-shaped curve with the following properties:

• Outcome between 0 and 1 [0,1]: As the linear predictor approaches negative infinity, the probability approaches 0. As the linear predictor approaches positive infinity, the probability approaches 1. This means that the predicted probabilities will always be between 0 and 1, making it suitable for binary outcomes.
• Symmetric around 0.5: The logistic function is symmetric around the point where the linear predictor equals 0.5. This means that when the linear predictor is equal to 0, the predicted probability is 0.5, indicating an equal chance of the event occurring or not occurring. This midpoint is useful in binary classification as it helps interpret probabilities around 50%.

Overall, the logistic function is the type of translator function we are ;ooking for. It works well for binary dependent variable because it restricts the predicted probabilities to the range of 0 and 1, which is essential for binary outcomes. Each coefficient \(\beta_i\) in the logistic regression model represents the effect of a predictor on the log-odds of the outcome, which can be transformed to interpret effects on probability.

Logistic Regression Hypothesis

For each predictors in our logistic regression model, we test whether the predictor has a statistically significant effect on the outcome variable.

The Null Hypothesis (\(H_0\)), which states that the predictor has no effect on the outcome variable. In other words, the coefficient of the predictor is equal to zero (\(\beta_i = 0\)), writtern as:

\[ H_0 : \beta_i = 0 \] The Alternative Hypothesis (\(H_a\)), which states that the predictor has an effect on the outcome variable. In other words, the coefficient of the predictor is not equal to zero (\(\beta_i \neq 0\)), writtern as:

\[ H_0 : \beta_i = 0 \]

To test these hypotheses, we use the Wald test, which is a statistical test used to assess the significance of individual coefficients(\(\beta_i\)) in a regression model. The Wald test statistic is calculated as:

\[ W_i = \frac{\hat{\beta}_i}{\sigma(\hat{\beta}_i)} \]

Where:

• \(\hat{\beta}_i\) is the estimated coefficient for predictor \(i\)
• \(\sigma(\hat{\beta}_i)\) is the standard error of the estimated coefficient

Under the null hypothesis, the Wald test statistic follows a standard normal distribution \(N(0, 1)\). Large values of the Wald test statistic indicate that \(\beta_i\) is significantly different from zero, suggesting that the predictor has a significant effect on the outcome variable. We can safely reject the null hypothesis.

However, rather than looking at the estimated beta coefficients (\(\beta\)), most statisticans still prefer to looks at odds Ratios, which is calculated by exponentiating the estimated coefficients (\(\beta\)). For example, if the odds ratio for a predictor is 1.6, it means that the odds of the outcome are 1.6 times higher for a one-unit increase in the predictor, holding all other predictors constant. This makes it easier to interpret the effect of each predictor on the outcome variable.

In addition, an odds ratio greater than1 indicates that the predictor is associated with an increase likelihood of the outcome, while an odds ratio less than 1 indicates that the predictor is associated with a decrease likelihood of the outcome. An odds ratio equal to 1 indicates no effect.

Logistic Regression Model Assessment

R-squared

In logistic regression, there are several ways to access the quality of the model fit. An R-squared value can be computed for logistic regression, but it differs from the R-squared used in ordinary least squares (OLS) regression. In OLS, R-squared represents the proportion of variance in the outcome explained by the model, and a higher R-squared generally indicates better model fit. However, this interpretation does not apply to logistic regression. This is because logistic regression models a binary outcome and uses maximum likelihood estimation (MLE) rather than minimizing the sum of squared residuals. In another word, there is no direct concept of variance in the dependent variable to be explained in the same way. Therefore, R-squared values in logistic regression are not directly comparable to OLS R² and are rarely used as primary indicators of model fit. Notably, the glm() function in R does not report an R-squared value by default, which reflects its limited interpretability. If an R-squared value is desired, one option is to use the rms package in R, which provides a version of Nagelkerke R². While this adjusted R-squared can offer some insight into model fit, it should be interpreted with caution, and the only safe conclusion is higher values are better.

Akaike Information Criterion (AIC)

The Akaike Information Criterion (AIC) is commonly used for model comparison in logistic regression as an alternative of R-squared. AIC is an estimator of prediction error and the relative quality of statistical models for a given dataset. It provides a principled approach to comparing candidate models by balancing goodness of fit with model complexity, penalizing models with more parameters to reduce the risk of overfitting. AIC is defined as:

\[ \text{AIC} = 2k - 2\ln(\hat{L}) \]

where:

• \(k\) is the number of estimated parameters in the model
• \(\hat{L}\) is the maximum value of the likelihood function for the model

A lower AIC value suggests a better fit, as they balance the trade-off between model complexity and goodness of fit.

Specificity, Sensitivity and Misclassification Rate

In logistic regression, residuals (Errors) are defined as the difference between the observed and predicted values, similar to linear regression. However, in logistic regression, the predicted values are probabilities, not direct values of the dependent variable. Specifically, the predicted value \(\hat{y}_i\) is the estimated probability that \(Y = 1\) for the \(i\)-th observation.

The residuals in logistic regression are defined as:

\[ \epsilon_i = y_i - \hat{y}_i \]

Where:

• \(y_i\) is the actual observed value (either 0 or 1).

• \(\hat{y}_i\) is the predicted probability that \(Y = 1\) for the \(i\)-th observation.

The predicted probabilities \(\text{P(y=1)}\) =\(\hat{y}_i\) are calculated using the logistic function:

\[ \text{P(y=1)}=\hat{y}_i = \frac{e^{(\hat{\beta}_0 + \hat{\beta}_1 x_{1i} + \hat{\beta}_2 x_{2i} + \cdots + \hat{\beta}_p x_{pi})}}{1 + e^{(\hat{\beta}_0 + \hat{\beta}_1 x_{1i} + \hat{\beta}_2 x_{2i} + \cdots + \hat{\beta}_p x_{pi})}}= \frac{1}{1 + \exp(-\hat{\beta}_0 - \hat{\beta}_1x_1 - \dots - \hat{\beta}_kx_k)} \]

Where:

• \(\hat{y}_i\) is the estimated probability that \(Y = 1\) for the \(i\)-th observation.

• \(\hat{\beta}_0, \hat{\beta}_1, \dots, \hat{\beta}_p\) are the estimated coefficients obtained from the logistic regression model.

• \(x_{1i}, x_{2i}, \dots, x_{pi}\) are the values of the predictor variables for observation \(i\).

After obtaining the predicted probabilities, we apply a threshold (cut-off) value to classify predictions into binary outcomes (e.g. \(\hat{y}\)>0.5 as y= 1, and \(\hat{y}\)<0.5 as y=0). The choice of cut-off value can significantly impact the model’s performance. The default cut-off value is 0.5, but it can be adjusted based on the specific context and goals of the analysis.

To assess predictive performance and different cut-off values, we use metrics such as sensitivity, specificity and misclassification rate. These metrics are derived from the predicted values of \(\hat{y}\), which in logistic regression represent the estimated probability that \(Y = 1\).

• Sensitivity (also called recall or the true positive rate) measures the proportion of actual positives correctly predicted and is complementary to the false negative rate:

\[ \text{Sensitivity} = \frac{TP}{TP + FN} \]

• Higher sensitivity is better, especially when failing to detect true positives is costly.

• Specificity (the true negative rate) measures the proportion of actual negatives correctly predicted and is complementary to the false positive rate:

\[ \text{Specificity} = \frac{TN}{TN + FP} \]

• Higher specificity is better when false positives are particularly problematic.

• The correct classification rate is the proportion of all predictions that are correct:

\[ \text{Correct Classification Rate} = \frac{TP + TN}{\text{Total}} \]

• A higher correct classification rate indicates better model performance.

• The misclassification rate is the proportion of all predictions that are incorrect:

\[ \text{Misclassification Rate} = 1 - \text{Correct Classification Rate} \]

• Lower values of the misclassification rate indicate better performance, as it reflects the lower proportion of incorrect predictions.

Using different cut-off values allows for a trade-off between sensitivity and specificity. For example, a lower cut-off value may increase sensitivity but decrease specificity, while a higher cut-off value may increase specificity but decrease sensitivity. For example, if the cutoff value = 0.7, the model will classify a crash as involving a drinking driver only if the predicted probability is greater than 0.7. This may lead to fewer false positives (higher specificity) but more false negatives (lower sensitivity). Conversely, if the cut-off value = 0.3, the model will classify a crash as involving a drinking driver if the predicted probability is greater than 0.3, leading to more false positives (lower specificity) but fewer false negatives (higher sensitivity). The choice of cut-off value should be based on the specific context and goals of the analysis, considering the consequences of false positives and false negatives.

Receiver Operating Characteristic (ROC) Curve

One way to plot sensitivity (true positive rate) against \(1 -\) specificity (false positive rate) is through the Receiver Operating Characteristic (ROC) curve. The ROC curve is a graphical representation that helps assess the predictive quality of the model by comparing sensitivity and specificity across various thresholds. ROC Curves will be to the left of the 45 degree line.

A best cut-off value may be determined by optimizing the balance between sensitivity and specificity. There are a couple of ways to identify the optimal probability cut-off based on the ROC curve:

• We can calculate the Youden Index, which defined as:

\[ J = \text{Sensitivity} + \text{Specificity} - 1 \] The best cut-off can be found where the Youden Index is maximized.

• The second approaches involve calculate the Minimum Distance. A cut-off where the ROC curve has the minimum distance from the upper-left corner of the graph, where both specificity and sensitivity are equal to 1. In this case, we are minimizing the distance from the upper left corner, where \(y = \text{Sensitivity} = \text{True Positive Rate} = \text{tpr} = 1\) and \(x = 1 - \text{Specificity} = \text{False Positive Rate} = \text{fpr} = 0\). The formula used to calculate the distance is:

\[ d = (x - 0)^2 + (y - 1)^2 \]

For this report, we will be using the minimum distance approach to determine the optimal cut-off value.

Area Under the ROC Curve (AUC)

Based on the ROC curve, we can calculate the Area Under Curve (AUC) to quantify the model’s performance. Area Under the ROC Curve (AUC) means how well the model predicts \(Y = 1\) responses as 1’s and \(Y = 0\) responses as 0’s. Higher AUC values indicate that we can find a cut-off value for which both sensitivity and specificity are relatively high. The possible values for AUC range between 0.5 (Area under the 45-degree line, means no discrimination between classes, equivalent to random guessing) and 1 (Area of the entire box means perfect prediction). The AUC value may be interpreted as the probability that the \(\hat{y}\) for observation 1 (where \(y = 1\)) will be higher than the \(\hat{y}\) for observation 2 (where \(y = 0\)). A rough guide for classifying the accuracy based on AUC is:

• AUC = 0.90–1.00: Excellent
• AUC = 0.80–0.89: Good
• AUC = 0.70–0.79: Fair
• AUC = 0.60–0.69: Poor
• AUC < 0.60: Failing

We can generally consider an AUC value greater than 0.7 as acceptable.

Logistic Regression Assumptions

Logistic regression shares some assumptions with OLS regression, but it also has unique assumptions that accommodate the binary nature of the dependent variable. The key assumptions of logistic regression include:

• First, just like OLS regression, logistic regression assumes that observations are independent of each other. This means that each data point should not influence another, without spatial and time autocorrelation.
• Second, similar to OLS regression, logistic regression assumes that there is no severe multicollinearity among the independent variables. Multicollinearity occurs when two or more independent variables are highly correlated, making it difficult to determine the individual effect of each variable on the dependent variable and instability in the coefficient estimates.
• Third, unless the OLS regression, logistic regression does not assumes homoscedasticity. In OLS regression, the residuals should have constant variance across all levels of the independent variables (homoscedasticity). However, in logistic regression, this assumption is not required because the model uses a different error structure (binomial distribution).
• Forth, in OLS, we assume that the residuals are normally distributed. In logistic regression, however, the errors do not need to follow a normal distribution because the model uses a binomial distribution instead of normal distribution. This is one of the key differences between OLS and logistic regression.
• Lastly, logistic regression does not assume a linear relationship between the independent variables and the dependent variable.
• Instead, dependent variable for logistic regression must be binary (0 or 1). Logistic regression assumes that the log odds of the dependent variable being 1 is linearly related to the each independent variable. This is a unique assumptions of logistic regression that does not apply to OLS regression.
• In addition, logistic regression requires a sample size of at least 50 per predictors, as logistic regression uses Maximum Likelihood Estimation (MLE) to estimate the coefficients.

Exploratory Analysis

Before running logistic regression, we need to conduct exploratory analyses to better understand the relationship between the dependent variable (which is the drinking driver indicator DRINKING_D in the report) and potential predictors. These analyses help identify statistically significant associations and guide the selection of appropriate variables for the regression model.

For binary predictors, we examine the distribution across groups using cross-tabulations and assess significance using the Chi-Square test. For continuous predictors, we compare group means and standard deviations and use the independent samples t-test to determine whether the differences are statistically meaningful.

Binary Predictors

We first examine the relationship between the dependent variable and binary predictors through cross-tabulations. Cross-tabulations allow us to explore the distribution of observations across the different dependent variables.

To formally test whether a statistically significant association exists between the dependent variable and each binary predictor, we use the Chi-Square test. The Chi-Square statistic is derived from the sum of squared differences between observed and expected frequencies, normalized by the expected frequencies. The degrees of freedom (df) for the test can be calculated as \((R-1)(C-1)\), where \(R\) is the number of rows and \(C\) is the number of columns in the cross-tabulation table, representing the number of categories in each of the two variables being analyzed. The hypotheses are as follows:

• Null Hypothesis (H0): The proportions of different category of predictors are the same for drinking and non-drinking drivers. There is no association between the binary predictor and the dependent variable (DRINKING_D)

• Alternative Hypothesis (Ha): The proportions of different category of predictos differ between drinking and non-drinking drivers. There is an association between the binary predictor and the dependent variable (DRINKING_D)

If the p-value is below 0.05, we reject the null hypothesis and conclude that there is a statistically significant association between the two categorical variables.

Continuous Predictors

For continuous predictors, we first examine the continuous predictors by comparing their means and standard deviations for the two groups of the dependent variable (in this report, those with and without drinking drivers in crashes). This allows us to compare the distributions of the predictor across groups.

Then, we use the independent samples t-test to assess whether the observed differences in means are statistically significant. The hypotheses for the t-test are:

• Null Hypothesis (H0): The average values of the continuous predictor are the same for crashes that involve drinking drivers and those that do not. There is no significant difference in the means of the continuous predictor for crashes that involved alcohol and those that did not.

• Alternative Hypothesis (Ha): The average values of the continuous predictor are different for crashes that involve drinking drivers and those that do not. There is a significant difference in the means of the continuous predictor for crashes that involved alcohol and those that did not.

If the p-value of t-test is below 0.05, we reject the null hypothesis in favor of the alternative, indicating that the mean of the continuous predictor significantly differs between the two groups.

By performing these exploratory analyses, we can identify which predictors exhibit significant associations with the dependent variable and should be considered in the logistic regression model.

Results

Exploratory Analyses

Dependent Variable

In the report, the dependent variable is DRINKING_D, indicating whether a crash involved a drinking driver. Of the 53,260 recorded crashes in Philadelphia from 2008 to 2012, we excluded 9,896 crashes that occurred in non-residential block groups, where median household income and vacancy rates are 0. This leaves 43,364 crashes that took place in Philadelphia’s residential block groups.

depetab <- prop.table(table(crashes$DRINKING_D))
deptab <- as.data.frame(depetab)
colnames(deptab) <- c("DRINKING_D", "Proportion")
deptab %>%
  kable("html", escape = FALSE, align = "cc") %>%
  kable_styling(full_width = FALSE, position = "center", fixed_thead = TRUE)

DRINKING_D	Proportion
0	0.943
1	0.057

We see that there are 2,485 (5.73%) crashes that involved a drunk driver, and 40,879 (94.27%) crashes that did not involve a drunk driver. That is, the probability of a crash involving drinking drivers can be calculated using the formula:

\[ \text{Probability(Drinking Drivers)} = \frac{\text{Number of Crashes Involving Drinking Driver}}{\text{Total Number of Crashes}} = \frac{2485}{43364} \approx 0.0573 \]

The odds of a crash involving drinking drivers can be calculated using the formula:

\[ \text{Odds(Drinking Drivers)} = \frac{\text{Number of Crashes Involving Drinking Driver}}{\text{Number of Crashes Not Involving Drinking Drivers}} = \frac{2485}{40879} \approx 0.0608 \]

This suggests that while the overall probability of a crash being alcohol-related is relatively low, the odds still warrant attention, especially given the serious consequences often associated with such incidents.

Independent Variables - Binary

Then, let’s take a look at binary independent variables FATAL_OR_M,OVERTURNED, CELL_PHONE, SPEEDING, AGGRESSIVE, DRIVER1617, and DRIVER65PLUS.

# Load required library
library(gmodels)

# Ensure binary variables are factors
crashes$DRINKING_D <- as.factor(crashes$DRINKING_D)
binary_vars <- c("FATAL_OR_M", "OVERTURNED", "CELL_PHONE", "SPEEDING",
                 "AGGRESSIVE", "DRIVER1617", "DRIVER65PLUS")
crashes[binary_vars] <- lapply(crashes[binary_vars], as.factor)

# Define function to extract Chi-Square test results and percentages
extract_chi2_results <- function(dep_var, pred_var) {
  crosstab <- table(dep_var, pred_var)

  # Total counts per group
  total_no_alcohol <- sum(crosstab[1, ])
  total_alcohol <- sum(crosstab[2, ])

  # Count of "1"s
  no_alcohol_1_count <- crosstab[1, "1"]
  alcohol_1_count <- crosstab[2, "1"]

  # Percentages
  percent_no_alcohol <- 100 * crosstab[1, "1"] / total_no_alcohol
  percent_alcohol <- 100 * crosstab[2, "1"] / total_alcohol

  # Chi-square test without Yates correction
  chi2_test <- suppressWarnings(chisq.test(crosstab, correct = FALSE))

  # Output list
  return(list(
    crosstab = crosstab,
    no_alcohol_1_count = no_alcohol_1_count,
    alcohol_1_count = alcohol_1_count,
    percent_no_alcohol = percent_no_alcohol,
    percent_alcohol = percent_alcohol,
    total_n = sum(crosstab[, "1"]),
    chi2_stat = chi2_test$statistic,
    df = chi2_test$parameter,
    p_value = chi2_test$p.value
  ))
}

# Initialize results table
all_results <- data.frame(
  Variable = binary_vars,
  Description = c(
    "Crash resulted in fatality or major injury",
    "Crash involved an overturned vehicle",
    "Driver was using cell phone",
    "Crash involved speeding car",
    "Crash involved aggressive driving",
    "Crash involved at least one driver aged 16–17",
    "Crash involved at least one driver aged 65+"
  ),
  No_Alcohol_N = NA,
  No_Alcohol_Percent = NA,
  Alcohol_N = NA,
  Alcohol_Percent = NA,
  Total_N = NA,
  Chi2_Stat = NA,
  DF = NA,
  P_Value = NA
)

# Loop through variables and populate result table
for (i in seq_along(binary_vars)) {
  var <- binary_vars[i]
  result <- extract_chi2_results(crashes$DRINKING_D, crashes[[var]])

  all_results$No_Alcohol_N[i] <- result$no_alcohol_1_count
  all_results$No_Alcohol_Percent[i] <- round(result$percent_no_alcohol, 1)
  all_results$Alcohol_N[i] <- result$alcohol_1_count
  all_results$Alcohol_Percent[i] <- round(result$percent_alcohol, 1)
  all_results$Total_N[i] <- result$total_n
  all_results$Chi2_Stat[i] <- round(result$chi2_stat, 2)
  all_results$DF[i] <- result$df
  all_results$P_Value[i] <- format.pval(result$p_value, digits = 3, eps = 0.05)
}

## print

library(knitr)
library(kableExtra)

kable(all_results[, c("Variable", "Description",
                      "No_Alcohol_N", "No_Alcohol_Percent",
                      "Alcohol_N", "Alcohol_Percent",
                      "Total_N", "Chi2_Stat", "DF", "P_Value")],
      format = "html",
      col.names = c("Variable", "Description",
                    "N", "%", "N", "%", "N",
                    "Chi2 Stat", "DF", "χ² p-value"),
      caption = "Cross-Tabulation and Chi-Square Test Results") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = F,
                position = "center") %>%
  column_spec(1, bold = TRUE, width = "10em") %>%
  column_spec(2, bold = FALSE, width = "40em") %>%
  column_spec(3:7, width = "7em") %>%
  column_spec(8:10, width = "8em") %>%
  add_header_above(c(" " = 2,
                     "No Alcohol Involved (DRINKING_D = 0)" = 2,
                     "Alcohol Involved (DRINKING_D = 1)" = 2,
                     "Total" = 1,
                     "Chi-Square Test Results" = 3))

Cross-Tabulation and Chi-Square Test Results

		No Alcohol Involved (DRINKING_D = 0)		Alcohol Involved (DRINKING_D = 1)		Total	Chi-Square Test Results
Variable	Description	N	%	N	%	N	Chi2 Stat	DF	χ² p-value
FATAL_OR_M	Crash resulted in fatality or major injury	1181	2.9	188	7.6	1369	167.56	1	<0.05
OVERTURNED	Crash involved an overturned vehicle	612	1.5	110	4.4	722	122.79	1	<0.05
CELL_PHONE	Driver was using cell phone	426	1.0	28	1.1	454	0.16	1	0.687
SPEEDING	Crash involved speeding car	1261	3.1	260	10.5	1521	376.78	1	<0.05
AGGRESSIVE	Crash involved aggressive driving	18522	45.3	916	36.9	19438	67.60	1	<0.05
DRIVER1617	Crash involved at least one driver aged 16–17	674	1.6	12	0.5	686	20.45	1	<0.05
DRIVER65PLUS	Crash involved at least one driver aged 65+	4237	10.4	119	4.8	4356	80.60	1	<0.05

The table above provides a cross-tabulation of the dependent variable with each of the binary predictors. We use the Chi-Square test, which is a statistical method used to evaluate whether there is a significant association between two categorical variables. This test compares the observed frequencies in each category with the expected frequencies under the assumption that the variables are independent. As usual, a high value of the \(\chi^2\) statistic, and a p-value lower than 0.05, suggest that there is evidence to reject the null hypothesis in favor of the alternative, indicating an association between drunk driving and crash outcomes.

Here, we apply the Chi-Square test to explore the association between drinking driver indicators in crashes DRINKING_D and various binary crash-related characteristics. The hypotheses and results are summarized below:

Crash resulted in fatality or major injury FATAL_OR_M

• Null Hypothesis (H0): The proportions of crashes resulting in fatality or major injury are the same for drinking and non-drinking drivers.

• Alternative Hypothesis (Ha): The proportions of crashes resulting in fatality or major injury differ between drinking and non-drinking drivers.

Among crashes involving drinking drivers, 7.6% resulted in fatality or major injury, compared to only 2.9% of crashes involving non-drinking drivers. The P-value of Chi-Square statistic is less than 0.05, indicating a significant association between drinking drivers and crash severity. Therefore, we reject the null hypothesis in favor of the alternative hypothesis.

Crash involved an overturned vehicle OVERTURNED

• Null Hypothesis (H0): The proportions of crashes involving overturned vehicles are the same for drinking and non-drinking drivers.

• Alternative Hypothesis (Ha): The proportions of crashes involving overturned vehicles differ between drinking and non-drinking drivers.

Among crashes involving drinking drivers, 4.4% involved an overturned vehicle, compared to only 1.5% of crashes involving non-drinking drivers. The P-value of Chi-Square statistic is less than 0.05,, indicating a significant association between drinking drivers and overturned crashes. Therefore, we reject the null hypothesis in favor of the alternative hypothesis.

Driver was using cell phone CELL_PHONE

• Null Hypothesis (H0): The proportions of crashes involving cell phone use are the same for drinking and non-drinking drivers.

• Alternative Hypothesis (Ha): The proportions of crashes involving cell phone use differ between drinking and non-drinking drivers.

Among crashes involving drinking drivers, 1.0% involved cell phone use, compared to 1.1% of crashes involving non-drinking drivers. The P-value of Chi-Square statistic is more than 0.05,, suggesting no significant association with drinking driver indicators. Therefore, we fail to reject the null hypothesis.

Crash involved speeding car SPEEDING

• Null Hypothesis (H0): The proportions of crashes involving speeding cars are the same for drinking and non-drinking drivers.

• Alternative Hypothesis (Ha): The proportions of crashes involving speeding cars differ between drinking and non-drinking drivers.

Among crashes involving drinking drivers, 10.5% involved speeding, compared to only 3.1% of crashes involving non-drinking drivers. The P-value of Chi-Square statistic is less than 0.05,, indicating a strong link between drinking drivers and speeding. Therefore, we reject the null hypothesis in favor of the alternative hypothesis.

Crash involved aggressive driving AGGRESSIVE

• Null Hypothesis (H0): The proportions of crashes involving aggressive driving are the same for drinking and non-drinking drivers.

• Alternative Hypothesis (Ha): The proportions of crashes involving aggressive driving differ between drinking and non-drinking drivers.

Among crashes involving drinking drivers, 63.1% involved aggressive driving, compared to 54.7% of crashes involving non-drinking drivers. The P-value of Chi-Square statistic is less than 0.05,indicating he association is statistically significant. Therefore, we reject the null hypothesis in favor of the alternative hypothesis.

Crash involved at least one driver aged 16–17 DRIVER1617

• Null Hypothesis (H0): The proportions of crashes involving at least one driver aged 16-17 are the same for drinking and non-drinking drivers.

• Alternative Hypothesis (Ha): The proportions of crashes involving at least one driver aged 16-17 differ between drinking and non-drinking drivers.

Among crashes involving drinking drivers, 0.5% involved at least one teenage driver (16–17 years old), compared to 1.6% of crashes involving non-drinking drivers. The Chi-Square statistic is 19.71 (p < 0.001), indicating a significant association. Therefore, we reject the null hypothesis in favor of the alternative hypothesis.

Crash involved at least one driver aged 65 or older DRIVER65PLUS

• Null Hypothesis (H0): The proportions of crashes involving at least one driver aged 65 or older are the same for drinking and non-drinking drivers.

• Alternative Hypothesis (Ha): The proportions of crashes involving at least one driver aged 65 or older differ between drinking and non-drinking drivers.

Among crashes involving drinking drivers, 5.2% involved at least one driver aged 65 or older, compared to 10.4% of crashes involving non-drinking drivers. The Chi-Square statistic is 79.99 (p < 0.001), indicating a statistically significant relationship between alcohol involvement and older drivers.. Therefore, we reject the null hypothesis in favor of the alternative hypothesis.

The Chi-Square test results indicate significant associations between drinking driver indicators and most crash-related variables, with the exception of cell phone use, where no significant relationship was found. For all other variables, we reject the null hypothesis in favor of the alternative hypothesis, confirming the presence of significant associations.

Independent Variables - Continuous

Then, we examine the relationship between our binary dependent variable DRINKING_D and continuous predictors PCTBACHMOR and MEDHHINC.

continuous_vars <- c("PCTBACHMOR", "MEDHHINC")

results_list <- list()

for (var in continuous_vars) {
  means <- tapply(crashes[[var]], crashes$DRINKING_D, mean, na.rm = TRUE)
  sds   <- tapply(crashes[[var]], crashes$DRINKING_D, sd, na.rm = TRUE)
  t_test_result <- t.test(crashes[[var]] ~ crashes$DRINKING_D)

  description <- ifelse(var == "PCTBACHMOR", "Percentage of population with bachelor's degree or more",
                        ifelse(var == "MEDHHINC", "Median household income", NA))

  results_list[[var]] <- data.frame(
    Variable         = var,
    Description      = description,
    Mean_No_Alcohol  = means["0"],
    SD_No_Alcohol    = sds["0"],
    Mean_Alcohol     = means["1"],
    SD_Alcohol       = sds["1"],
    T_Statistic      = t_test_result$statistic,
    DF               = t_test_result$parameter,
    P_Value          = t_test_result$p.value
  )
}

continuous_summary <- do.call(rbind, results_list)
row.names(continuous_summary) <- NULL

# print
kable(continuous_summary[, c("Variable", "Description", "Mean_No_Alcohol", "SD_No_Alcohol",
                             "Mean_Alcohol", "SD_Alcohol",
                             "T_Statistic", "DF", "P_Value")],
      format = "html",
      col.names = c("Variable", "Description",
                    "Mean", "SD",
                    "Mean", "SD",
                    "t Stat", "DF", "p-value"),
      caption = "Summary of Continuous Variables") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = F,
                position = "center") %>%
  column_spec(1, bold = TRUE, width = "10em") %>%
  column_spec(2, bold = FALSE, width = "40em") %>%
  column_spec(3:4, width = "7em") %>%
  column_spec(5:6, width = "7em") %>%
  column_spec(7:9, width = "8em") %>%
  add_header_above(c(" " = 2,
                     "No Alcohol Involved (DRINKING_D = 0)" = 2,
                     "Alcohol Involved (DRINKING_D = 1)" = 2,
                     "T-Test Results" = 3))

Summary of Continuous Variables

		No Alcohol Involved (DRINKING_D = 0)		Alcohol Involved (DRINKING_D = 1)		T-Test Results
Variable	Description	Mean	SD	Mean	SD	t Stat	DF	p-value
PCTBACHMOR	Percentage of population with bachelor’s degree or more	16.6	18.2	16.6	18.7	-0.108	2778	0.914
MEDHHINC	Median household income	31483.1	16930.1	31998.8	17810.5	-1.405	2764	0.160

We first look at the means and standard deviations (SD) of these variables for crashes that involve alcohol and those that do not. Then, we use the independent samples t-test to assess whether there is a significant association between the dependent variable and each of the continuous predictors. The t-test results include the t-statistic, degrees of freedom, and p-value.

Percentage of Population with Bachelor’s Degree or More PCTBACHMOR

• Null Hypothesis (H0): The average values of the variable PCTBACHMOR are the same for crashes that involve drunk drivers and crashes that don’t.

• Alternative Hypothesis (Ha): The average values of the variable PCTBACHMOR are different for crashes that involve drunk drivers and crashes that don’t.

Since the p-value is much greater than the 0.05 significance threshold (p = 0.914), we fail to reject the null hypothesis. This suggests there is no significant difference in the mean percentage of the population with a bachelor’s degree or more between crashes involving drinking and non-drinking drivers.

Median Household Income MEDHHINC

• Null Hypothesis (H0): The average values of the variable MEDHHINC are the same for crashes that involve drunk drivers and crashes that don’t..

• Alternative Hypothesis (Ha): The average values of the variable MEDHHINC are different for crashes that involve drunk drivers and crashes that don’t.

Since the p-value is greater than 0.05 (p = 0.160), we fail to reject the null hypothesis. There is no significant difference in median household income between crashes involving drinking and non-drinking drivers.

For both continuous predictors, PCTBACHMOR and MEDHHINC, the p-values from the t-tests are above the 0.05 significance threshold, indicating that there is no statistically significant difference between the two groups (drinking and non-drinking drivers) for these continuous variables. Therefore, we fail to reject the null hypothesis in both cases, meaning there is no significant association between alcohol involvement in crashes and these continuous predictors.

Logistic Regression Assumptions

Based on the correlation matrix and exploitative analysis result, most logistic regression assumptions are met. First, The binary dependent variable , DRINKING_D, fulfills the requirement for a binary variable (0 and 1). Second, the independence of observation is satisfied, as each crash is considered an independent event. Third, the assumption of no perfect multicollinearity is also met, as the correlation matrix shows low pairwise pearson correlations between predictors, with only the highest correlation being 0.48 between percent with bachelor’s degree and median household income. These correlation coefficients are all below 0.5, indicating that there is no perfect multicollinearity among the predictors. These correlation are not high enough to cause instability in the model. Additionally, the large sample size of over 43,300 observations ensure the model meet the assumption for sufficient sample size.

However, there is a potential violation in the assumption of linearity of log odds for continuous predictors, as high p-values in t-test suggest weal or non-linear relationships between percent with bachelor’s degree and median household income and the dependent variable.

Using Pearson correlation to measure associations between binary predictors have several limitations, as Pearson’s correlation is most appropriate for continuous normally distributed variables. For binary variables which only has two values (0 and 1), Person correlation can lead to less meaningful or interpretable correlation values, as the measure does not fully capture the association between two binary variables. Furthermore, Pearson correlation may underestimate or fail to reflect non-linear relationships that could potentially exist among binary variable.

Multicollinearity occurs when two or more predictor variables are highly correlated, potentially inflating standard errors and destabilizing coefficient estimates. An examination of our dataset’s correlation matrix reveals only low pairwise correlations with the highest of 0.48 between the percentage of residents with a bachelor’s degree and median household income. This value falls well below the conventional threshold of 0.7, indicating that our model is unlikely to suffer from coefficient instability due to collinear predictors. While Pearson’s correlation is not the optimal metric for binary variables, it nonetheless confirms the absence of any problematic multicollinearity in this analysis.

crashes <- crashes %>%
  mutate(across(c(FATAL_OR_M, OVERTURNED, CELL_PHONE, SPEEDING,
                  AGGRESSIVE, DRIVER1617, DRIVER65PLUS),
                ~ as.numeric(as.character(.))))

cor_matrix <- cor(crashes %>% dplyr::select(-c(CRN, DRINKING_D, AREAKEY, COLLISION_)), use = "pairwise.complete.obs", method = "pearson")


custom_label <- c(
  "Crash resulted in fatality or major injury" = "FATAL_OR_M",
  "Crash involved an overturned vehicle" = "OVERTURNED",
  "Driver was using cell phone" = "CELL_PHONE",
  "Crash involved speeding car" = "SPEEDING",
  "Crash involved aggressive driving" = "AGGRESSIVE",
  "Crash involved at least one driver 16 or 17 years old" = "DRIVER1617",
  "Crash involved at least one driver over 65 years old" = "DRIVER65PLUS",
  "% with bachelor’s degree" = "PCTBACHMOR",
  "Median household income" = "MEDHHINC"
)


rownames(cor_matrix) <- names(custom_label)
colnames(cor_matrix) <- names(custom_label)

ggcorrplot(cor_matrix,
           type = "lower",
           lab = TRUE,
           lab_size = 3,
           colors = c("#001219", "white", "#9b2226")) +
  labs(title = "Correlation Matrix for all Predictor Variables") +
  theme(plot.subtitle = element_text(size = 9, face = "italic"),
        plot.title = element_text(size = 12, face = "bold"),
        axis.text.x = element_text(size = 7),
        axis.text.y = element_text(size = 7),
        axis.title = element_text(size = 8))

logistic regression results

Our first logistic regression model includes all the predictors and presented the results of it below:

• FATAL_OR_M: The presence of a fatality or major injury significant increases the odds of a drinking related crash, with an odds ratio of 2.257 (95% confident interval: 1.91 to 2.65) and a highly significant p-value (p < 0.001). This suggests that these types of crashes are more than double the odds of being associated with drinking drivers.
• OVERTURNED: Overturned vehicles significantly increase the likelihood of it being a drinking related crash with and odd ratio of 2.532 (95% confident interval: 2.04 to 3.12) and a highly significant p-value (p < 0.001). This result implies that overturned vehicles are strongly associated with a drinking-driver incurred crashes.
• CELL_PHONE: The use of celphone during a crash shows no statistically significent effect on the odds of a drinking-related crash, as the p-value is greater than 0.05 (P=0.881). The odd ratio close to 1 (OR=1.03, 95% confident interval: 0.68 to 1.49) indicates that use of cellphone no significantly meaning affect the odds of a drinking-related crash.
• SPEEDING: Speeding is a highly significant predictor, with an odd ratio of 4.66 (95% confident interval: 3.97 to 5.45) and a p-value less than 0.001. This indicates that crashes involving speeding are more than four times as likely to involve alcohol.
• AGGRESSIVE: Aggressive driving is associated with lower odds of drinking-related crashes, with an odds ratio of 0.551 (95% confident interval: 0.50 to 0.60) and a highly significant p-value (p < 0.001). This suggests that aggressive driving is less likely to be associated with drinking drivers.
• DRIVER1617: Incidents involving drivers aged 16-17 are significantly less likely to associate with drinking-related crashes, with an odds ratio of 0.278 (95% confident interval: 0.148 to 0.471) and a highly significant p-value (p < 0.001). This probably caused by the fact that the drinking age in the US is 21, and drivers aged 16-17 are not legally allowed to drink alcohol.
• DRIVER65PLUS: Incidents involving drivers aged 65 or older are also show significantly lower odds of drinking-related, with an odds ratio of 0.461 (95% confident interval: 0.38 to 0.55) and a highly significant p-value (p < 0.001). This suggests that older drivers are less likely to be involved in drinking-related crashes.
• PCTBACHMOR: The percentage of the population with a bachelor’s degree or more shows no significant effect on the odds of a drinking-related crash, with an odds ratio of 1.000 (95% confident interval: 0.99 to 1.00) and a p-value greater than 0.05 (p = 0.775). This suggests that the educational attainment of the population does not significantly affect the likelihood of drinking-related crashes.
• MEDHHINC: The median household income shows has a very slight statistically significant effect on the likelihood of drinking-related crash, with p-value <0.05 (p= 0.036) and an odds ratio of 1.000 (95% confident interval: 1.00 to 1.00). This suggests that higher median household income is associated with a very slight negligible effect in the odds of drinking-related crashes.

In summary, FATAL_OR_M, OVERTURNED, SPEEDING, AGGRESSIVE, DRIVER1617, and DRIVER65PLUS are the most significant predictors of drinking-related crashes (P<0.001), while CELL_PHONE and PCTBACHMOR show no significant effect. Among those, FATAL_OR_M, OVERTURNED, and SPEEDING are the most significant predictors associated with drinking-related crashes, with odds ratios greater than 2.5. The predictor MEDHHINC is marginally significant and has minimal effect on the odds of drinking-related crashes.

logit <- glm(DRINKING_D ~ FATAL_OR_M + OVERTURNED + CELL_PHONE + SPEEDING + AGGRESSIVE + DRIVER1617 + DRIVER65PLUS + PCTBACHMOR + MEDHHINC, data = crashes, family = binomial)

logitoutput <- summary(logit)
logitcoeff <- logitoutput$coefficients
or_ci <- exp(cbind(OR = coef(logit), confint(logit)))

## Waiting for profiling to be done...

logitcoeff[, "Pr(>|z|)"] <- round(logitcoeff[, "Pr(>|z|)"], 4)
final_output <- cbind(logitcoeff, or_ci)
final_output %>%
kable("html", escape = FALSE, align = "cccccc") %>%
kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width        = FALSE,
    position          = "center"
  )

	Estimate	Std. Error	z value	Pr(>|z|)	OR	2.5 %	97.5 %
(Intercept)	-2.733	0.046	-59.563	0.000	0.065	0.059	0.071
FATAL_OR_M	0.814	0.084	9.713	0.000	2.257	1.910	2.653
OVERTURNED	0.929	0.109	8.509	0.000	2.532	2.035	3.122
CELL_PHONE	0.030	0.198	0.149	0.881	1.030	0.684	1.488
SPEEDING	1.539	0.081	19.107	0.000	4.660	3.974	5.450
AGGRESSIVE	-0.597	0.048	-12.493	0.000	0.551	0.501	0.604
DRIVER1617	-1.280	0.293	-4.367	0.000	0.278	0.148	0.471
DRIVER65PLUS	-0.775	0.096	-8.081	0.000	0.461	0.380	0.553
PCTBACHMOR	0.000	0.001	-0.286	0.775	1.000	0.997	1.002
MEDHHINC	0.000	0.000	2.091	0.036	1.000	1.000	1.000

Then, we test different cutoff threshold to see the model overall performance through the specificity, sensitivity, and misclassficiantion rate. The table below presents the specificity, sensitivity, and misclassification rates for various probability cut-offs in predicting the likelihood of a drinking-related crash.

The probability cut-off value of 0.50 yields the lowest misclassfication rate of 0.057. At this threshold, specificity is maximized at 1, but sensitivity is relative low at 0.02. This indicate the model capture relatively few true positive cases of drinking-related crashes, but it is very good at identifying non-drinking related crashes.

The highest misclassfication rate occurs at the cut-off of 0.02 with a rate of 0.889. At this threshold, sensitivity is high at 0.899, capturing a high proportion of true positive cases, but specificity is very low at 0.058. This indicates that the model misclassifies a large number of non-drinking related crashes as drinking-related crashes.

cutoffs <- c(0.02, 0.03, 0.05, 0.07, 0.08, 0.09, 0.10, 0.15, 0.20, 0.50)
results <- data.frame(
  Cutoff = cutoffs,
  Sensitivity = NA_real_,
  Specificity = NA_real_,
  MisclassificationRate = NA_real_
)

#–– 3. Pull out true labels and predicted probabilities
labels      <- crashes$DRINKING_D
predictions <- fitted(logit)   # or logit$fitted.values

#–– 4. Metric‐calculation function
calculate_metrics <- function(cutoff, labels, preds) {
  preds_bin <- ifelse(preds >= cutoff, 1, 0)
  tp  <- sum(preds_bin == 1 & labels == 1)
  tn  <- sum(preds_bin == 0 & labels == 0)
  fp  <- sum(preds_bin == 1 & labels == 0)
  fn  <- sum(preds_bin == 0 & labels == 1)
  sens <- tp / (tp + fn)
  spec <- tn / (tn + fp)
  mis  <- (fp + fn) / length(labels)
  c(Sensitivity = sens,
    Specificity = spec,
    MisclassificationRate = mis)
}

#–– 5. Populate results table
for (i in seq_along(cutoffs)) {
  results[i, 2:4] <- calculate_metrics(cutoffs[i], labels, predictions)
}

#–– 6. Render as a styled kable
results %>%
  kable(
    caption = "Model Performance Metrics Across Probability Cutoffs",
    digits  = 3,
    align   = c("c","c","c","c")
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width        = FALSE,
    position          = "center"
  )

Model Performance Metrics Across Probability Cutoffs

Cutoff	Sensitivity	Specificity	MisclassificationRate
0.02	0.984	0.058	0.889
0.03	0.981	0.064	0.884
0.05	0.735	0.469	0.516
0.07	0.221	0.914	0.126
0.08	0.185	0.939	0.105
0.09	0.168	0.946	0.099
0.10	0.164	0.948	0.097
0.15	0.104	0.972	0.078
0.20	0.023	0.995	0.060
0.50	0.002	1.000	0.057

The ROC curve below show the trade of between sensitivity and specificity across various cut-offs threshold for the logistic regression model. An ideal model would reach the top-left corner of the ROC space, where sensitivity and specificity all equal to one, which would mean perfect classifications. Our ROC curve lies beneath the ideal point, indicating trade-offs between sensitivity and specificity at each cut-off.

ggplot(roc, aes(d = labels, m = predictions)) +
  geom_roc(n.cuts = 50, labels = FALSE, colour = "#283d3b") +
  labs(title = "ROC Curve") +
  style_roc(theme = theme_grey) +
  geom_abline(slope = 1, intercept = 0, size = 1, color = "#c44536") +
  theme(
    plot.subtitle = element_text(size = 9, face = "italic"),
    plot.title = element_text(size = 12, face = "bold"),
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    axis.title = element_text(size = 9),
    panel.background = element_blank(),
    panel.border = element_rect(colour = "grey", fill = NA, linewidth = 0.8)
  )

We calculated the optimal cut-off by minimizing the Euclidean distance from the ROC curve to the top-left corner of the ROC space, which represents the ideal point where both sensitivity and specificity are equal to 1. This cut-off of 0.064 provides a balance approaches between specificity and sensitivity, with a sensitivity of 0.661 and specificity of 0.545.

pred <- prediction(roc$predictions, roc$labels)
roc.perf <- performance(pred, measure = "tpr", x.measure = "fpr")

opt.cut <- function(perf, pred) {
  cut.ind <- mapply(FUN = function(x, y, p) {
    d <- (x - 0)^2 + (y - 1)^2
    ind <- which(d == min(d))
    c(sensitivity = y[[ind]], specificity = 1 - x[[ind]], cutoff = p[[ind]])
  }, perf@x.values, perf@y.values, pred@cutoffs)

  cut.df <- data.frame(
    sensitivity = cut.ind[1, ],
    specificity = cut.ind[2, ],
    cutoff = cut.ind[3, ]
  )
  rownames(cut.df) <- NULL
  return(cut.df)
}

optimal_cutoffs_df <- opt.cut(roc.perf, pred)
optimal_cutoffs_df %>%
  kable("html", escape = FALSE, align = "ccc") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width        = FALSE,
    position          = "center"
  )

sensitivity	specificity	cutoff
0.661	0.545	0.064

In the above section, the optimal cut-off based on the minimum misclassification rate was 0.50, which yielded a misclassfication rate of 0.057. This threshold focused more on prioritizing specificity (1.00), but resulting in low sensitivity (0.02). In contrast, the optimal cut-off based on the ROC curve was 0.064, which provided a balance between sensitivity (0.661) and specificity (0.545). This cut-off is more suitable for practical applications where both false positives and false negatives are important to consider.

The AREA Under the ROC Curve (AUC) for our model is 0.640 suggests that the model has poor discrimination based on the metrix provided in the method section. The AUC value mean that if we randomly select one positive case and negative case, the model has a 64% chance of correctly ranking the positive case higher than the negative case. In practical, this indicates that while the model is able to distinguish between crashes with and without alcohol involvement to some extent, its performance is limited and may misclassify a substantial number of observations.

auc.perf <-  performance(pred, measure ="auc")
auc.value <- auc.perf@y.values
cat(sprintf("AUC: %.3f\n", auc.value))

## AUC: 0.640

The second logistic model we ran only includes the binary predictors. PCTBACHMOR and MEDHHINC were not included as they were not statistically significant predictors in the first model. However, in the second model, the exclusion of these two predictors did not lead to any predictors becoming significant or any previous significant predictors becoming non-significant. The odds ratios (OR) for the binary predictors in both models are very similar, suggesting that removing PCTBACHMOR and MEDHHINC had little impact on the estimated effects of the other variables.

logit2 <- glm(DRINKING_D ~ FATAL_OR_M + OVERTURNED + CELL_PHONE + SPEEDING + AGGRESSIVE + DRIVER1617 + DRIVER65PLUS, data = crashes, family = binomial)

logitoutput2 <- summary(logit2)
logitcoeff2 <- logitoutput2$coefficients
or_ci2 <- exp(cbind(OR = coef(logit2), confint(logit2)))

## Waiting for profiling to be done...

logitcoeff2[, "Pr(>|z|)"] <- round(logitcoeff2[, "Pr(>|z|)"], 4)
final_output2 <- cbind(logitcoeff2, or_ci2)
final_output2 %>%
  kable("html", escape = FALSE, align = "cccccc") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width        = FALSE,
    position          = "center"
  )

	Estimate	Std. Error	z value	Pr(>|z|)	OR	2.5 %	97.5 %
(Intercept)	-2.652	0.028	-96.324	0.000	0.071	0.067	0.074
FATAL_OR_M	0.809	0.084	9.662	0.000	2.246	1.901	2.640
OVERTURNED	0.940	0.109	8.619	0.000	2.559	2.057	3.156
CELL_PHONE	0.031	0.198	0.157	0.875	1.032	0.685	1.491
SPEEDING	1.540	0.081	19.128	0.000	4.666	3.980	5.457
AGGRESSIVE	-0.594	0.048	-12.433	0.000	0.552	0.503	0.606
DRIVER1617	-1.272	0.293	-4.338	0.000	0.280	0.149	0.475
DRIVER65PLUS	-0.766	0.096	-8.004	0.000	0.465	0.383	0.558

The Akaike Information Criterion (AIC) for both models is 18360. AIC is a measure of model quality that balances goodness of fit with model complexity. Since the AIC is the same for both models, it suggests that the second model (with only binary predictors) is just as good as the first model (with all predictors) in terms of explaining the data. This indicates that the additional complexity of including PCTBACHMOR and MEDHHINC in the first model does not provide a significant improvement in model fit, and thus, the simpler model with only binary predictors is preferred.

aic_table <- AIC(logit, logit2)
colnames(aic_table) <- c("Degree of Freedom", "AIC")
rownames(aic_table) <- c("Model 1: All Predictors", "Model 2: Binary Predictors Only")
aic_table %>%
  kable("html", escape = FALSE, align = "cc") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width        = FALSE,
    position          = "center"
  )

	Degree of Freedom	AIC
Model 1: All Predictors	10	18360
Model 2: Binary Predictors Only	8	18360

Discussion

In this paper, we examined over 43,000 motor vehicle crashes in Philadelphia from 2008 to 2012 to identify which factors were associated with alcohol involvement using logistic regression. We explored a mix of binary crash-level indicators and neighborhood-level continuous predictors to assess their relationship with the likelihood that a crash involved a drinking driver.

Several variables emerged as strong predictors of alcohol-involved crashes. Crashes that resulted in fatalities or major injuries, involved overturned vehicles, or involved speeding were significantly more likely to include a drinking driver. On the other hand, the presence of aggressive driving, teenage drivers, and elderly drivers were all significantly associated with lower odds of alcohol involvement. Cell phone use showed no significant association, nor did the neighborhood-level variables—percentage of residents with a bachelor’s degree and median household income.

While many results aligned with expectations—particularly the strong links between alcohol involvement and speeding or severe crashes—the lack of significance for socioeconomic predictors was somewhat unexpected. It suggests that while neighborhood context is often important in shaping health and safety outcomes, it may play a more limited role in predicting whether an individual crash involved alcohol. Similarly, the inverse relationship between aggressive driving and alcohol involvement may reflect different behavioral risk profiles, where alcohol impairment substitutes, rather than coincides with, certain forms of recklessness.

Logistic regression was an appropriate modeling choice given the binary outcome and the large, independent dataset. However, since only 5.7% of crashes involved a drinking driver, the problem qualifies as a rare event scenario. In such cases, alternative approaches like rare event logistic regression, as proposed by Paul Allison, may offer improved parameter estimates and better classification performance, especially for policy applications that require accurate identification of low-probability events.

This analysis has several limitations. The dataset does not include contextual factors such as time of day, proximity to bars, or day of the week, which could provide additional explanatory power. Pearson correlation was used to assess multicollinearity, but it is not ideal for binary variables and may overlook more complex interactions. Although multicollinearity was not a concern, the model’s AUC of 0.64 indicates limited predictive ability—adequate for identifying general associations but insufficient for high-stakes prediction. Lastly, underreporting of alcohol use or misclassification in crash records could affect the reliability of the outcome variable.