Predicting Car Crashes Caused by Alcohol Using Logistic Regression
Emily Zhou, Ziyi Guo, Emma Jiang
2024-11-11
Version 2.0 | First Created Nov 4, 2024 | Updated Nov 11, 2024
Keywords: t-test, chi-square, multiple logistic regression, ROC, sensitivity, specificity, odds ratio
GitHub Repository: CPLN671-Logit-Regression
Introduction
Alcohol-impaired driving is a significant public safety issue in the United States, with nearly 30 fatalities each day attributed to it, along with countless injuries and an annual economic burden exceeding $59 billion. Understand the back reason of drunk driving, and what factors can be associated with drunk driving is essential to tackling this problem and ease the safety issue. Here we are using Philadelphia as our study area.
In this analysis, we aim to identify predictors associated with drunk driving. Here, our dependent variable is whether or not the driver was drunk, and we use several independent variables to predict this outcome. These predictors include factors such as whether the crash resulted in a fatality or major injury, involved an overturned vehicle, or involved a driver using a cell phone. Additional behavioral factors like speeding, aggressive driving, and the presence of an immature driver (defined as at least one driver aged 16 or 17) or a senior driver (defined as at least one driver aged 65 or older) are also included. Furthermore, socio-economic indicators from the crash location, such as the percentage of individuals with at least a bachelor’s degree and the median household income within the census block group, are considered. By analyzing these variables, we aim to better understand the factors that correlate with drunk driving incidents.
Some of the predictors in this analysis include crashes involving fatalities or major injuries, overturned vehicles, cell phone use, speeding, and aggressive driving. These factors may be associated with drunk driving because they reflect behaviors and outcomes often connected to the impaired judgment and slowed reaction times that result from alcohol consumption. For instance, alcohol impairment increases the likelihood of severe crashes resulting in fatalities or major injuries, as impaired drivers are less able to respond to sudden changes on the road. Overturned vehicles are also more common in drunk-driving incidents, as alcohol can lead to a loss of vehicle control. Additionally, cell phone use, which is already a risky behavior, becomes even more dangerous when combined with alcohol impairment, as drunk drivers may be more prone to engage in distractions. Speeding is another common factor, as alcohol reduces inhibitions, leading drivers to ignore speed limits and traffic laws. Aggressive driving is also heightened by alcohol, which can increase aggression and lead to risky behaviors like tailgating or cutting off other vehicles. Altogether, these predictors capture patterns of behavior that are more likely when alcohol impairs a driver’s mental and physical abilities, leading to a greater likelihood of severe accidents.
These predictors, \(\text{DRIVER1617}\) (crash involving a driver aged 16 or 17) and \(\text{DRIVER65PLUS}\) (crash involving a driver aged 65 or older), might be associated with drunk driving due to age-related factors in behavior and ability. Younger drivers, especially those aged 16 or 17, may be more prone to risky behaviors, including experimenting with alcohol, due to lower experience levels and a tendency toward risk-taking. This age group may lack the experience and judgment to handle alcohol’s effects, making them more susceptible to involvement in drunk-driving crashes. On the other hand, older drivers, represented by DRIVER65PLUS, may also be affected by alcohol impairment differently. While they may not engage in risky behavior as frequently, alcohol can significantly impair their reaction times and decision-making abilities, which may already be slower due to age. This age group might also be more vulnerable to the effects of even small amounts of alcohol, leading to a higher risk of accidents when impaired. Together, these predictors help capture how different age groups may be uniquely vulnerable to the risks associated with drunk driving.
The predictors \(\text{PCTBACHMOR}\) (the percentage of individuals aged 25 or older with at least a bachelor’s degree) and \(\text{MEDHHINC}\) (median household income) might be associated with drunk driving as they reflect socio-economic characteristics of the area where crashes occur. Areas with higher education levels, indicated by a greater percentage of individuals with bachelor’s degrees, may have lower instances of drunk driving due to increased awareness of the risks associated with impaired driving and potentially stronger adherence to safety practices. Conversely, areas with lower education levels might see higher rates of alcohol-related incidents due to various socio-economic stressors or a lack of access to alternative transportation options. Similarly, median household income can play a role in influencing drunk driving patterns. Lower-income areas may experience higher rates of alcohol-impaired driving if public transportation options are limited, making driving the primary mode of transport regardless of impairment. On the other hand, higher-income areas might have better access to ridesharing services or designated driver options, potentially reducing the likelihood of drunk driving. By examining these socio-economic indicators, we can gain insight into how community characteristics might influence the prevalence of drunk driving incidents.
To conduct this analysis, we will use logistic regression in R to evaluate the relationship between these predictors and the likelihood of drunk driving.
Methods
Issues with OLS
Using Ordinary Least Squares (OLS) regression for a binary dependent variable poses several issues. First, recall below the equation of OLS regression:
\[ y = \beta_0 + \beta_1x_1 + \beta_2x_2 + ... + \beta_kx_k + \epsilon \]
Here, \(\beta_1\) is interpreted as the amount by which the dependent variable \(y\) changes for a one-unit change in the independent variable \(x_1\), 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 if we use OLS regression for binary dependent variables.
Second, OLS regression assumes a linear relationship between the independent variables and the dependent variable. With a binary dependent variable, this assumption does not hold, as binary outcomes (0 or 1) are inherently nonlinear. Third, OLS assumes normally distributed residuals (errors), but with a binary outcome, the residuals are not normally distributed. Instead, they follow a binomial distribution. This makes OLS estimation inefficient, as it cannot account for the specific distribution of binary outcomes. On top of these, with binary data, the variance of the error term is not constant; it varies with the predicted probability, leading to heteroscedasticity. This affects the reliability of standard errors and confidence intervals in OLS regression, leading to potentially incorrect inferences. The assumpstions that apply to logistic regression will be elaborated below.
Logistic Regression Concepts
Logistic regression provides a way for us to get around the above-mentioned issues. Instead of predicting the \(y\) value, we predict the probability of \(y\) being 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. Mathematically, odds are calculated as:
\[ \text{Odds} = \frac{\text{Probability of Event Happening}}{\text{Probability of Event Not Happening}} = \frac{p}{1 - p} \] For example, if there is a 70% chance of an event happening (probability \(p = 0.7\)), the odds would be:
\[ \text{Odds} = \frac{0.7}{1 - 0.7} = \frac{0.7}{0.3} = 2.33 \]
This means the event is 2.33 times more likely to happen than not happen. In logistic regression, instead of predicting the probability directly, we predict the log-odds (also called the logit) of the outcome. The logit model with multiple predictors is given by:
\[ \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 that drinking was involved in a crash. \(\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, representing the log-odds of drinking being involved when all predictors are zero. \(\beta_1, \beta_2, \dots, \beta_9\) are the coefficients for each predictor variable, indicating the change in log-odds for a one-unit increase in each predictor.
The predictor variable in the model is defined as follows:
FATAL_OR_M: Indicator for whether the crash involved a fatality.
OVERTURNED: Indicator for whether the crash involved an overturned vehicle.
CELL_PHONE: Indicator for whether a cell phone was in use.
SPEEDING: Indicator for whether speeding was a factor in the crash.
AGGRESSIVE: Indicator for whether aggressive driving was involved.
DRIVER1617: Indicator for whether the driver was between 16-17 years old.
DRIVER65PLUS: Indicator for whether the driver was 65 years or older.
PCTBACHMOR: Percentage of people in the area with a bachelor’s degree or more.
MEDHHINC: Median household income in the area.
The Logistic Function: The logistic function is the inverse of the logit function. It takes the linear predictor (log-odds) and converts it back to a probability, ensuring that the output is within the range \([0, 1]\). This transformation is ideal for modeling binary outcomes because it naturally maps any linear predictor onto a probability scale. The logistic function is given by:
\[ P(\text{DRINKING_D} = 1) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 \cdot \text{FATAL_OR_M} + \dots + \beta_9 \cdot \text{MEDHHINC})}} \]
The logistic function is the type of translator function we are looking for. It has an S-shaped (sigmoid) curve with the following properties:
Bounded Output: As the linear predictor approaches \(+\infty\), the probability \(P(\text{DRINKING\_D} = 1)\) approaches 1. As the predictor approaches \(-\infty\), the probability approaches 0, ensuring the output is within the valid probability range \([0, 1]\).
Symmetry: The logistic function is symmetric around 0.5, meaning that when the linear predictor is 0 (i.e., \(\beta_0 + \beta_1 \cdot \text{FATAL_OR_M} + \dots + \beta_9 \cdot \text{MEDHHINC} = 0\)), the probability is 0.5. This midpoint is useful in binary classification as it helps interpret probabilities around 50%.
Overall, it works well for binary models because it restricts predicted probabilities to a range of 0 and 1, which is required in binary classification tasks.Each coefficient \(\beta_i\) represents the effect of a predictor on the log-odds of the outcome, which can be transformed to interpret effects on probabilities.
Logistic Regression Hypothesis
For each predictor in our logistic regression model, we want to test whether the predictor has a statistically significant effect on the likelihood of the outcome.
The Null Hypothesis (\(H_0\)), which states that the predictor has no effect on the outcome, can be written as:
\[ H_0 : \beta_i = 0 \]
The Alternative Hypothesis (\(H_a\)) which states that the predictor has significant effect on the outcome, can be written as:
\[ H_a : \beta_i \neq 0 \]
To test these hypotheses, we use the Wald statistic for each coefficient \(\beta_i\). The Wald statistic is calculated as:
\[ W_i = \frac{\hat{\beta}_i}{\sigma(\hat{\beta}_i)} \]
where \(\hat{\beta}_i\) is the estimated coefficient for predictor, and \(\sigma(\hat{\beta}_i)\) is the standard error of this estimate. Under the null hypothesis, the Wald statistic approximately follows a standard normal distribution \(N(0, 1)\). Large values of \(|W_i|\) indicate that \(\beta_i\) is significantly different from zero, leading us to reject the null hypothesis.
However, rather than looking at the estimated \(\beta\) coefficients, most statisticians still prefer to look at odds ratios. The odds ratio is calculated by exponentiating the coefficients. For example, if the odds ratio for a predictor is 1.5, it means that the odds of the outcome happening are 1.5 times higher for a one-unit increase in the predictor. An odds ratio greater than 1 indicates that the predictor is associated with an increased likelihood of the outcome, while an odds ratio less than 1 indicates a decreased likelihood.
Logistic Regression Model Assessment
In logistic regression, there are several ways to assess the quality of the model fit. While \(R^2\) is commonly used in OLS to measure the proportion of variance explained by the model, it is not directly applicable to logistic regression. Although a pseudo-\(R^2\) statistic can be calculated for logistic regression (e.g., McFadden’s \(R^2\)), it does not have the same interpretation as in OLS regression and is not as widely used for model assessment.
A better alternative for comparing models is the Akaike Information Criterion (AIC). The AIC is a measure that penalizes models for complexity (i.e., the number of predictors) and rewards goodness of fit. It is defined as:
\[ \text{AIC} = 2k - 2\ln(\hat{L}) \]
where \(k\) is the number of parameters in the model, and \(\hat{L}\) is the maximum likelihood of the model. Lower AIC values indicate a better fit, as they balance the trade-off between model complexity and goodness of fit.
An even better way to assess the quality of the model is to compute the specificity, sensitivity and the misclassification rate, and then the AUC value. To do that, we first need to calculate the predicted values of \(y\). The predicted values \(\hat{y}\) are calculated as:
\[ \text{P(y=1)} = \hat{y} = \frac{1}{1 + \exp(-\hat{\beta}_0 - \hat{\beta}_1x_1 - \dots - \hat{\beta}_kx_k)} \] where \(\hat{y}\) represents the predicted probability of the outcome being \(y = 1\).
After obtaining the predicted probabilities, we apply a threshold (cut-off) to classify predictions into binary outcomes (e.g. classify \(\hat{y} > 0.5\) as \(y = 1\), and \(\hat{y} \leq 0.5\) as \(y = 0\)). It is important to note that the choice of cut-off value (i.e., threshold for classifying \(y = 1\)) impacts the specificity, sensitivity, and misclassification rate. The default threshold of 0.5 may not always be optimal. Using different cut-offs allows for a trade-off between sensitivity and specificity. For example, a lower cut-off might increase sensitivity at the expense of specificity, and vice versa.
Then, we could calculate the sensitivity and specificity of the model. Sensitivity is the proportion of true positives (actual positives correctly predicted as positives), while specificity is the proportion of true negatives (actual negatives correctly predicted as negatives). The misclassification rate is the proportion of incorrect predictions made by the model. They are calculated as:
\[ \text{Sensitivity} = \frac{TP}{TP + FN} \] where \(TP\) is the number of true positives and \(FN\) is the number of false negatives. Higher sensitivity values are better, as they indicate that the model is good at identifying the positive class.
\[ \text{Specificity} = \frac{TN}{TN + FP} \]
where \(TN\) is the number of true negatives and \(FP\) is the number of false positives. Higher specificity values are better, as they indicate that the model is good at identifying the negative class.
\[ \text{Misclassification Rate} = \frac{FP + FN}{TP + TN + FP + FN} \]
Lower misclassification rates are better, as they indicate fewer errors in classification.
Different cutoff can lead to different sensitivity, specificity and misclassification rates.To determine the optimal cut-off for classification, we can calculate the Youden’s Index, which is defined as:
\[ J = \text{Sensitivity} + \text{Specificity} - 1 \] Alternatively, we can use the ROC curve to visualize the trade-off between sensitivity and specificity at different cut-off values, and find a cut-off for which the ROC curve has the minimum distance from the upper left corner of the graph (the point at which sensitivity and specificity are both 1). We implement this
To provide more cotext, the ROC curve is then a graphical representation of a classifier’s performance, plotting sensitivity (true positive rate) on the y-axis against 1 - specificity (false positive rate) on the x-axis for different threshold values. The ROC curve helps visualize the trade-off between sensitivity and specificity at various cut-offs.
Based on the ROC curve, we can calculate the Area Under the Curve (AUC) to quantify the model’s performance. Area under ROC Curve is a measure of prediction accuracy of the model (how well a model predicts 1 responses as 1’s and 0 responses as 0’s). The AUC ranges from 0 to 1, with higher values indicating better model performance. An AUC of 0.5 suggests that the model is no better than random guessing, while an AUC of 1 indicates perfect classification. Higher AUCs mean that we can find a cut-off value for which both sensitivity and specificity of the model are relatively high. Typically:
- \(\text{AUC} = 1\): Perfect model
- \(\text{AUC} = 0.9\): Excellent model
- \(\text{AUC} = 0.8\): Good model
- \(\text{AUC} = 0.7\): Fair model
- \(\text{AUC} = 0.5\): Poor model (random classifier)
- \(\text{AUC} < 0.5\): Failing model (worse than random)
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. Below, we outline which assumptions of OLS regression also hold for logistic regression and which do not.
First, just like in OLS regression, logistic regression assumes that observations are independent of each other. This means that each data point (or case) should not influence another.
Second, similar to OLS regression, logistic regression assumes that there is no perfect multicollinearity among the independent variables. Perfect multicollinearity means that one predictor is an exact linear combination of others, which can lead to instability in the model coefficients.
Third, OLS regression assumes homoscedasticity, meaning the residuals (errors) have constant variance across all levels of the independent variables. In logistic regression, this assumption does not apply because the error terms are not normally distributed; instead, they follow a binomial distribution.
Fourth, 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 for the response variable.
And lastly, logistic regression does not assume a linear relationship between the dependent variable and each independent variable, unlike OLS regression. Instead, dependent variable must be binary. Logistic regression assumes that the log odds of the dependent variable being 1 is linearly related to each independent variable. This is a unique assumption of logistic regression that does not apply to OLS regression.
As a special note, logistic regression requires a sample size of at least 50. This is because logistic regression uses Maximum Likelihood Estimation (MLE) to estimate regression coefficients, which requires a larger sample size than OLS regression.
Exploratory Analysis Ideas
Before running logistic regression, we also need to conduct exploratory analyses to understand the relationships between the dependent variable and the predictors.
For binary predictors, we run cross-tabulations with the dependent variable, which in this case is \(\text{DRINKING_D}\). This provides us with counts for each combination of values for \(\text{DRINKING_D}\) and a binary predictor, allowing us to observe potential associations. The Chi-Square test is the appropriate statistical test for assessing associations between two categorical variables. The null hypothesis \(H_0\) for the Chi-Square test is that:
\[ H_0: \text{There is no association between the binary predictor and DRINKING_D.} \]
The alternative hypothesis is that:
\[ H_a: \text{There is an association between the binary predictor and DRINKING\_D.} \]
For continuous predictors, we calculate the means of the continuous predictors for both values of the dependent variable. We compare the mean values of \(\text{PCTBACHMOR}\) and \(\text{MEDHHINC}\) for crashes involving alcohol versus those that didn’t. The independent samples t-test is used to test whether there is a significant difference in mean values of a continuous predictor between two groups of a binary outcome. The null hypothesis \(H_0\) for the t-test is that:
\[ H_0: \text{There is no difference in the means of the continuous predictor for crashes that involved alcohol and those that did not.} \]
The alternative hypothesis is that:
\[ H_a: \text{There is a significant difference in the means of the continuous predictor for crashes that involved alcohol and those that did not.} \]
Results
Exploratory Analysis
We begin our analysis by examining the distribution of the dependent variable, \(\text{DRINKING_D}\). The tabulation shows that the majority of crashes do not involve alcohol impairment. The tabulation of the dependent variable, \(\text{DRINKING_D}\), shows that the majority of crashes do not involve alcohol impairment. Specifically, 94.3% of crashes (proportion = 0.943) are categorized as not involving drunk driving, while only 5.7% of crashes (proportion = 0.057) involve a drunk driver. This indicates that alcohol-related crashes constitute a relatively small proportion of the overall crash data. The skewed distribution of the dependent variable suggests that drunk driving incidents are relatively rare compared to non-alcohol-related crashes, which is an important consideration for model performance and evaluation, as it may impact the model’s ability to accurately predict alcohol involvement.
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 |
The table below provides a cross-tabulation of the dependent variable with each of the binary predictors. It shows that among crashes without alcohol involvement, 2.89% resulted in fatality or major injury. Yet, the percentage is higher at 7.57% when alcohol is involved. While 1.50% of non-alcohol-related crashes involved an overturned vehicle, 4.43% of alcohol-related crashes did. Speeding is involved in 3.08% of non-alcohol-related crashes and 10.46% of alcohol-related crashes. Similarly, aggressive driving is more common in alcohol-related crashes (36.86%) than in non-alcohol-related crashes (45.31%). As for driving age, 1.65% of non-alcohol-related crashes involve very young drivers, while this drops to 0.48% in alcohol-related crashes. 10.37% of non-alcohol-related crashes involve older drivers, compared to 4.79% of alcohol-related crashes. However, cell phone use is similar across both alcohol-involved and non-alcohol-involved crashes.
extract_crosstab_values <- function(var1, var2) {
crosstab <- CrossTable(var1, var2, prop.r = FALSE, prop.chisq = FALSE, chisq = FALSE, prop.t = FALSE)
value_1 <- crosstab$t["1", "1"] # Count where DRINKING_D = 1, other variable = 1
value_2 <- (crosstab$t["1", "1"] / sum(crosstab$t["1",])) * 100 # Pct where DRINKING_D = 1, other variable = 1
value_3 <- crosstab$t["0", "1"] # Count where DRINKING_D = 0, other variable = 1
value_4 <- (crosstab$t["0", "1"] / sum(crosstab$t["0",])) * 100 # Pct where DRINKING_D = 0, other variable = 1
value_5 <- crosstab$t["1", "1"] + crosstab$t["0", "1"] # Total count where other variable = 1
return(c(Value1 = value_1, Value2 = value_2, Value3 = value_3, Value4 = value_4, Value5 = value_5))
}
tabulation <- data.frame(
Variable = character(),
Value1 = numeric(),
Value2 = numeric(),
Value3 = numeric(),
Value4 = numeric(),
Value5 = numeric()
)
descriptions <- c(
FATAL_OR_M = "Crash resulted in fatality or major injury",
OVERTURNED = "Crash involved an overturned vehicle",
CELL_PHONE = "Driver was using cell phone",
SPEEDING = "Crash involved speeding car",
AGGRESSIVE = "Crash involved aggressive driving",
DRIVER1617 = "Crash involved at least one driver 16 or 17 years old",
DRIVER65PLUS = "Crash involved at least one driver over 65 years old"
)
variables <- c("FATAL_OR_M", "OVERTURNED", "CELL_PHONE", "SPEEDING", "AGGRESSIVE", "DRIVER1617", "DRIVER65PLUS")
for (var in variables) {
values <- extract_crosstab_values(crashes$DRINKING_D, crashes[[var]])
tabulation <- rbind(
tabulation,
data.frame(Variable = paste0(var, ": ", descriptions[[var]]), t(values))
)
}
colnames(tabulation) <- c("", "N", "%", "N ", "% ", "N ")tabulation %>%
kable("html", escape = FALSE, align = "lccccc") %>%
kable_styling(full_width = FALSE, position = "center", fixed_thead = TRUE) %>%
column_spec(1, width = "30em") %>%
column_spec(2, width = "5em") %>%
column_spec(3, width = "5em") %>%
column_spec(4, width = "5em") %>%
column_spec(5, width = "5em") %>%
column_spec(6, width = "5em") %>%
add_header_above(c(" " = 1, "DRINKING_D = 1" = 2, "DRINKING_D = 0" = 2, " " = 1)) %>%
add_header_above(c(" " = 1, "Alcohol Involved" = 2, "No Alcohol" = 2, "Total" = 1))| N | % | N | % | N | |
|---|---|---|---|---|---|
| FATAL_OR_M: Crash resulted in fatality or major injury | 188 | 7.565 | 1181 | 2.89 | 1369 |
| OVERTURNED: Crash involved an overturned vehicle | 110 | 4.427 | 612 | 1.50 | 722 |
| CELL_PHONE: Driver was using cell phone | 28 | 1.127 | 426 | 1.04 | 454 |
| SPEEDING: Crash involved speeding car | 260 | 10.463 | 1261 | 3.08 | 1521 |
| AGGRESSIVE: Crash involved aggressive driving | 916 | 36.861 | 18522 | 45.31 | 19438 |
| DRIVER1617: Crash involved at least one driver 16 or 17 years old | 12 | 0.483 | 674 | 1.65 | 686 |
| DRIVER65PLUS: Crash involved at least one driver over 65 years old | 119 | 4.789 | 4237 | 10.37 | 4356 |
The results of the Chi-Square test indicate significant associations between most of the binary predictors and alcohol involvement in crashes, with the exception of \(\text{CELL_PHONE}\). For example, \(\text{FATAL_OR_M}\) showed a highly significant association, with a Chi-Square statistic of 167.56 and a highly significant \(p < 0.0001\), suggesting that severe crashes are more likely to involve alcohol. Likewise, \(\text{OVERTURNED}\) also had a significant association (Chi-Square = 122.79, df =1) with \(p < 0.0001\), indicating that overturned vehicles are more frequently associated with alcohol impairment. The remaining predictors, \(\text{SPEEDING}\), \(\text{AGGRESSIVE}\), \(\text{DRIVER1617}\), and \(\text{DRIVER65PLUS}\) also showed significant associations with alcohol involvement, with Chi-Square statistics ranging from 20.45 to 376.78 and \(p < 0.0001\). These results suggest that crashes involving fatalities, overturned vehicles, speeding, aggressive driving, and certain driver age groups are more likely to involve alcohol impairment. Therefore, we could reject the null hypothesis for these predictors and conclude that they are significantly associated with alcohol involvement in crashes.
However, the predictor \(\text{CELL_PHONE}\) did not show a significant association with alcohol involvement, with a Chi-Square statistic of 0.16 and \(p > 0.5\). This suggests that cell phone use is not strongly associated with alcohol-related crashes in this dataset. We can retain the null hypothesis for this predictor and conclude that there is no significant association between cell phone use and alcohol involvement in crashes.
tabulation$Chi_squared <- NA
tabulation$p_value <- NA
for (i in seq_along(variables)) {
var <- variables[i]
chi_test <- chisq.test(table(crashes$DRINKING_D, crashes[[var]]), correct = FALSE)
tabulation$Chi_squared[i] <- round(chi_test$statistic, 2)
tabulation$p_value[i] <- round(chi_test$p.value, 5)
}
tabulation %>%
kable("html", escape = FALSE, align = "lcccccc") %>%
kable_styling(full_width = FALSE, position = "center", fixed_thead = TRUE) %>%
column_spec(1, width = "50em") %>%
column_spec(2, width = "5em") %>%
column_spec(3, width = "5em") %>%
column_spec(4, width = "5em") %>%
column_spec(5, width = "5em") %>%
column_spec(6, width = "5em") %>%
column_spec(7, width = "5em") %>%
column_spec(8, width = "5em") %>%
add_header_above(c(" " = 1, "DRINKING_D = 0" = 2, "DRINKING_D = 1" = 2, " " = 1, " " = 2)) %>%
add_header_above(c(" " = 1, "No Alcohol" = 2, "Alcohol Involved" = 2, "Total" = 1, "χ2 Test" = 2))| N | % | N | % | N | Chi_squared | p_value | |
|---|---|---|---|---|---|---|---|
| FATAL_OR_M: Crash resulted in fatality or major injury | 188 | 7.565 | 1181 | 2.89 | 1369 | 167.56 | 0.000 |
| OVERTURNED: Crash involved an overturned vehicle | 110 | 4.427 | 612 | 1.50 | 722 | 122.79 | 0.000 |
| CELL_PHONE: Driver was using cell phone | 28 | 1.127 | 426 | 1.04 | 454 | 0.16 | 0.687 |
| SPEEDING: Crash involved speeding car | 260 | 10.463 | 1261 | 3.08 | 1521 | 376.78 | 0.000 |
| AGGRESSIVE: Crash involved aggressive driving | 916 | 36.861 | 18522 | 45.31 | 19438 | 67.60 | 0.000 |
| DRIVER1617: Crash involved at least one driver 16 or 17 years old | 12 | 0.483 | 674 | 1.65 | 686 | 20.45 | 0.000 |
| DRIVER65PLUS: Crash involved at least one driver over 65 years old | 119 | 4.789 | 4237 | 10.37 | 4356 | 80.60 | 0.000 |
We also looked at the mean and standard deviation of the continuous predictors, \(\text{PCTBACHMOR}\) and \(\text{MEDHHINC}\), for crashes involving and not involving alcohol. The mean percentage of people with a bachelor’s degree or more is identical for both non-alcohol and alcohol-involved crashes, at 16.6%. The standard deviation is slightly higher in the alcohol-involved group (18.7 vs. 18.2). This suggests that educational attainment, as measured by the percentage of individuals with a bachelor’s degree or more, does not vary significantly between areas where crashes involve alcohol and those where they do not.The median household income is slightly higher in alcohol-involved crashes with a mean of $31,998.8, compared to $31,483.1 in non-alcohol-related crashes. The standard deviation is also higher for the alcohol-involved group (17,810.5 vs. 16,930.1). This also means that income level may not be a strong differentiator for alcohol involvement in crashes.
mean_pctbachmor <- tapply(crashes$PCTBACHMOR, crashes$DRINKING_D, mean)
sd_pctbachmor <- tapply(crashes$PCTBACHMOR, crashes$DRINKING_D, sd)
mean_medhhinc <- tapply(crashes$MEDHHINC, crashes$DRINKING_D, mean)
sd_medhhinc <- tapply(crashes$MEDHHINC, crashes$DRINKING_D, sd)
tabulation_cont <- data.frame(
Variable = c("PCTBACHMOR: % with bachelor’s degree or more", "MEDHHINC: Median household income"),
Value1 = c(mean_pctbachmor["0"], mean_medhhinc["0"]), # Mean for DRINKING_D = 0
Value2 = c(sd_pctbachmor["0"], sd_medhhinc["0"]), # SD for DRINKING_D = 0
Value3 = c(mean_pctbachmor["1"], mean_medhhinc["1"]), # Mean for DRINKING_D = 1
Value4 = c(sd_pctbachmor["1"], sd_medhhinc["1"]) # SD for DRINKING_D = 1
)
colnames(tabulation_cont) <- c("", "Mean", "SD", "Mean ", "SD ")
tabulation_cont %>%
kable("html", escape = FALSE, align = "lcccccc") %>%
kable_styling(full_width = FALSE, position = "center", fixed_thead = TRUE) %>%
add_header_above(c(" " = 1, "DRINKING_D = 0" = 2, "DRINKING_D = 1" = 2)) %>%
add_header_above(c(" " = 1, "No Alcohol" = 2, "Alcohol Involved" = 2))| Mean | SD | Mean | SD | |
|---|---|---|---|---|
| PCTBACHMOR: % with bachelor’s degree or more | 16.6 | 18.2 | 16.6 | 18.7 |
| MEDHHINC: Median household income | 31483.1 | 16930.1 | 31998.8 | 17810.5 |
The independent samples t-test results for the continuous predictors provide more insight into whether these socio-economic factors are significantly associated with alcohol involvement in crashes.
For \(\text{PCTBACHMOR}\), the t-test yields a p-value of 0.914, which is far greater than the conventional significance level of 0.05. This indicates that there is no statistically significant difference in the mean educational attainment between crashes involving alcohol and those that do not. Therefore, we fail to reject the null hypothesis, suggesting that educational attainment does not have a significant association with alcohol involvement in crashes.
Similarly, for \(\text{MEDHHINC}\), the t-test produces a p-value of 0.160, which also exceeds the 0.05 threshold. Although there is a small difference in mean income between the two groups, this difference is not statistically significant. Consequently, we again fail to reject the null hypothesis, implying that median household income is not significantly associated with alcohol involvement in crashes.
As such, the t-test results demonstrate that neither \(\text{PCTBACHMOR}\) nor \(\text{MEDHHINC}\) shows a significant association with alcohol involvement in crashes. Thus, we cannot conclude that socio-economic factors, as measured by educational attainment and income levels, have a meaningful impact on the likelihood of alcohol involvement in traffic incidents. This lack of significance suggests that other variables, such as driver behavior or crash-specific characteristics, may be more relevant in understanding alcohol-related crashes.
pval_pctbachmor <- t.test(crashes$PCTBACHMOR ~ crashes$DRINKING_D)$p.value
pval_medhhinc <- t.test(crashes$MEDHHINC ~ crashes$DRINKING_D)$p.value
tabulation_cont$`t-test p-value` <- c(pval_pctbachmor, pval_medhhinc)
tabulation_cont %>%
kable("html", escape = FALSE, align = "lcccccc") %>%
kable_styling(full_width = FALSE, position = "center", fixed_thead = TRUE) %>%
add_header_above(c(" " = 1, "DRINKING_D = 0" = 2, "DRINKING_D = 1" = 2, " " = 1)) %>%
add_header_above(c(" " = 1, "No Alcohol" = 2, "Alcohol Involved" = 2, " " = 1))| Mean | SD | Mean | SD | t-test p-value | |
|---|---|---|---|---|---|
| PCTBACHMOR: % with bachelor’s degree or more | 16.6 | 18.2 | 16.6 | 18.7 | 0.914 |
| MEDHHINC: Median household income | 31483.1 | 16930.1 | 31998.8 | 17810.5 | 0.160 |
Regression Assumptions
Based on the correlation matrix and prior information, most logistic regression assumptions for analyzing drunk driving crashes are met. The binary outcome variable, \(\text{DRINKING_D}\), fulfills the requirement for a binary dependent variable, and the independence of observations is largely satisfied, as each crash is considered a separate event. The assumption of no perfect multicollinearity is also met, as the correlation matrix shows low pairwise correlations between predictor variables, with the highest being 0.48 between % with Bachelor’s Degree and Median Household Income, and a mild correlation (0.21) between Crash Involved Speeding Car and Crash Involved Aggressive Driving. These correlations are not high enough to cause instability in the model. Additionally, the large sample size of over 43,000 observations ensures stable estimates, meeting 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-tests suggest weak or non-linear relationships between % with Bachelor’s Degree and Median Household Income and the outcome variable.
Using Pearson correlations to measure associations between binary predictors can have limitations, as Pearson’s correlation coefficient is most appropriate for continuous, normally distributed data. Binary variables, which only take on two values (e.g., 0 and 1), can lead to less meaningful or interpretable correlation values, as the measure does not fully capture the potential association between two binary variables. Additionally, Pearson correlations may underestimate or fail to reflect non-linear relationships that could exist among binary predictors.
In terms of multicollinearity, which we define as the presence of high correlation between two or more predictor variables that could lead to instability in the model’s coefficient estimates, there is no evidence of multicollinearity in this dataset. The correlation matrix shows low pairwise correlations between the predictors, with the highest correlation being 0.48 between % with Bachelor’s Degree and Median Household Income. This value falls below commonly accepted thresholds for multicollinearity (typically 0.7 or higher), suggesting that the model is unlikely to suffer from instability due to collinear predictors. Therefore, while the Pearson correlation may not be the most precise measure for binary predictors, it does not indicate any problematic levels of multicollinearity in this analysis.
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("#283d3b", "white", "#c44536")) +
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))
Regression Analysis
Our first logistic regression model includes all predictors and the results of it is presented below.
FATAL_OR_M: The presence of a fatality or major injury significantly increases the odds of a drinking-related crash, with an odds ratio of 2.257 (95% CI: 1.910 to 2.653) and a highly significant \(p < 0.0001\), suggesting that these types of incidents more than double the odds of it there being a drunk driver.
OVERTURNED: Overturned vehicles significantly increase the likelihood of it being a drinking-related crash, with an odds ratio of 2.532 (95% CI: 2.035 to 3.122) and a highly significant \(p < 0.0001\). This result implies that overturned vehicles are strongly associated with a drinking-driver incurred crashes.
CELL_PHONE: The use of a cellphone during a crash shows no statistically significant effect on the odds of a drinking-related crash (\(p = 0.881\)), with an odds ratio close to 1 (OR = 1.030, 95% CI: 0.684 to 1.488), indicating no significantly meaningful association.
SPEEDING: Speeding is a highly significant predictor, with an odds ratio of 4.660 (95% CI: 3.974 to 5.450), and \(p < 0.0001\), indicating that crashes involving speeding are more than four times as likely to involve alcohol.
AGGRESSIVE: Aggressive driving is associated with lower odds of a drinking-related crash, with an odds ratio of 0.551 (95% CI: 0.501 to 0.604), and \(p < 0.0001\), suggesting that aggressive driving incidents are less likely to involve alcohol.
DRIVER1617: Incidents involving drivers aged 16-17 are significantly less likely to be drinking-related (OR = 0.278, 95% CI: 0.148 to 0.471), with a significant \(p < 0.0001\), possibly reflecting the influence of drinking laws for younger drivers.
DRIVER65PLUS: Incidents involving drivers aged 65 and above also show significantly lower odds of being drinking-related, with an odds ratio of 0.461 (95% CI: 0.380 to 0.553), and \(p < 0.0001\), indicating that older drivers are less likely to be involved in alcohol-related crashes.
PCTBACHMOR: The percentage of people with a bachelor’s degree or higher in the area does not significantly impact the odds of a drinking-related crash (\(p = 0.775\)) that occur there, with an odds ratio near 1 (OR = 1.000, 95% CI: 0.997 to 1.002).
MEDHHINC: Median household income has a very slight, statistically significant positive effect on the likelihood of a drinking-related crash there, with \(p = 0.036\), but the effect size is negligible, as shown by the confidence interval limits around 1 (OR = 1.000).
logit <- glm(DRINKING_D ~ FATAL_OR_M + OVERTURNED + CELL_PHONE + SPEEDING + AGGRESSIVE + DRIVER1617 + DRIVER65PLUS + PCTBACHMOR + MEDHHINC, data = crashes, family = binomial)
summary(logit)##
## Call:
## glm(formula = DRINKING_D ~ FATAL_OR_M + OVERTURNED + CELL_PHONE +
## SPEEDING + AGGRESSIVE + DRIVER1617 + DRIVER65PLUS + PCTBACHMOR +
## MEDHHINC, family = binomial, data = crashes)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.73250662 0.04587566 -59.56 < 0.0000000000000002 ***
## FATAL_OR_M 0.81401380 0.08380692 9.71 < 0.0000000000000002 ***
## OVERTURNED 0.92892138 0.10916632 8.51 < 0.0000000000000002 ***
## CELL_PHONE 0.02955008 0.19777782 0.15 0.881
## SPEEDING 1.53897567 0.08054589 19.11 < 0.0000000000000002 ***
## AGGRESSIVE -0.59691595 0.04777924 -12.49 < 0.0000000000000002 ***
## DRIVER1617 -1.28029596 0.29314717 -4.37 0.00001257244712793 ***
## DRIVER65PLUS -0.77466464 0.09585832 -8.08 0.00000000000000064 ***
## PCTBACHMOR -0.00037063 0.00129639 -0.29 0.775
## MEDHHINC 0.00000280 0.00000134 2.09 0.036 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 19036 on 43363 degrees of freedom
## Residual deviance: 18340 on 43354 degrees of freedom
## AIC: 18360
##
## Number of Fisher Scoring iterations: 6As a quick summary, predictor that are highly significant \(p < 0.001\) are FATAL_OR_M, OVERTURNED, SPEEDING, AGGRESSIVE, DRIVER1617, and DRIVER65PLUS. Among which those that increase the odds of an incident that includes drunk driving are FATAL_OR_M, OVERTURNED, and SPEEDING. On the other hand, predictors that are not significant are CELL_PHONE and PCTBACHMOR. The predictor MEDHHINC is marginally significant (\(p = 0.036\)).
logitoutput <- summary(logit)
logitcoeff <- logitoutput$coefficients
or_ci <- exp(cbind(OR = coef(logit), confint(logit)))
logitcoeff[, "Pr(>|z|)"] <- round(logitcoeff[, "Pr(>|z|)"], 4)
final_output <- cbind(logitcoeff, or_ci)
final_output %>%
kable("html", escape = FALSE, align = "cccccc") %>%
kable_styling(full_width = FALSE, position = "center", fixed_thead = TRUE)| 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 |
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 lowest and highest misclassification rates highlight the effectiveness of each cut-off.
The probability cut-off of 0.50 yields the lowest misclassification rate of 0.057. At this threshold, specificity is maximized at 1.000, but sensitivity is very low (0.002), meaning that the model captures very few true positives, sacrificing sensitivity for high specificity.
The highest misclassification rate occurs at the cut-off of 0.02 with a rate of 0.889. At this threshold, sensitivity is high (0.984), capturing nearly all true positives, but specificity is very low (0.058), resulting in a high rate of false positives.
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,
Specificity = NA,
MisclassificationRate = NA
)
calculate_metrics <- function(cutoff, labels, predictions) {
predicted_labels <- ifelse(predictions >= cutoff, 1, 0)
tp <- sum(predicted_labels == 1 & labels == 1) # True positives
tn <- sum(predicted_labels == 0 & labels == 0) # True negatives
fp <- sum(predicted_labels == 1 & labels == 0) # False positives
fn <- sum(predicted_labels == 0 & labels == 1) # False negatives
# Sensitivity: True positive rate
sensitivity <- tp / (tp + fn)
# Specificity: True negative rate
specificity <- tn / (tn + fp)
# Misclassification rate
misclassification_rate <- (fp + fn) / length(labels)
return(c(sensitivity, specificity, misclassification_rate))
}
for (i in seq_along(cutoffs)) {
metrics <- calculate_metrics(cutoffs[i], roc$labels, roc$predictions)
results$Sensitivity[i] <- metrics[1]
results$Specificity[i] <- metrics[2]
results$MisclassificationRate[i] <- metrics[3]
}
results %>%
kable("html", escape = FALSE, align = "cccc") %>%
kable_styling(full_width = FALSE, position = "center", fixed_thead = TRUE)| 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 |
We plotted the ROC curve to show the tradeoff between sensitivity (True Positive Rate) and specificity (1 - False Positive Rate) across various probability cut-offs for the logistic regression model. An ideal model would reach the top-left corner of the ROC space (sensitivity = 1, specificity = 1), which would mean perfect discrimination between the two classes. Our ROC curve, however, lies beneath this 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)
)
Using the function opt.cut, we calculated the optimal
cut-off by minimizing the Euclidean distance from the ROC curve to the
point (0, 1) — the top-left corner of the ROC space, which represents
perfect sensitivity and specificity. This cut-off of 0.064 provides a
balanced approach, achieving a sensitivity of 0.661 and specificity of
0.545.
In the above section, the optimal cut-off selected based on the minimum misclassification rate was 0.50, which yielded a misclassification rate of 0.057. This threshold focused on reducing the overall number of misclassified observations, prioritizing high specificity (1.000) but resulting in low sensitivity (0.002). In contrast, the cut-off of 0.064 identified through the ROC analysis aims to balance both sensitivity and specificity, achieving a moderate level for each metric rather than prioritizing one over the other. This cut-off might be more suitable for applications where both true positive and true negative rates are equally important.
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(full_width = FALSE, position = "center", fixed_thead = TRUE)| sensitivity | specificity | cutoff |
|---|---|---|
| 0.661 | 0.545 | 0.064 |
The Area Under the ROC Curve (AUC) for our model is 0.640. An AUC of 0.640 suggests that the model has poor discrimination (generally, a model with an AUC between 0.6 and 0.7 is considered to have poor discrimination). Specifically, the AUC value means that if we randomly select one positive case and one negative case, the model has a 64% chance of correctly ranking the positive case higher in terms of predicted probability than the negative case. In practical terms, 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.640The second logistic regression model we ran includes only the binary predictors. In the first model, \(\text{PCTBACHMOR}\) and \(\text{MEDHHINC}\) were not statistically significant predictors of \(\text{DRINKING_D}\), as their \(p < 0.05\). However, in the reduced model, the exclusion of these two predictors did not lead to any new predictors becoming significant or any previously significant predictors becoming insignificant.
logit2 <- glm(DRINKING_D ~ FATAL_OR_M + OVERTURNED + CELL_PHONE + SPEEDING + AGGRESSIVE + DRIVER1617 + DRIVER65PLUS, data = crashes, family = binomial)
summary(logit2)##
## Call:
## glm(formula = DRINKING_D ~ FATAL_OR_M + OVERTURNED + CELL_PHONE +
## SPEEDING + AGGRESSIVE + DRIVER1617 + DRIVER65PLUS, family = binomial,
## data = crashes)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.6519 0.0275 -96.32 < 0.0000000000000002 ***
## FATAL_OR_M 0.8093 0.0838 9.66 < 0.0000000000000002 ***
## OVERTURNED 0.9398 0.1090 8.62 < 0.0000000000000002 ***
## CELL_PHONE 0.0311 0.1978 0.16 0.88
## SPEEDING 1.5403 0.0805 19.13 < 0.0000000000000002 ***
## AGGRESSIVE -0.5936 0.0477 -12.43 < 0.0000000000000002 ***
## DRIVER1617 -1.2716 0.2931 -4.34 0.0000143637414327 ***
## DRIVER65PLUS -0.7665 0.0958 -8.00 0.0000000000000012 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 19036 on 43363 degrees of freedom
## Residual deviance: 18344 on 43356 degrees of freedom
## AIC: 18360
##
## Number of Fisher Scoring iterations: 6The odds ratios (OR) for the binary predictors in both models are very similar, suggesting that removing \(\text{PCTBACHMOR}\) and \(\text{MEDHHINC}\) had little impact on the estimated effects of the other variables.
logitoutput2 <- summary(logit2)
logitcoeff2 <- logitoutput2$coefficients
or_ci2 <- exp(cbind(OR = coef(logit2), confint(logit2)))
logitcoeff2[, "Pr(>|z|)"] <- round(logitcoeff2[, "Pr(>|z|)"], 4)
final_output2 <- cbind(logitcoeff2, or_ci2)
final_output2 %>%
kable("html", escape = FALSE, align = "cccccc") %>%
kable_styling(full_width = FALSE, position = "center", fixed_thead = TRUE)| 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 . AIC is a measure of model quality that balances goodness of fit with model complexity (penalizing additional parameters). Since the AIC is the same in both models, neither model is superior in terms of fit. Given that the first model includes two additional predictors that are not statistically significant and does not improve the AIC, the reduced model (binary predictors only) is preferable for its simplicity.
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(full_width = FALSE, position = "center", fixed_thead = TRUE)| Degree of Freedom | AIC | |
|---|---|---|
| Model 1: All Predictors | 10 | 18360 |
| Model 2: Binary Predictors Only | 8 | 18360 |
Discussion
In this study, a logistic regression analysis was conducted to examine factors associated with alcohol involvement in vehicular crashes in Philadelphia. The analysis included a variety of predictors, encompassing crash characteristics (such as whether the crash resulted in a fatality or major injury, involved an overturned vehicle, or involved cell phone use), driver behaviors (including speeding, aggressive driving, and the age of drivers involved), and socio-economic variables from the crash locations (such as the percentage of individuals with a bachelor’s degree or higher and median household income). The primary objective was to identify significant predictors of alcohol-related crashes and gain insight into the relationships between these factors and the likelihood of drunk driving incidents.
The results indicate that several crash characteristics and driver behaviors are significant predictors of alcohol involvement. Specifically, \(\text{FATAL_OR_M}\) (whether the crash resulted in a fatality or major injury), \(\text{OVERTURNED}\) (whether the crash involved an overturned vehicle), and ( ) are all positively and significantly associated with alcohol involvement. For instance, crashes involving speeding were found to be more than four times as likely to involve a drunk driver, highlighting the strong connection between risky driving behaviors and alcohol impairment. In contrast, ( ) (use of a cell phone during the crash) and \(\text{PCTBACHMOR}\) (percentage of individuals with a bachelor’s degree or higher) did not show a significant association with alcohol involvement, suggesting that these factors do not substantially influence the likelihood of a crash involving alcohol. \(\text{MEDHHINC}\) (median household income) demonstrated a marginally significant effect, but the impact was negligible, indicating that income levels do not play a meaningful role in alcohol-related crash risk.
The findings are somewhat surprising. While the significance of variables like \(\text{SPEEDING}\) and \(\text{FATAL_OR_M}\) matches our expectations, the lack of a significant association for \(\text{CELL_PHONE}\) use is unexpected. It was anticipated that cell phone use, an already high-risk behavior, would have an amplified effect when combined with alcohol impairment, but the data did not support this assumption. Similarly, the limited impact of socio-economic variables was unforeseen. Although it was hypothesized that areas with lower educational attainment or income levels might experience higher rates of alcohol-related crashes due to various socio-economic stressors, the analysis did not find meaningful associations. Nevertheless, the significant predictors, such as speeding and crash severity, align with expectations, reinforcing the understanding that alcohol impairment exacerbates the risk of severe and reckless driving behaviors.
The application of logistic regression is appropriate in this context, given its suitability for modeling a binary dependent variable. However, considering the rarity of the event of interest—only 5.7% of crashes involved alcohol—there may be limitations to this approach. Methods for modeling rare events, as proposed by Paul Allison, such as penalized likelihood approaches, could be more effective. These methods may yield more reliable and efficient estimates when dealing with infrequent outcomes, addressing potential biases and enhancing the model’s performance.
This analysis has several limitations. The relatively small proportion of alcohol-involved crashes may affect the robustness and stability of the model estimates. Employing rare event modeling techniques could improve the reliability of the findings. Additionally, unobserved confounding variables, such as time of day, road conditions, or proximity to alcohol-serving establishments, may influence crash outcomes but were not included in the model. The binary nature of many predictors may also oversimplify complex driving behaviors or socio-economic conditions, limiting the model’s ability to capture nuances. Finally, while socio-economic factors did not emerge as significant predictors, this does not rule out their potential importance in a broader context or when analyzed with more granular data. Future research could address these limitations by incorporating additional contextual variables, utilizing alternative modeling methods, and expanding the analysis to other geographic areas to validate and extend these findings.