--- title: "Traffic Accident Analysis 2025" output: html_document --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE) ``` # Load Libraries ```{r libraries} library(ggplot2) library(dplyr) library(plotly) library(caret) library(rpart) library(rpart.plot) library(randomForest) library(gbm) library(themis) library(pROC) library(glmnet) library(fastshap) library(shapviz) library(tidyverse) library(MASS) library(e1071) ``` # 1. PREPROCESSING AND VISUALIZATION ### 1 Quick study and inspection of the dataset ```{r inspection} readLines("2025_Accidentalidad.csv", 2) db = read.csv(file ="2025_Accidentalidad.csv", header = TRUE, sep = ";", dec=".", stringsAsFactors = TRUE) summary(db) str(db) glimpse(db) ``` ## 1.1 Remove useless rows, NAs and duplicated rows **Remove duplicates** ```{r duplicates} sum(duplicated(db)) db <- db[!duplicated(db), ] ``` **Num and num_expediente, are the same informations, but we do not need them, the value of the row is enough** ```{r remove_ids} db$num_expediente = NULL db$numero = NULL ``` **See which variables have NAs** ```{r check_nas} names(db)[colSums(is.na(db)) > 0] ``` **The geographic coordinates were removed as they were redundant with the district variable, contained a high number of missing values, and were not needed for the analysis** ```{r remove_coords} db$coordenada_x_utm = NULL db$coordenada_y_utm = NULL ``` **In the lesividad variable, blank values indicate cases with no medical assistance. In cod_lesividad, these cases appear as NA, but according to the documentation, the corresponding code for "Sin asistencia sanitaria" is 14.** ```{r clean_lesividad} db$cod_lesividad[is.na(db$cod_lesividad)] = "14" db$lesividad[db$lesividad == ""] = "Sin asistencia sanitaria" ``` **In the variable positiva_droga, the NAs are equal to 0, i.e. negative** ```{r clean_drugs} db$positiva_droga[is.na(db$positiva_droga)] = 0 ``` **In this dataset, missing information is often represented by blank values rather than NA; given its large size, observations with missing values in key variables were removed to ensure consistency without significant information loss** ```{r clean_blanks} colSums(db == "", na.rm = TRUE) vars_blank <- c("estado_meteorológico", "tipo_vehiculo", "positiva_alcohol") db <- db[rowSums(db[vars_blank] == "") == 0, ] rm(vars_blank) ``` ## 1.2 Convert the variables to a better type ```{r conversions} db$fecha = as.Date(db$fecha, format = "%d/%m/%Y") db$cod_distrito = as.factor(db$cod_distrito) db$positiva_alcohol = factor(db$positiva_alcohol, levels=c("N","S"), labels=c("No","Yes")) db$cod_lesividad = as.factor(db$cod_lesividad) db$positiva_droga = factor(db$positiva_droga, levels = c(0, 1), labels = c("No", "Yes")) ``` ### 2. Outliers No outlier analysis was performed, as the dataset does not contain meaningful continuous numerical variables ### 3. Feature engineering In this point we will create new useful variables from the ones that already exist. The first ones are based on the date. Creating month and day ```{r fe_date} db$day = as.factor(format(db$fecha, "%d")) db$month = factor(months(db$fecha)) meses_es <- c("enero","febrero","marzo","abril","mayo","junio", "julio","agosto","septiembre","octubre","noviembre","diciembre") db$month <- factor(months(db$fecha), levels = meses_es) rm(meses_es) ``` Now we will create time slots, separating by Morning, Afternoon, Evening and Night. To do this we start creating an auxiliary variable only with the hour of the accidents. Then we create the final variable using the auxiliary one to separate the accidents by the labels we defined before. The last line of code is to reorder the labels following the logical daily sequence with the morning as the first one ```{r fe_timeslot} hora_num <- as.integer(sub(":.*", "", as.character(db$hora))) db$time_slot <- cut(hora_num, breaks = c(-1, 5, 11, 17, 23), labels = c("Night", "Morning", "Afternoon", "Evening")) rm(hora_num) db$time_slot <- factor(db$time_slot, levels = c("Morning","Afternoon","Evening","Night")) ``` Injury severity was grouped into four meaningful categories (Ileso, Leve, Grave, Mortal) to improve interpretability while preserving clinically relevant distinctions ```{r fe_severity} db$lesividad_cat = factor(ifelse(db$cod_lesividad=="4","Mortal", ifelse(db$cod_lesividad=="3","Grave", ifelse(db$cod_lesividad %in% c("1","2","5","6","7"),"Leve", ifelse(db$cod_lesividad=="14","Ileso","Desconocido")))), levels=c("Ileso","Leve","Grave","Mortal","Desconocido")) ``` Now we will order the rango_edad in order of ages ```{r fe_age_order} orden_edad <- c("Menor de 5 años", "De 6 a 9 años", "De 10 a 14 años", "De 15 a 17 años", "De 18 a 20 años", "De 21 a 24 años", "De 25 a 29 años", "De 30 a 34 años", "De 35 a 39 años", "De 40 a 44 años", "De 45 a 49 años", "De 50 a 54 años", "De 55 a 59 años", "De 60 a 64 años", "De 65 a 69 años", "De 70 a 74 años", "Más de 74 años", "Desconocido") db$rango_edad <- factor(db$rango_edad, levels = orden_edad) rm(orden_edad) ``` Now we create a new variable that group types of vehicles ```{r fe_vehicle_cat} unique(db$tipo_vehiculo) grupos_vehiculos = list( Dos_Ruedas = c("Ciclomotor", "Motocicleta hasta 125cc", "Bicicleta", "Motocicleta > 125cc", "VMU eléctrico", "Bicicleta EPAC (pedaleo asistido)", "Patinete no eléctrico", "Ciclomotor de dos ruedas L1e-B", "Ciclo de motor L1e-A", "Moto de tres ruedas > 125cc", "Ciclo"), Turismo = c("Turismo", "Todo terreno", "Autocaravana", "Cuadriciclo ligero", "Cuadriciclo no ligero"), Pesados_Carga = c("Furgoneta", "Vehículo articulado", "Camión rígido", "Tractocamión", "Semiremolque", "Caravana"), Transporte_Publico = c("Autobús articulado", "Autobús", "Autobus EMT", "Tren/metro"), Otros_Emergencias = c("Maquinaria de obras", "Otros vehículos con motor", "Sin especificar", "Otros vehículos sin motor", "Ambulancia SAMUR") ) db$tipo_vehiculo_cat = db$tipo_vehiculo levels(db$tipo_vehiculo_cat) = grupos_vehiculos rm(grupos_vehiculos) ``` Now we remove the value 'desconocido' in sex and rango_edad ```{r clean_unknown} db = subset(db, sexo != "Desconocido" & rango_edad != "Desconocido" & lesividad_cat != "Desconocido") ``` # VISUALIZATION **Question 1. How does the number of traffic accidents vary by month?** ```{r plot1} p1 = db %>% ggplot(aes(x = month)) + geom_bar(fill = "darkseagreen", alpha = 0.8) + theme_minimal() + labs(title = "Number of accidents by month", x = "Month", y = "Number of accidents") ggplotly(p1) ``` Accident counts remain fairly stable throughout the year, with noticeable peaks in March and October. A clear drop is observed in August, which is likely related to reduced traffic during the holiday period. November and December show zero observations because the dataset corresponds to the current year and, at the time of data collection, information for those months was not yet available. **Question 2. How does the type of traffic accident vary across different weather conditions?** ```{r plot2} p2 = db %>% filter(estado_meteorológico != "Se desconoce") %>% count(estado_meteorológico, tipo_accidente) %>% ggplot(aes(x = estado_meteorológico, y = tipo_accidente, fill = n)) + geom_tile(color = "white", alpha = 0.9) + scale_fill_viridis_c(trans = "log10") + theme_minimal() + labs(title = "Heatmap of accident types by weather condition (log scale)", x = "Weather condition", y = "Accident type", fill = "Number of accidents (log)") + theme(plot.title = element_text(size = 16, face = "bold"), axis.text.x = element_text(angle = 45, hjust = 1)) ggplotly(p2) ``` A logarithmic scale was used to improve readability, as accident counts differ greatly across categories. Most accidents occur under clear weather conditions, which is consistent with DGT reports: when conditions seem safe, drivers tend to relax, increasing accident risk. In contrast, very few accidents occur during snow or hail, likely because drivers are more cautious. Under rain, especially intense rain, accidents are fewer than in clear weather, suggesting that increased attention reduces risk, even though these conditions may seem more dangerous. **Question 3. How does injury severity vary by type of person involved?** ```{r plot3} p3 = db %>% ggplot(aes(x = tipo_persona, fill = lesividad_cat)) + geom_bar(position = "fill", alpha = 0.85) + theme_minimal() + labs(title = "Injury severity distribution by type of person", x = "Type of person", y = "Proportion", fill = "Injury severity") + theme(plot.title = element_text(size = 16, face = "bold")) ggplotly(p3) ``` Injury severity clearly differs by type of person involved. Drivers and passengers are mostly uninjured or lightly injured, while pedestrians show a much higher proportion of severe and fatal outcomes. This highlights the greater vulnerability of pedestrians in traffic accidents, especially in pedestrian–vehicle collisions. **Question 4. How does injury severity vary by type of vehicle involved?** ```{r plot4} p4 = db %>% count(tipo_vehiculo, lesividad_cat) %>% ggplot(aes(x = tipo_vehiculo, y = n, color = lesividad_cat)) + geom_point(size = 3, alpha = 0.8) + theme_minimal() + labs(title = "Injury severity by vehicle type (dot plot)", x = "Type of vehicle", y = "Number of accidents", color = "Injury severity") + theme(plot.title = element_text(size = 16, face = "bold"), axis.text.x = element_text(angle = 45, hjust = 1)) ggplotly(p4) ``` Differences in accident severity can be observed depending on the type of vehicle, although mild or no injuries are predominant in almost all cases. Most accidents involve passenger cars, which concentrate the highest number of incidents. **Question 5. How do different types of traffic accidents vary across time slots?** ```{r plot5} p5 = db %>% count(time_slot, tipo_accidente) %>% ggplot(aes(x = time_slot, y = n, color = tipo_accidente, group = tipo_accidente)) + geom_line(linewidth = 1) + geom_point(size = 2) + theme_minimal() + labs(title = "Accident types across different time slots", x = "Time slot", y = "Number of accidents", color = "Accident type") + theme(plot.title = element_text(size = 16, face = "bold")) ggplotly(p5) ``` Accident frequency varies clearly across time slots. Most accident types peak during the afternoon and decrease sharply at night, suggesting that higher traffic volume during daytime hours leads to more accidents, regardless of the accident type. **Question 6. How are traffic accidents distributed across districts?** ```{r plot6} p6 = db %>% count(distrito) %>% ggplot(aes(x = reorder(distrito, n), y = n)) + geom_segment(aes(xend = distrito, y = 0, yend = n), color = "grey") + geom_point(color = "steelblue", size = 3) + coord_flip() + theme_minimal() + labs(title = "Number of accidents by district", x = "District", y = "Number of accidents") + theme(plot.title = element_text(size = 16, face = "bold")) ggplotly(p6) ``` Accident distribution varies across districts, with some areas showing a much higher concentration of incidents, likely related to traffic density and urban activity. **Question 7. How does alcohol involvement vary across age groups?** ```{r plot7} p7 = db %>% filter(rango_edad != "Desconocido") %>% group_by(rango_edad) %>% summarise(prop_yes = mean(positiva_alcohol == "Yes"), .groups = "drop") %>% ggplot(aes(x = rango_edad, y = prop_yes, group = 1)) + geom_line(linewidth = 1, color = "darkslateblue") + geom_point(size = 2, color = "darkslateblue") + theme_minimal() + labs(title = "Proportion of alcohol-positive accidents by age group", x = "Age group", y = "Proportion (Alcohol = Yes)") + theme(plot.title = element_text(size = 16, face = "bold")) ggplotly(p7) ``` The proportion of alcohol-positive accidents increases from adolescence and peaks among young adults, then gradually decreases with age. Very low values are observed in children and older age groups, indicating a clear age-related pattern in alcohol involvement. **Question 8. Distribution of types of persons involved by gender** ```{r plot8} p8 = db %>% filter(sexo != "Desconocido") %>% ggplot(aes(x = sexo, fill = tipo_persona)) + geom_bar(position = "fill", alpha = 0.85) + theme_minimal() + labs(title = "Distribution of types of persons involved by gender", x = "Gender", y = "Proportion", fill = "Type of person") + theme(plot.title = element_text(size = 16, face = "bold")) ggplotly(p8) ``` In both genders, most individuals involved in accidents are drivers. However, women show a higher proportion of passengers and pedestrians compared to men, indicating different mobility patterns between genders. # 2. CLASSIFICATION (emphasis on prediction) In this point we are going to use as the base variable lesividad_cat to study if with the other variables of the accident we are able to predict the severity of this one and what patterns seems to appear ```{r part2_setup} set.seed(123) ## Because we are using lesividad_cat we are going to create a subset without ## cod_lesividad and lesividad and other variables that we kept for other analyses, ## but that here are useless to speed up computations as the full dataset made ## some methods take unreasonably long to execute. We are also going to create the ## training and test df = subset(db) df$cod_lesividad = NULL df$lesividad = NULL df$hora = NULL df$fecha = NULL df$localizacion = NULL df$cod_distrito = NULL df$tipo_vehiculo = NULL df$day = NULL df = droplevels(df) idx = createDataPartition(df$lesividad_cat, p = 0.80, list = FALSE) training = df[idx, ] testing = df[-idx, ] ``` ## MACHINE LEARNING ### Decision tree ```{r decision_tree} control = rpart.control(minsplit = 40, maxdepth = 4, cp = 0.001) dtFit = rpart(lesividad_cat~., data = training, method = "class", control = control) rpart.plot(dtFit, digits = 2) ``` #### Prediction: **Define the cost matrix for 4 categories (4x4). As in this case the variable that we want to predict is not binary, the cost matrix has to be n*n, in this case 4*4** ```{r dt_prediction} cost.unit = matrix(c( 0, 50, 200, 500, 10, 0, 100, 300, 50, 20, 0, 100, 100, 50, 20, 0 ), nrow = 4, ncol = 4, byrow = TRUE) #### Predict using the testing data dtPred = predict(dtFit, testing, type = "class") #### Confusion Matrix CM = confusionMatrix(dtPred, testing$lesividad_cat)$table cost = sum(as.vector(CM) * as.vector(cost.unit)) / sum(CM) CM cost ``` The model is biased because it never predicts Grave or Mortal cases. It only predicts the most common categories (Ileso and Leve). Although the average cost is 7.33, the model is not useful for safety because it fails to detect the most dangerous accidents. **Now using Caret** ```{r caret_tree} caret.fit <- train(lesividad_cat~., data = training, method = "rpart", control=rpart.control(minsplit = 40, maxdepth = 12), trControl = trainControl(method = "cv", number = 5), tuneLength=10) caret.fit #### Visualization rpart.plot(caret.fit$finalModel) #### Prediction dtPredCaret = predict(caret.fit, testing) CM_Caret = confusionMatrix(dtPredCaret, testing$lesividad_cat)$table cost_caret = sum(as.vector(CM_Caret) * as.vector(cost.unit)) / sum(CM_Caret) #### Results CM_Caret cost_caret ``` The Caret model reduced the cost to 7.31, but it still fails to predict Grave or Mortal cases. The tree remains biased toward the majority classes (Ileso and Leve) because dangerous accidents are too rare for the algorithm to prioritize ### Random forest ```{r random_forest} rf.train = randomForest(lesividad_cat ~ ., data = training) #### Basic Prediction rf.pred = predict(rf.train, newdata = testing) confusionMatrix(rf.pred, testing$lesividad_cat) ``` **Now let's create a function so Caret does not look for the model with the best accuracy, but the one with the lowest cost** ```{r rf_economic_cost} EconomicCost = function(data, lev = NULL, model = NULL) { y.pred = data$pred y.true = data$obs CM = confusionMatrix(y.pred, y.true)$table out = sum(as.vector(CM) * as.vector(cost.unit)) / sum(CM) names(out) = c("EconomicCost") out } #### Now include this function in the Caret control ctrl = trainControl(method = "cv", number = 5, classProbs = TRUE, summaryFunction = EconomicCost, verboseIter = TRUE) #### Now train a RF using Caret with the specific metric rf.caret = train(lesividad_cat ~ ., method = "rf", data = training, preProcess = c("center", "scale"), ntree = 100, tuneGrid = expand.grid(mtry = c(2, 4, 6)), metric = "EconomicCost", maximize = F, trControl = ctrl) ``` **Variable importance:** ```{r rf_importance} rf_imp = varImp(rf.train, scale = F) plot(rf_imp, scales = list(y = list(cex = .95))) rf_imp ``` Vehicle type is the most critical factor, while alcohol and drugs have minimal impact on the model. Despite this, the Random Forest remains biased, failing to predict any Grave or Mortal cases due to the data imbalance **Final prediction and cost** ```{r rf_final} rfPred = predict(rf.caret, newdata = testing) CM = confusionMatrix(rfPred, testing$lesividad_cat)$table cost = sum(as.vector(CM) * as.vector(cost.unit)) / sum(CM) cost ``` The cost increased to 8.17 because the Random Forest takes more risky decisions, making more expensive mistakes without solving the main issue, it still fails to predict Grave or Mortal cases due to data imbalance. ### Gradient Boosting ```{r gbm} GBM.train = gbm(lesividad_cat ~ ., data = training, distribution = "multinomial", n.trees = 250, shrinkage = 0.01, interaction.depth = 2, n.minobsinnode = 8) #### Prediction and cost gbmProb = predict(GBM.train, newdata = testing, n.trees = 250, type = "response") gbmPred_idx = apply(gbmProb, 1, which.max) levels_order = levels(training$lesividad_cat) gbmPred = factor(levels_order[gbmPred_idx], levels = levels_order) CM_gbm = confusionMatrix(gbmPred, testing$lesividad_cat)$table cost_gbm = sum(as.vector(CM_gbm) * as.vector(cost.unit)) / sum(CM_gbm) cost_gbm ``` The GBM model achieved a cost of 7.40, performing better than Random Forest but still worse than the simple Decision Tree. While Gradient Boosting is more refined, it remains focused on the majority classes and fails to identify Grave or Mortal accidents ### Subsampling techniques In this case we will use SMOTE, because this methods is used to handle imbalanced datasets. And that is the case of this dataset, where the 'ileso' cases dominates the model's learning ```{r smote} ctrl$sampling <- "smote" rf.smote <- train(lesividad_cat ~ ., method = "rf", data = training, preProcess = c("center", "scale"), ntree = 200, tuneGrid = expand.grid(mtry = c(2, 4, 6)), importance = TRUE, metric = "EconomicCost", maximize = FALSE, trControl = ctrl) #### Prediction using the SMOTE-trained model smotePred = predict(rf.smote, newdata = testing) #### Confusion Matrix and costs CM_smote = confusionMatrix(smotePred, testing$lesividad_cat)$table CM_smote cost_smote = sum(as.vector(CM_smote) * as.vector(cost.unit)) / sum(CM_smote) cost_smote ``` SMOTE reduced class imbalance and helped the model detect Grave and Mortal cases that were previously missed. Even though the economic cost increased to 7.73 due to more false positives, this is preferable in road safety, where missing fatal events is riskier than over-predicting severity ## PROBABILISTIC LEARNING. Bayes classifiers ### Naive Bayes ```{r naive_bayes} #### Training Naive Bayes nb.model = naiveBayes(lesividad_cat ~ ., data = training) #### Predicting nb_pred = predict(nb.model, newdata = testing) #### Performance cm_nb = confusionMatrix(nb_pred, testing$lesividad_cat) cm_nb$table cm_nb$overall[1] ``` Naive Bayes achieved the highest accuracy (80.84%) and showed stable performance due to its independence assumption. It avoided the mathematical issues observed in QDA. The model performed better than LDA for Leve injuries but struggled with the Mortal class because of data imbalance ### LDA ```{r lda} #### Training the LDA model lda.model = lda(lesividad_cat ~ ., data = training) lda.model #### Predicting probabilities probability = predict(lda.model, newdata = testing)$posterior head(probability) #### Predicting classes prediction = predict(lda.model, newdata = testing)$class head(prediction) #### Confusion matrix confusionMatrix(prediction, testing$lesividad_cat)$table print(paste("Accuracy:", confusionMatrix(prediction, testing$lesividad_cat)$overall[1])) ``` The model has a high accuracy (0.8), because it tries to minimize the error and tends to classify all as 'Ileso' or 'Leve'. This leads towards a high accuracy model, but that is not able to predict not a single Mortal case ### QDA Training the QDA model: `qda.model <- qda(lesividad_cat ~ ., data = training)` We were unable to perform QDA because the 'Grave' and 'Mortal' classes have very little data. This makes it impossible to calculate their covariance matrices, so we decided to use LDA, Benchmark naive and Naive Bayes ### Benchmark model ```{r benchmark} majority_class = names(which.max(table(training$lesividad_cat))) benchmark_pred = factor(rep(majority_class, nrow(testing)), levels = levels(testing$lesividad_cat)) cm_benchmark = confusionMatrix(benchmark_pred, testing$lesividad_cat) cm_benchmark$table print(paste("Accuracy del Benchmark:", round(cm_benchmark$overall[1], 4))) cost_benchmark = sum(as.vector(cm_benchmark$table) * as.vector(cost.unit)) / sum(cm_benchmark$table) cost_benchmark ``` The Benchmark model achieves a 71.48% accuracy by always predicting the 'Ileso' class. However, as the Confusion Matrix shows, it fails to detect any 'Leve', 'Grave', or 'Mortal' cases. This confirms that accuracy is a misleading metric for imbalanced data and proves that our LDA model is necessary to actually identify accident risks ## Model comparison of all algorithms ```{r comparison} model_comparison <- data.frame( Model = c("Benchmark", "Decision Tree", "Random Forest", "Gradient Boosting (GBM)", "SMOTE RF", "Naive Bayes", "LDA"), Accuracy = c( cm_benchmark$overall["Accuracy"], confusionMatrix(dtPredCaret, testing$lesividad_cat)$overall["Accuracy"], confusionMatrix(rf.pred, testing$lesividad_cat)$overall["Accuracy"], confusionMatrix(gbmPred, testing$lesividad_cat)$overall["Accuracy"], confusionMatrix(smotePred, testing$lesividad_cat)$overall["Accuracy"], cm_nb$overall["Accuracy"], confusionMatrix(prediction, testing$lesividad_cat)$overall["Accuracy"] ), Economic_Cost = c( cost_benchmark, cost_caret, cost, cost_gbm, cost_smote, NA, NA ) ) print(model_comparison) ``` In Machine Learning and in Benchmark model, we use an Economic Cost because accuracy is not enough for road safety. If a model predicts "Ileso" correctly but misses a "Mortal" case, the human cost is very high. Our cost matrix penalizes these dangerous mistakes. In contrast, models like LDA and Naive Bayes are probabilistic baselines used to measure pure statistical accuracy, so they don't include this specific cost optimization. **Conclusions from the Comparison Table:** The Accuracy Paradox: All models have high accuracy (around 80-82%), but their utility is different. The Benchmark has 71% accuracy but the highest cost (16.3), proving that doing nothing is very expensive for society. The Decision Tree and SMOTE RF achieve the best balance and the best performance. They have the lowest costs (7.71 and 7.73), meaning they make the "least dangerous" mistakes. As the final Choice we recommend the SMOTE RF model. Even if the accuracy is slightly lower than the Decision Tree, it is the only model that effectively identifies the most severe accidents by balancing the data. In road safety, identifying a fatal risk is more important than having a high general accuracy # 3. CLASSIFICATION (emphasis on interpretation) We encountered some problems while running this section with the entire database, so we had to take a sample of 12k rows to be able to work comfortably with the data, avoiding crashes and infinite periods of waiting. To address the class imbalance (as severe accidents are much rarer than non-severe ones), we use stratified sampling to ensure our 12000-row subset contains a representative proportion of both classes. This is crucial for the model to actually learn the patterns of "Severe" accidents. ```{r part3_setup} db_clean <- db %>% filter(lesividad_cat != "Desconocido") %>% mutate(target = factor(ifelse(lesividad_cat %in% c("Grave", "Mortal"), "Severe", "NonSevere"), levels = c("NonSevere", "Severe"))) %>% dplyr::select(target, tipo_persona, positiva_alcohol, rango_edad, time_slot, estado_meteorológico, tipo_vehiculo, distrito) %>% na.omit() set.seed(123) # Stratified subset of 12000 observations strat_idx <- createDataPartition(db_clean$target, p = 12000/nrow(db_clean), list = FALSE) db_small <- db_clean[strat_idx, ] # Splitting the small dataset into training (70%) and testing (30%) trainIndex <- createDataPartition(db_small$target, p = 0.7, list = FALSE) train_set <- db_small[trainIndex, ] test_set <- db_small[-trainIndex, ] test_set <- test_set %>% filter(tipo_persona %in% unique(train_set$tipo_persona), tipo_vehiculo %in% unique(train_set$tipo_vehiculo), rango_edad %in% unique(train_set$rango_edad), distrito %in% unique(train_set$distrito)) %>% mutate( tipo_persona = factor(as.character(tipo_persona), levels = levels(train_set$tipo_persona)), tipo_vehiculo = factor(as.character(tipo_vehiculo), levels = levels(train_set$tipo_vehiculo)), rango_edad = factor(as.character(rango_edad), levels = levels(train_set$rango_edad)), distrito = factor(as.character(distrito), levels = levels(train_set$distrito)) ) ``` We fit a logistic regression to interpret the impact of individual predictors. We focus on odds ratios to quantify risk, which is a key requirement for this section. ```{r logit_interp} logit_model <- glm(target ~ tipo_persona + positiva_alcohol + time_slot + rango_edad, data = train_set, family = "binomial") summary(logit_model) # Odds Ratios exp(coef(logit_model)) ``` For the prediction part, we implement a Random Forest with basic hyperparameter tuning using the caret package. We tune 'mtry' (the number of variables randomly sampled at each split) to optimize the model's performance. We will use 'downsampling' so that the model will learn from a balanced data sample, with many 'severe' and 'non severe' cases. If we didn't use this technique, the model would ignore the rare 'severe' cases and just predict 'non-severe' most of the time to get a high accuracy, but it wouldn't be useful for finding the real risks. ```{r rf_tuned} control <- trainControl(method = "cv", number = 5, classProbs = TRUE, summaryFunction = twoClassSummary, sampling = "down") rf_tuned <- train(target ~ ., data = train_set, method = "rf", trControl = control, metric = "ROC", tuneGrid = expand.grid(mtry = c(2, 3, 4))) ``` We evaluate using the area under the curve (AUC) and sensitivity instead of just accuracy. In traffic safety, being able to detect severe cases (sensitivity) is much more important than overall precision. ```{r evaluate_auc} probs_logit <- predict(logit_model, test_set, type = "response") probs_rf <- predict(rf_tuned, test_set, type = "prob")[, "Severe"] roc_logit <- roc(test_set$target, probs_logit) roc_rf <- roc(test_set$target, probs_rf) print(paste("Logit AUC:", auc(roc_logit))) print(paste("RF AUC:", auc(roc_rf))) # Confusion Matrix to check sensitivity for 'severe' cases preds_rf <- predict(rf_tuned, test_set) confusionMatrix(preds_rf, test_set$target) ``` ## SHAP VARIABLE IMPORTANCE (Global Interpretation) **Prepare the data: We take the original predictors from train_set** ```{r shap} X_shap <- train_set[, -1] # Define the prediction wrapper pfun <- function(object, newdata) { predict(object, newdata = newdata, type = "prob")[, "Severe"] } # Calculate SHAP values set.seed(123) shap_val <- fastshap::explain(rf_tuned, X = X_shap, pred_wrapper = pfun, nsim = 30) # Create shapviz object shp_viz <- shapviz(shap_val, X = X_shap) # Plot Global Importance sv_importance(shp_viz, kind = "bar", fill = "darkslateblue") + ggtitle("Global Feature Importance (SHAP Values)") ``` **INTERPRETATION OF THE PLOT:** The SHAP global importance plot reveals that 'tipo_vehiculo' is the most critical feature for our Random Forest model, followed by 'rango_edad' and 'distrito'. Regarding 'positiva_alcohol', its near-zero global importance suggests that while it may be a high-risk factor in specific cases, it does not drive the majority of the model's predictions across the entire dataset. This is likely because the number of accidents where alcohol is positive is very small compared to the total volume of incidents. Therefore, the model relies more on features that are present and vary across all observations, such as the type of vehicle involved or the age of the person. ## CONCLUSIONS In this section, we used stratified sampling because the number of severe accidents is very small compared to the others. The Random Forest model, after tuning the mtry parameter, gives better results for prediction. However, we also used Logistic Regression because it is easier to interpret. To bridge the gap, we implemented SHAP values, which is a key technique from the course to explain "black-box" models. While the Logistic Regression gives us a baseline, the SHAP analysis reveals that 'tipo_vehiculo' and 'rango_edad' are the most important factors for the Random Forest's decisions. This demonstrates how Interpretable Machine Learning allows us to understand complex models beyond simple coefficients. By using AUC and Sensitivity instead of Accuracy, we can better detect the risk in these rare cases and avoid the problems of unbalanced data. ``` ``` # 4. FEATURE SELECTION (and comparison of predictor sets) ## Feature Selection based on SHAP Based on our SHAP importance plot, we select the top 3 predictors: **Vehicle type**, **Age group**, and **District**. These variables show the highest average impact on the model's output. ## Define the Reduced Model ```{r reduced_model} reduced_logit <- glm(target ~ tipo_vehiculo + rango_edad + distrito, data = train_set, family = "binomial") ``` ## Force Test Set to match Training Set Levels ```{r filter_test_levels} test_set_filtered <- test_set %>% filter(tipo_persona %in% unique(train_set$tipo_persona), tipo_vehiculo %in% unique(train_set$tipo_vehiculo), rango_edad %in% unique(train_set$rango_edad), distrito %in% unique(train_set$distrito)) %>% mutate( tipo_persona = factor(as.character(tipo_persona), levels = levels(train_set$tipo_persona)), tipo_vehiculo = factor(as.character(tipo_vehiculo), levels = levels(train_set$tipo_vehiculo)), rango_edad = factor(as.character(rango_edad), levels = levels(train_set$rango_edad)), distrito = factor(as.character(distrito), levels = levels(train_set$distrito)) ) ``` ## Predictions ```{r predictions_comparison} probs_full <- predict(rf_tuned, test_set_filtered, type = "prob")[, "Severe"] probs_red <- predict(reduced_logit, test_set_filtered, type = "response") ``` ## ROC Calculation ```{r roc_comparison} roc_full <- roc(test_set_filtered$target, probs_full) roc_red <- roc(test_set_filtered$target, probs_red) # Print AUC Comparison print(paste("Full Model (Random Forest) AUC:", round(auc(roc_full), 4))) print(paste("Reduced Model (Logistic) AUC:", round(auc(roc_red), 4))) ``` ## Plotting ROC Comparison ```{r plot_roc_comparison} plot(roc_full, col = "blue", main = "ROC Comparison: Full RF vs. Reduced Logit") plot(roc_red, col = "red", add = TRUE) legend("bottomright", legend = c("Full RF", "Reduced Logit"), col = c("blue", "red"), lwd = 2) ``` ## CONCLUSIONS In this section, we performed feature selection using SHAP importance values. Our analysis revealed that **'tipo_vehiculo'** and **'rango_edad'** are the primary predictors of severity. By comparing the complex Random Forest with a reduced Logistic Regression, we can observe if a simpler model maintains high predictive power. This approach is highly valuable for road safety authorities, as it allows them to focus on a few critical factors (like vehicle safety and age-specific vulnerability) instead of managing dozens of irrelevant variables. ```