Classification and Regression Analysis using British Seatbelt Data

Where to find datasets?: The “datasets” library

library(datasets)

seatbelts <- datasets::Seatbelts

df <- as.data.frame(seatbelts)
# df$Date <- rownames(seatbelts)

df

Load time series modeling library (XTS and Zoo) and using them to decode the data column in the “datasets::Seatbelts” dataset

library(zoo)
library(xts)
# df$Date <- as.yearmon(time(seatbelts))
df$Date <- as.Date(time(seatbelts))

df$law <- as.factor(df$law)
sapply(df,class)

DriversKilled       drivers         front          rear           kms 
    "numeric"     "numeric"     "numeric"     "numeric"     "numeric" 
  PetrolPrice     VanKilled           law          Date 
    "numeric"     "numeric"      "factor"        "Date" 

Exploratory analysis of the difference between classes of law = 1 and law = 0 w.r.t. DriversKilled and VanKilled

boxplot(DriversKilled ~ law, data = df)

t.test(DriversKilled ~ law, data = df) #Establish the statistical significance of the difference of means between law = 1 and law = 0 w.r.t. DriversKilled

    Welch Two Sample t-test

data:  DriversKilled by law
t = 5.1253, df = 29.609, p-value = 1.693e-05
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 15.39892 35.81899
sample estimates:
mean in group 0 mean in group 1 
       125.8698        100.2609 

boxplot(VanKilled ~ law, data = df)

t.test(VanKilled ~ law, data = df) #Establish the statistical significance of the difference of means between law = 1 and law = 0 w.r.t. VanKilled

    Welch Two Sample t-test

data:  VanKilled by law
t = 9.4624, df = 47.965, p-value = 1.508e-12
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 3.474406 5.349366
sample estimates:
mean in group 0 mean in group 1 
       9.585799        5.173913 

Logistic Regression. Use “binomial” to set it as categorical data

model <- glm (formula = law ~ DriversKilled + VanKilled, family = "binomial", data = df)
summary(model)

exp(model$coefficients) # 

Create predictions of the fit logistic model against our original data - each row being an independent evalutation of the model

predictionsOfLawInEffect <- predict(model, newdata = df)
# Categorize the predicted 
predictionsOfLawInEffect_CATEGORICAL <- predictionsOfLawInEffect > -2

print(table(df$law, predictionsOfLawInEffect_CATEGORICAL)) #Confusion matrix

   predictionsOfLawInEffect_CATEGORICAL
    FALSE TRUE
  0   133   36
  1     4   19

Visualize a probability density plot of the predicted odds of the response of law = 1 in the predicted set of evaluated records

ODDS = p / (1 - p); where ODDS = exp(predictionsOfLawInEffect)

=> (1 - p)*ODDS = p

=> p = ODDS/(1+ODDS)

plot(density(exp(predictionsOfLawInEffect))) 

plot(density(exp(predictionsOfLawInEffect)/(1+exp(predictionsOfLawInEffect))))  # plot of the "ODDS" of each row being classified for the law being in effect

print(summary(exp(predictionsOfLawInEffect)/(1+exp(predictionsOfLawInEffect))))

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.0000343 0.0030743 0.0338102 0.1197917 0.1595388 0.8213848 

Boxplot to find the optimal probability cutoff for dichotimization of the law = 1 status

boxplot(exp(predictionsOfLawInEffect)/(1+exp(predictionsOfLawInEffect)) ~ df$law)

Perform Linear Regression as an example use-case to predicted drivers that died as a function of the law status being in effect

# Linear Regression for REGRESSION of the number of people that died ON the Xs
#  /regressors law-being-in-effect and total number of drivers
# Y = number of people that died
# X = [law, drivers, PetrolPrice]

model2 <- lm (DriversKilled ~ law, data = df)
summary (model2)

Call:
lm(formula = DriversKilled ~ law, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-46.870 -17.870  -5.565  14.130  72.130 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  125.870      1.849  68.082  < 2e-16 ***
law1         -25.609      5.342  -4.794 3.29e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 24.03 on 190 degrees of freedom
Multiple R-squared:  0.1079,    Adjusted R-squared:  0.1032 
F-statistic: 22.98 on 1 and 190 DF,  p-value: 3.288e-06

Export df to a file:

write.csv(df, row.names = F, file = "britishSeatBeltStudy.csv")