Analysis and Modelling

Dataset

“Fatalities” dataset is from the AER R package which originally came from US Department of Transportation Fatal Accident Reporting System and is dated from 1982-1988 (The data includes all the states except Alaska and Hawaii)

Link to Fatalities data documentation

Variable Variable Type Description
state Categorical A categorical variable that represents the state
year Categorical The year is indicated by the categorical factor
spirits Continuous A numerical value indicating the amount of spirits consumed.
unemp Continuous The unemployment rate expressed as a numerical value.
income Continuous A numerical value in dollars that represents the per capita personal income. 
emppop Continuous The employment to population ratio expressed as a numerical value.
beertax Continuous The tax on a case of beer expressed as a numerical value.
baptist Continuous A number that indicates the proportion of Southern Baptists in the general population.
mormon Continuous A numerical value that indicates the proportion of Mormons in the general population.
drinkage Continuous A number that represents the minimum age at which alcohol is permitted. 
dry Continuous A numerical number that indicates the proportion of people living in “dry” counties. 
youngdrivers Continuous A numerical value indicating the proportion of drivers between the ages of 15 and 24. 
miles Continuous A numerical value indicating the average number of miles driven by each driver.
breath Boolean A categorical factor that denotes the existence of a law requiring a preliminary breath test. 
jail Boolean A categorical factor that indicates whether a jail sentence is required. 
service Boolean A categorical factor that denotes the existence of a community service requirement. 
fatal Continuous A number that indicates the overall number of people killed in cars. 
nfatal Continuous A numerical value that indicates the quantity of car deaths that occur at night.
sfatal Continuous A numerical value that indicates how many people have died in just one vehicle.
fatal1517 Continuous A numerical value that indicates how many teenagers between the ages of 15 and 17 were killed in cars.
nfatal1517 Continuous A numerical value that indicates how many teenagers, aged 15 to 17, are killed in cars at night.
fatal1820 Continuous A numerical value that indicates the total number of 18–20-year-olds who die in cars.
nfatal1820 Continuous A numerical value that indicates how many teenagers, aged 18 to 20, are killed in cars at night.
fatal2124 Continuous A numerical value that indicates how many people aged 21 to 24 are killed in cars.
nfatal2124 Continuous A numerical value that indicates how many people aged 21 to 24 are killed in cars.
afatal Continuous A numerical value that indicates the quantity of drunk driving deaths.
pop Continuous A numerical figure that indicates the entire population.
pop1517 Continuous A numerical value that indicates the population of people aged 15 to 17.
pop1820 Continuous A numerical value that denotes the 18–20 year old population.
pop2124 Continuous A numerical value that denotes the population of people aged 21 to 24.
milestot Continuous A number (in millions) that indicates the total number of miles driven by a vehicle.
unemppopus Continuous The US unemployment rate expressed as a numerical value. 
Emppopus Continuous The US employment to population ratio expressed as a numerical value. 
gsp Continuous A numerical value that indicates the Gross State Product’s (GSP) rate of change. 
head(df)
  state year spirits unemp   income   emppop  beertax baptist  mormon drinkage
1    al 1982    1.37  14.4 10544.15 50.69204 1.539379 30.3557 0.32829    19.00
2    al 1983    1.36  13.7 10732.80 52.14703 1.788991 30.3336 0.34341    19.00
3    al 1984    1.32  11.1 11108.79 54.16809 1.714286 30.3115 0.35924    19.00
4    al 1985    1.28   8.9 11332.63 55.27114 1.652542 30.2895 0.37579    19.67
5    al 1986    1.23   9.8 11661.51 56.51450 1.609907 30.2674 0.39311    21.00
6    al 1987    1.18   7.8 11944.00 57.50988 1.560000 30.2453 0.41123    21.00
      dry youngdrivers    miles breath jail service fatal nfatal sfatal
1 25.0063     0.211572 7233.887     no   no      no   839    146     99
2 22.9942     0.210768 7836.348     no   no      no   930    154     98
3 24.0426     0.211484 8262.990     no   no      no   932    165     94
4 23.6339     0.211140 8726.917     no   no      no   882    146     98
5 23.4647     0.213400 8952.854     no   no      no  1081    172    119
6 23.7924     0.215527 9166.302     no   no      no  1110    181    114
  fatal1517 nfatal1517 fatal1820 nfatal1820 fatal2124 nfatal2124  afatal
1        53          9        99         34       120         32 309.438
2        71          8       108         26       124         35 341.834
3        49          7       103         25       118         34 304.872
4        66          9       100         23       114         45 276.742
5        82         10       120         23       119         29 360.716
6        94         11       127         31       138         30 368.421
      pop  pop1517  pop1820  pop2124 milestot unempus emppopus         gsp
1 3942002 208999.6 221553.4 290000.1    28516     9.7     57.8 -0.02212476
2 3960008 202000.1 219125.5 290000.2    31032     9.6     57.9  0.04655825
3 3988992 197000.0 216724.1 288000.2    32961     7.5     59.5  0.06279784
4 4021008 194999.7 214349.0 284000.3    35091     7.2     60.1  0.02748997
5 4049994 203999.9 212000.0 263000.3    36259     7.0     60.7  0.03214295
6 4082999 204999.8 208998.5 258999.8    37426     6.2     61.5  0.04897637

  • No of rows and columns

    dim(df)
    [1] 336  34

  • Null and missing values

    print(c(sum(is.na(df)),sum(is.null(df))))
    [1] 2 0

Data Preprocessing

  1. Spirits

  2. Beertax

  3. Baptist

  4. Mormon

  5. Afatal

We believe that the above variables helps to cluster the car crashes due to alcohol consumption

df=df[c(3,7,8,9,26)]

Convert the percentage of Baptists and Mormons to a count

df$baptist=df$baptist/100
df$mormon=df$mormon/100
df$baptist=df$baptist*df$pop
df$mormon=df$mormon*df$pop

Check if there are any missing or null values in our data

print(c(sum(is.na(df)),sum(is.null(df))))
[1] 0 0
describe(df)
        vars   n      mean        sd   median   trimmed      mad    min
spirits    1 336      1.75      0.68     1.67      1.67     0.53   0.79
beertax    2 336      0.51      0.48     0.35      0.42     0.24   0.04
baptist    3 336 353164.70 557450.81 58039.21 239427.34 84461.06   0.00
mormon     4 336  60798.39 160724.27 17809.97  24459.64 17199.78 930.33
afatal     5 336    293.33    303.58   211.59    234.89   200.14  24.60
               max      range skew kurtosis       se
spirits       4.90       4.11 2.20     6.63     0.04
beertax       2.72       2.68 2.17     5.14     0.03
baptist 2879464.51 2879464.51 2.26     6.24 30411.44
mormon  1049096.92 1048166.59 5.04    26.72  8768.23
afatal     2094.90    2070.30 2.70     9.28    16.56

As shown in the boxplot below there are some outliers that need to be removed in our dataset

boxplot(df)

First, check if there are any outliers in our data and remove them before scaling the data. We need to remove the outliers since it can affect the centroid in the k-means clsutering process. In addition, the data needs to be scaled since k-means clustering is reliant on Euclidean distance and a large difference in values can adversly affect it.

detection=function(x) {
  q1=quantile(x,probs=.25) 
  q3=quantile(x,probs=.75) 
  iqr=q3-q1 
  x>q3+(iqr*1.5)|x<q1-(iqr*1.5)}

elimination= function(df,cols=names(df)) {
  for(col in cols) { 
    df=df[!detection(df[[col]]),]}
  return(df)} 
column_names=c('spirits','baptist','mormon', 'beertax', 'afatal')
df2=df
cleaned_df_2=elimination(df2,column_names)
cleaned_df2=as.data.frame(scale(cleaned_df_2))
boxplot(cleaned_df2)

detection=function(x) {
  q1=quantile(x,probs=.25) 
  q3=quantile(x,probs=.75) 
  iqr=q3-q1 
  x>q3+(iqr*1.5)|x<q1-(iqr*1.5)}

elimination= function(df,cols=names(df)) {
  for(col in cols) { 
    df=df[!detection(df[[col]]),]}
  return(df)} 
column_names=c('spirits','baptist','mormon', 'beertax', 'afatal')
df2=df
cleaned_df_2=elimination(df2,column_names)
cleaned_df2=as.data.frame(scale(cleaned_df_2))
boxplot(cleaned_df2)
boxplot(cleaned_df2)

Optimal Number of Clusters

The elbow method is performed and as shown by the graph, the optimal number of clusters is 6.

fviz_nbclust(cleaned_df2, kmeans, method = "wss")

According to this graph the optimal number of clusters could be either 10 or 6. These numbers were chosen since they are the closest to 1 for average silhouette width.

fviz_nbclust(cleaned_df2, kmeans, method = "silhouette")

Lastly we are using the gap statistic as a final check to determine the best number of clusters. This graph also shows 6 clusters as being the most optimal.

fviz_nbclust(cleaned_df2, kmeans, method = "gap_stat")

Modelling with K-means Clustering

Building the model using our variables of interest with the selected number of clusters as 6

custom_colors=c("red","green", "yellow", "orange", "blue", "cyan")
km <- kmeans(cleaned_df2, centers= 6, nstart = 25)
fviz_cluster(km,data=cleaned_df2, geom = "point")+scale_color_manual(values = custom_colors)+scale_fill_manual(values = custom_colors)