Homework4_EK
2023-01-20
mydata <- read_excel("C:/Users/Eneja/Desktop/MVA/Electircal_Cars.xlsx")
head(mydata)## # A tibble: 6 × 9
## Brand Model `Capacity (kWh)` `Acceleration (sec)` TopSpeed (km…¹ Range…² Effic…³ Drive Price…⁴
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 Opel Ampera-e 58 7.3 150 335 173 Fron… 42990
## 2 Renault Kangoo Maxi ZE 33 31 22.4 130 160 194 Fron… 0
## 3 Nissan Leaf 36 7.9 144 220 164 Fron… 29990
## 4 Audi e-tron Sportback 55 quattro 86.5 5.7 200 375 231 All … 0
## 5 Porsche Taycan Turbo S 83.7 2.8 260 390 215 All … 186336
## 6 Nissan e-NV200 Evalia 36 14 123 165 218 Fron… 43433
## # … with abbreviated variable names ¹`TopSpeed (km/h)`, ²`Range (km)`, ³`Efficiency (Wh/km)`, ⁴`PriceinGer (€)`colnames(mydata) <- c("Brand","Model","Capacity","Acceleration","Top_Speed","Range","Efficiency","Drive", "Price_Germany", "ID")mydata$DriveF <- factor(mydata$Drive,
levels = c( "Front Wheel Drive", "Rear Wheel Drive", "All Wheel Drive"),
labels = c("Front","Rear","All"))mydata1 <- mydata %>%
replace_with_na(replace = list(Price_Germany = c(0)))mydata1 <- tidyr::drop_na(mydata1)
mydata1$ID <- 1:nrow(mydata1)
mydata1 <- mydata1[c(11, 3, 4, 6, 7, 8 ,9 , 10)]
mydata1 <- mydata1[c(-17, -55), ]Description: -Brand: Vehicle manufacturer -Model: Vehicle model -Capacity: Battery storage capacity in kWh -Acceleration: time to reach 100 km/h in seconds -Top_Speed: Maximum travel speed in km/h -Range: Nominal value of maximal traveled distance in km -Efficiency: Ratio between capacity and range in Wh/km -Drive: type of drive front/rear/all wheel drive -Price_Germany: Estimated value of vehicle in Germany in € in 2022
Unit of observation: BEV car model in Germany in 2022
Sample size: 168 (without missing values)
Source of the data: Kaggle- Cheapest Electric Cars
The main goal of the data analysis is to see if there is a possibility to cluster 168 BEV var vehicles based on 4 characteristics and find out the differences.
psych::describeBy(mydata1[c(-1,-6,-8)])## Warning in psych::describeBy(mydata1[c(-1, -6, -8)]): no grouping variable requested## vars n mean sd median trimmed mad min max range skew kurtosis se
## Capacity 1 166 63.95 20.57 66.5 64.21 24.98 16.7 120 103.3 -0.08 -0.46 1.60
## Acceleration 2 166 7.84 2.98 7.5 7.75 3.04 2.1 15 12.9 0.31 -0.75 0.23
## Range 3 166 332.11 105.95 340.0 329.63 103.78 95.0 660 565.0 0.24 0.19 8.22
## Efficiency 4 166 194.26 32.73 188.5 191.90 30.39 104.0 281 177.0 0.61 -0.27 2.54
## Price_Germany 5 166 57685.93 30544.33 50207.5 52512.41 18558.45 18460.0 187746 169286.0 1.92 4.17 2370.70mydata1$Price_z <- scale(mydata1$Price_Germany)
mydata1$Range_z <- scale(mydata1$Range)
mydata1$Acceleration_z <- scale(mydata1$Acceleration)
mydata1$Capacity_z <- scale(mydata1$Capacity)library(Hmisc)
rcorr(as.matrix(mydata1[, c("Price_z", "Range_z", "Acceleration_z", "Capacity_z")]),
type = "pearson")## Price_z Range_z Acceleration_z Capacity_z
## Price_z 1.00 0.57 -0.55 0.64
## Range_z 0.57 1.00 -0.73 0.86
## Acceleration_z -0.55 -0.73 1.00 -0.64
## Capacity_z 0.64 0.86 -0.64 1.00
##
## n= 166
##
##
## P
## Price_z Range_z Acceleration_z Capacity_z
## Price_z 0 0 0
## Range_z 0 0 0
## Acceleration_z 0 0 0
## Capacity_z 0 0 0There are high correlations, however we will still continue with analysis. We will now calculate dissimilarities.
mydata1$Dissimilarity <- sqrt(mydata1$Price_z^2 +mydata1$Range_z^2 + mydata1$Acceleration_z^2 + mydata1$Capacity_z^2) #Finding outliershead(mydata1[order(-mydata1$Dissimilarity), ], 10) ## # A tibble: 10 × 13
## ID Capacity Acceleration Range Efficiency Drive Price_…¹ DriveF Price…² Range…³ Accel…⁴ Capac…⁵ Dissi…⁶
## <int> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 83 110 3.2 660 167 All Wheel Drive 140000 All 2.69 3.09 -1.56 2.24 4.93
## 2 126 83.7 2.9 380 220 All Wheel Drive 187746 All 4.26 0.452 -1.66 0.960 4.69
## 3 3 83.7 2.8 390 215 All Wheel Drive 186336 All 4.21 0.546 -1.69 0.960 4.67
## 4 140 108. 4.3 610 177 All Wheel Drive 135529 All 2.55 2.62 -1.19 2.13 4.40
## 5 139 108. 6.2 640 168 Rear Wheel Drive 106374 Rear 1.59 2.91 -0.551 2.13 3.98
## 6 16 60 10 575 104 All Wheel Drive 149000 All 2.99 2.29 0.726 -0.192 3.84
## 7 38 16.7 12.7 95 176 Rear Wheel Drive 19120 Rear -1.26 -2.24 1.63 -2.30 3.81
## 8 115 90 2.1 535 168 All Wheel Drive 126990 All 2.27 1.91 -1.93 1.27 3.76
## 9 125 83.7 3.3 385 217 All Wheel Drive 154444 All 3.17 0.499 -1.53 0.960 3.68
## 10 37 16.7 11.9 95 176 Rear Wheel Drive 21720 Rear -1.18 -2.24 1.36 -2.30 3.68
## # … with abbreviated variable names ¹Price_Germany, ²Price_z[,1], ³Range_z[,1], ⁴Acceleration_z[,1], ⁵Capacity_z[,1],
## # ⁶Dissimilarity[,1]library(factoextra)## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa#Calculating Euclidean distances
distance <- get_dist(mydata1[c("Price_z", "Range_z", "Acceleration_z", "Capacity_z")],
method = "euclidian")
distance2 <- distance^2
fviz_dist(distance2) 
get_clust_tendency(mydata1[c("Price_z", "Range_z", "Acceleration_z", "Capacity_z")], #Hopkins statistics
n = nrow(mydata1) - 1,
graph = FALSE)## $hopkins_stat
## [1] 0.8700753
##
## $plot
## NULLBased on Graph above and Hopkins statistics the data is clusterable. Hopkins statistics need to be bigger than 0.5.
library(dplyr)
WARD <- mydata1[c("Price_z", "Range_z", "Acceleration_z", "Capacity_z")] %>%
get_dist(method = "euclidean") %>%
hclust(method = "ward.D2")
WARD##
## Call:
## hclust(d = ., method = "ward.D2")
##
## Cluster method : ward.D2
## Distance : euclidean
## Number of objects: 166library(factoextra)
fviz_dend(WARD)## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
## Please report the issue at <]8;;https://github.com/kassambara/factoextra/issueshttps://github.com/kassambara/factoextra/issues]8;;>.
library(NbClust)
OptNumber <- mydata1[c("Price_z", "Range_z", "Acceleration_z", "Capacity_z")] %>%
NbClust(distance = "euclidean",
min.nc = 2, max.nc = 10,
method = "ward.D2",
index = "all") 
## *** : The Hubert index is a graphical method of determining the number of clusters.
## In the plot of Hubert index, we seek a significant knee that corresponds to a
## significant increase of the value of the measure i.e the significant peak in Hubert
## index second differences plot.
## 
## *** : The D index is a graphical method of determining the number of clusters.
## In the plot of D index, we seek a significant knee (the significant peak in Dindex
## second differences plot) that corresponds to a significant increase of the value of
## the measure.
##
## *******************************************************************
## * Among all indices:
## * 2 proposed 2 as the best number of clusters
## * 11 proposed 3 as the best number of clusters
## * 4 proposed 4 as the best number of clusters
## * 2 proposed 6 as the best number of clusters
## * 2 proposed 8 as the best number of clusters
## * 2 proposed 10 as the best number of clusters
##
## ***** Conclusion *****
##
## * According to the majority rule, the best number of clusters is 3
##
##
## *******************************************************************Based on the statistics and dendrogram we will use 3 clusters. Even majority of rules suggest that the best number of clusters is 3.
mydata1$ClusterWard <- cutree(WARD, #Assigning each unit to a specific cluster
k = 3)
head(mydata1[c("ID", "ClusterWard")])## # A tibble: 6 × 2
## ID ClusterWard
## <int> <int>
## 1 1 1
## 2 2 2
## 3 3 3
## 4 4 2
## 5 5 2
## 6 6 1Leaders_initial <- aggregate(mydata1,
by = list(mydata1$ClusterWard),
FUN = mean)## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical: returning NALeaders_initial## Group.1 ID Capacity Acceleration Range Efficiency Drive Price_Germany DriveF Price_z Range_z
## 1 1 93.70886 70.08608 6.943038 371.0127 188.9241 NA 49475.85 NA -0.2687925 0.3671971
## 2 2 71.75862 44.38793 10.800000 220.2586 200.0517 NA 41541.93 NA -0.5285434 -1.0556932
## 3 3 88.24138 86.34828 4.365517 449.8276 197.2069 NA 112339.34 NA 1.7893146 1.1110909
## Acceleration_z Capacity_z Dissimilarity ClusterWard
## 1 -0.3014927 0.2984124 1.092086 1
## 2 0.9944122 -0.9509834 2.106822 2
## 3 -1.1675167 1.0890502 2.855088 3K-means (non-hierarhical)
library(factoextra)
K_MEANS <- hkmeans(mydata1[c("Price_z", "Range_z", "Acceleration_z", "Capacity_z")],
k = 3,
hc.metric = "euclidean",
hc.method = "ward.D2")
K_MEANS## Hierarchical K-means clustering with 3 clusters of sizes 79, 63, 24
##
## Cluster means:
## Price_z Range_z Acceleration_z Capacity_z
## 1 -0.1502435 0.4179738 -0.3788982 0.4317738
## 2 -0.5755210 -1.0071964 0.9512133 -0.9810338
## 3 2.0052940 1.2680601 -1.2497283 1.1539582
##
## Clustering vector:
## [1] 1 2 3 2 2 1 1 2 1 2 2 2 3 2 2 3 2 2 2 2 1 1 3 2 2 2 2 1 1 1 2 2 1 3 2 2 2 2 1 1 2 3 3 2 1 2 1 1 1 1 1 1 1 2 1 1 2
## [58] 3 3 1 2 2 1 1 1 1 2 2 1 1 2 2 1 1 1 1 1 1 1 1 3 3 3 2 1 1 2 2 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 3 3 3 3
## [115] 3 2 2 3 1 2 1 3 3 3 1 1 3 1 1 1 1 1 1 1 1 1 3 3 1 1 1 1 2 1 1 1 1 1 1 1 2 2 2 2 1 1 1 1 1 1 2 2 1 1 1 1
##
## Within cluster sum of squares by cluster:
## [1] 77.78146 79.07623 54.53236
## (between_SS / total_SS = 68.0 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss" "betweenss" "size"
## [8] "iter" "ifault" "data" "hclust"The ratio of variability between sum of squares vs total sum squares is in our case 68%. The higher the ratio, the more separated the groups.
fviz_cluster(K_MEANS,
palette = "set1",
repel = FALSE,
ggtheme = theme_classic())
We had two outliers ID 17 and 55 in group 3- we removed them.
mydata1$ClusterK_Means <- K_MEANS$clustertable(mydata1$ClusterWard)##
## 1 2 3
## 79 58 29table(mydata1$ClusterK_Means)##
## 1 2 3
## 79 63 24table(mydata1$ClusterWard, mydata1$ClusterK_Means)##
## 1 2 3
## 1 72 7 0
## 2 2 56 0
## 3 5 0 24From the table1 above we can see what was the initial sample size in all three groups. Second table presents sample sizes after reclassification. Last table shows how the units were classified. From Group1 72 units stayed in the same group, 7 were reclassified to group 2. 56 units stayed in group 2 and 2 were classified in group 1. Lastly, 24 units stayed in group 3 and 5 were reclassified to group 1.
Centroids <- K_MEANS$centers
Centroids## Price_z Range_z Acceleration_z Capacity_z
## 1 -0.1502435 0.4179738 -0.3788982 0.4317738
## 2 -0.5755210 -1.0071964 0.9512133 -0.9810338
## 3 2.0052940 1.2680601 -1.2497283 1.1539582library(ggplot2)
library(tidyr)
Figure <- as.data.frame(Centroids)
Figure$id <- 1:nrow(Figure)
Figure <- pivot_longer(Figure, cols = c(Price_z, Range_z, Acceleration_z, Capacity_z))
Figure$Groups <- factor(Figure$id,
levels = c(1, 2, 3),
labels = c("1", "2", "3"))
Figure$nameFactor <- factor(Figure$name,
levels = c("Price_z", "Range_z", "Acceleration_z", "Capacity_z"),
labels = c("Price_z", "Range_z", "Acceleration_z", "Capacity_z"))
ggplot(Figure, aes(x = nameFactor, y = value)) +
geom_hline(yintercept = 0) +
theme_bw() +
geom_point(aes(shape = Groups, col = Groups), size = 3) +
geom_line(aes(group = id), linewidth = 1) +
ylab("Averages") +
xlab("Cluster variables")+
ylim(-3,3)
From the chart above we can see that Group 2 presents BEV vehicles which have a price, available range and capacity below average. However, they have acceleration above average. Reasons for that is usually being that vehicles like that have smaller battery, which makes the vehicle much lighter in weight and also cheaper. We can also see that Group 1 is the closest to average value for all variables researched compared to other two group. Group 1 is below average in price and acceleration and above average in range and capacity. Group 3 has above average values for all variables researched, except acceleration.
Anova
fit <- aov(cbind(Price_z, Range_z, Acceleration_z, Capacity_z) ~ as.factor(ClusterK_Means), data = mydata1)
summary(fit)## Response 1 :
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(ClusterK_Means) 2 119.159 59.580 211.85 < 2.2e-16 ***
## Residuals 163 45.841 0.281
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response 2 :
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(ClusterK_Means) 2 116.303 58.151 194.65 < 2.2e-16 ***
## Residuals 163 48.697 0.299
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response 3 :
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(ClusterK_Means) 2 105.828 52.914 145.76 < 2.2e-16 ***
## Residuals 163 59.172 0.363
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response 4 :
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(ClusterK_Means) 2 107.32 53.660 151.64 < 2.2e-16 ***
## Residuals 163 57.68 0.354
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1H0: µ(var, cluster1) = µ(var, cluster2) = µ(var, cluster3) H1: At least one µ is different
Based on Anova test we can reject null hypothesis at p < 0.001. Therefore, it means the mean values differ between all 3 clusters. Data is also standardized, which makes it comparable based on F-statistics. We found that Price seems to be the most successful with separating groups.
VALIDATION
Once we validated that clustering variables were well chosen, we use new variables to validate chosen variables and explain them even further.
aggregate(mydata1$Efficiency,
by = list(mydata1$ClusterK_Means),
FUN = "mean")## Group.1 x
## 1 1 194.1519
## 2 2 194.8095
## 3 3 193.1667fit1 <- aov(Efficiency ~ as.factor(ClusterK_Means),
data= mydata1)
summary(fit1)## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(ClusterK_Means) 2 49 24.3 0.022 0.978
## Residuals 163 176705 1084.1Ho: µ(Group 1, Price) = µ(Group 2, Price) = µ(Group 3, Price) H1: µ(Group 1, Price) ≠ µ(Group 2, Price) ≠ µ(Group 3, Price)
Based on p-value we cannot reject null hypothesis. We need to find at least one variable that separates the groups.
#Categorical Variable
chi_square <- chisq.test(mydata1$DriveF, as.factor(mydata1$ClusterK_Means))
chi_square##
## Pearson's Chi-squared test
##
## data: mydata1$DriveF and as.factor(mydata1$ClusterK_Means)
## X-squared = 86.622, df = 4, p-value < 2.2e-16addmargins(chi_square$observed)##
## mydata1$DriveF 1 2 3 Sum
## Front 16 48 0 64
## Rear 26 15 3 44
## All 37 0 21 58
## Sum 79 63 24 166addmargins(round(chi_square$expected, 2)) ##
## mydata1$DriveF 1 2 3 Sum
## Front 30.46 24.29 9.25 64
## Rear 20.94 16.70 6.36 44
## All 27.60 22.01 8.39 58
## Sum 79.00 63.00 24.00 166round(chi_square$res, 2) ##
## mydata1$DriveF 1 2 3
## Front -2.62 4.81 -3.04
## Rear 1.11 -0.42 -1.33
## All 1.79 -4.69 4.36library(effectsize)##
## Attaching package: 'effectsize'## The following object is masked from 'package:psych':
##
## phieffectsize::cramers_v(mydata1$DriveF, mydata1$ClusterK_Means)## Cramer's V (adj.) | 95% CI
## --------------------------------
## 0.50 | [0.40, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].Ho: There is no association between clusters of the categorical variable H1: There is association between clusters of the categorical variable
Based on p-value (p < 0.001) we can reject null hypothesis. We found a new variable and we found differences in drive type, there is a relationship. With standardized residuals we can check the direction of the relationship. In group 1 we can see that there is less than expected amount of Front wheel drive type BEV vehicles (alpha = 0.01). We can also see that there is more than expected amount of Rear and All wheel drive type, however based on the value we cannot say it is statistically different from 0 (less than 1.96). In group 2 we can see that there is more than expected amount of Front wheel drive type BEV vehicles (alpha = 0.001) and less than expected amount of All wheel drive type BEV vehicles (alpha = 0.001). We can also see that there is less than expected amount of Rear wheel drive type, however based on the value we cannot say it is statistically different from 0 (less than 1.96). In group 3 we can see that there is We can also see that there is less than expected amount of Rear wheel drive type, however based on the value we cannot say it is statistically different from 0 (less than 1.96).
Conclusion
Classification of 166 BEV vehicles was based on four standardized variables (Price, Range, Acceleration, Capacity).
For hierarchical clustering, Ward’s algorithm was used and based on the analysis of the dendrogram and the indices analyzing the increase in heterogeneity, it was decised to classify BEV vehicles in 3 groups. The classification was further optimized using the K-means Clustering.
Group 1 contains the most BEV vehicles (47,6%), characterized as being closest to the average in all 4 variables compared to other groups. Nevertheless, Group 1 is below average in price and acceleration and above average in range and capacity. There is also less than expected amount of Front wheel drive type BEV vehicles (alpha = 0.01).
Group 2 is the second biggest group of BEV vehicles (37,9%), characterized as below average in price, available range and capacity. However, they have acceleration above average. There is also more than expected amount of Front wheel drive type BEV vehicles (alpha = 0.001) and less than expected amount of All wheel drive type BEV vehicles (alpha = 0.001). BEV vehicles in this group have the best efficiency on average.
Group 3 is the smallest group of BEV vehicles (14,5%) characterized as above average in values for all variables researched, except acceleration.There is also less than expected amount of Front wheel drive type BEV vehicles (alpha = 0.01) and more than expected amount of All wheel drive type BEV vehicles (alpha = 0.001). BEV vehicles in this group have the worst efficiency on average.