Задание для АнДана
Филатов Михаил
19 мая 2018 г
Это всего лишь мой вариант решения задания для андана
Для начала загузим данные.
df <- na.omit(read.csv('https://raw.githubusercontent.com/AFenin/lulsha/master/cars.csv', sep = ';'))Подгузим библиотеки
library(ggplot2)
library(lmtest)
library(memisc)
library(caret)
library(car)
library(sandwich)
library(ggfortify)
library(dplyr)
library(stringr)
library(factoextra)
library(RColorBrewer)
library(psych)Задание 0.
Для начала почистите данные и уберите все строки с неполными данными.
Что ж. С неполными данными мы вроде справились. На самом деле, если посмотреть даже на первые строки в датасете, то можно заметить необычную цену у одного автомобиля.
knitr::kable(head(df, 5))| ID | make | model | version | months_old | power | sale_type | num_owners | gear_type | fuel_type | kms | price | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 97860 | Porsche | 911 | Carrera 4 S CoupпїЅ | 240 | 210 | classic | 3 | manual | gasoline | 202000 | 999999 |
| 2 | 27821 | Ford | Mustang | Gt500 Cabrio Vendido | 54 | 487 | used | 1 | manual | gasoline | 30000 | 685000 |
| 4 | 98251 | Porsche | 911 | R Unidad 343 De 991-Iva Deducible | 14 | 368 | used | 1 | manual | gasoline | 2800 | 470000 |
| 6 | 98222 | Porsche | 911 | 911 R | 17 | 368 | used | 1 | manual | gasoline | 145 | 445000 |
| 12 | 62294 | Bentley | Continental | Gt Convertible V8 S | 8 | 388 | demo | 1 | automatic | gasoline | 2500 | 263724 |
Цена за 1 автомобиль 999999 слишком странная. Посмотрим на boxplot’ы наших данных.
par(mfrow = c(2,2))
boxplot(df$months_old, main = 'months', col = 'blue')
boxplot(df$power, main = 'power', col = 'red')
boxplot(df$kms, main = 'kms', col = 'yellow')
boxplot(df$price, main = 'price', col = 'green')
Вроде везде выбросов хватает. Поэтому начнем чистить данные, начиная с целевой переменной. А потом займемся остальными.
ind_price <- which(df$price %in% boxplot.stats(df$price)$out) # Возьмем индексы наших выбросов
ndf <- df[-ind_price,] # удалим наблюдения с выбросамиАналогично с возрастом автомобилей.
ind_month <- which(ndf$months_old %in% boxplot.stats(ndf$months_old)$out)
ndf <- ndf[-ind_month,]Теперь с пробегом.
ind_kms <- which(ndf$kms %in% boxplot.stats(ndf$kms)$out)
ndf <- ndf[-ind_kms,]И с мощностью.
ind_power <- which(ndf$power %in% boxplot.stats(ndf$power)$out)
ndf <- ndf[-ind_power,]Посмотрим на категориальные переменные.
(levels(ndf$gear_type))## [1] "" "automatic" "manual" "semi-automatic"(levels(ndf$sale_type))## [1] "" "almost_new" "classic" "demo" "km_0"
## [6] "new" "used"(levels(ndf$fuel_type))## [1] "" "CNG" "diesel" "electric" "etanol" "gasoline"
## [7] "hybrid" "LPG"Уберем наблюдения, где есть значения ’’.
ndf <- ndf[!ndf$gear_type == '',]
ndf <- ndf[!ndf$sale_type == '',]
ndf <- ndf[!ndf$fuel_type == '',]Теперь можно приступить непосредственно к заданиям.
Задание 1.
Выберите подходящие переменные-предикторы и постройте на них линейную регрессию для предсказания цены машины. Аргументируйте, почему вы выбрали именно эти переменные. Интерпретируйте полученные данные. О чем они говорят? Визуализируйте модель с помощью функций пакета ggplot2. Поменяйте как минимум одну из переменных (можно просто удалить или добавить одну из переменных), постройте новую модель и сравните эту модель с предыдущей. Какая модель лучше? Почему?
Сделаем сначала дамми для коробок передач.
ndf <- data.frame(ndf, model.matrix(~gear_type-1,ndf)) # дамми для '' не должна была создаться
(sum(ndf$gear_type.1)) # тут даже не осталось значений с ''## [1] 0ndf <- ndf[-13] # поэтому просто удалим ееТеперь дамми для типа топлива.
ndf <- data.frame(ndf, model.matrix(~fuel_type-1,ndf)) # дамми для '' не должна была создаться
(sum(ndf$fuel_type.1)) # Все как и в предыдущем случае## [1] 0ndf <- ndf[-16] # поэтому просто удалим ееВыберем переменные. Для начала посмотрим на то, как остальные переменные коррелируют с зависимой переменной. Будем выбирать те, которые имеют коэффициент корреляции \(> |0.20|\)
cor(ndf[c(5,6,8,11,12)], method = 'spearman') # Спирман, чтобы игнорировать "масштаб" переменных## months_old power num_owners kms price
## months_old 1.00000000 0.01105506 0.18241095 0.84514052 -0.56612609
## power 0.01105506 1.00000000 0.03466709 0.04916158 0.63838773
## num_owners 0.18241095 0.03466709 1.00000000 0.16721504 -0.09591042
## kms 0.84514052 0.04916158 0.16721504 1.00000000 -0.50441845
## price -0.56612609 0.63838773 -0.09591042 -0.50441845 1.00000000И корреляция для коробок передач и цены.
cor(ndf[12:15])## price gear_typeautomatic gear_typemanual
## price 1.00000000 0.47990686 -0.47510334
## gear_typeautomatic 0.47990686 1.00000000 -0.99397152
## gear_typemanual -0.47510334 -0.99397152 1.00000000
## gear_typesemi.automatic -0.03094886 -0.02819548 -0.08156952
## gear_typesemi.automatic
## price -0.03094886
## gear_typeautomatic -0.02819548
## gear_typemanual -0.08156952
## gear_typesemi.automatic 1.00000000А также корреляция для типа топлива.
cor(ndf[c(12,16:ncol(ndf))])## price fuel_typeCNG fuel_typediesel
## price 1.000000000 -7.916937e-03 0.24331676
## fuel_typeCNG -0.007916937 1.000000e+00 -0.02246534
## fuel_typediesel 0.243316758 -2.246534e-02 1.00000000
## fuel_typeelectric 0.023195625 -2.461708e-04 -0.05269833
## fuel_typeetanol -0.005995149 -7.420382e-05 -0.01588498
## fuel_typegasoline -0.247391004 -4.600639e-03 -0.98486898
## fuel_typehybrid 0.004737796 -5.992573e-04 -0.12828434
## fuel_typeLPG -0.008436500 -1.285315e-04 -0.02751503
## fuel_typeelectric fuel_typeetanol fuel_typegasoline
## price 0.0231956248 -5.995149e-03 -0.247391004
## fuel_typeCNG -0.0002461708 -7.420382e-05 -0.004600639
## fuel_typediesel -0.0526983327 -1.588498e-02 -0.984868983
## fuel_typeelectric 1.0000000000 -1.740645e-04 -0.010792004
## fuel_typeetanol -0.0001740645 1.000000e+00 -0.003253058
## fuel_typegasoline -0.0107920038 -3.253058e-03 1.000000000
## fuel_typehybrid -0.0014057150 -4.237278e-04 -0.026271137
## fuel_typeLPG -0.0003015044 -9.088314e-05 -0.005634757
## fuel_typehybrid fuel_typeLPG
## price 0.0047377965 -8.436500e-03
## fuel_typeCNG -0.0005992573 -1.285315e-04
## fuel_typediesel -0.1282843424 -2.751503e-02
## fuel_typeelectric -0.0014057150 -3.015044e-04
## fuel_typeetanol -0.0004237278 -9.088314e-05
## fuel_typegasoline -0.0262711369 -5.634757e-03
## fuel_typehybrid 1.0000000000 -7.339566e-04
## fuel_typeLPG -0.0007339566 1.000000e+00Переменную sale_type, в расчет брать не стоит, так как она наверняка дублирует переменную months_old. Судя по тому, как потенциальные предикторы коррелируют с ценой машины, можно выделить такой набор:
months_old
power
kms
gear_typeautomatic
gear_typemanual
fuel_typediesel
fuel_typegasoline
Перед тем, как строить саму модель, взглянем на графики, может будет что-то интересное.
par(mfrow = c(1,3))
ggplot(data = ndf, aes(x = months_old, y = price))+
geom_point(alpha = 0.4)+
geom_smooth(col = 'red', method = lm)
ggplot(data = ndf, aes(x = power, y = price))+
geom_point(alpha = 0.4)+
geom_smooth(col = 'red', method = lm)
ggplot(data = ndf, aes(x = kms, y = price))+
geom_point(alpha = 0.4)+
geom_smooth(col = 'red', method = lm)
На графиках отчетливо видна гетероскедастичность, которая вроде выражена не так сильно. Поэтому попробуем сгладить ее последствия логарифмированием цены
ndf$lnprice <- log(ndf$price)Теперь попробуем построить саму модель.
model_1 = lm(lnprice ~ months_old + power + kms + gear_typeautomatic + gear_typemanual + fuel_typediesel + fuel_typegasoline, data = ndf)Посмотрим, что получилось.
summary(model_1)##
## Call:
## lm(formula = lnprice ~ months_old + power + kms + gear_typeautomatic +
## gear_typemanual + fuel_typediesel + fuel_typegasoline, data = ndf)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.54206 -0.14640 -0.00242 0.14304 1.77944
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.635e+00 4.459e-02 193.662 <2e-16 ***
## months_old -6.271e-03 6.812e-05 -92.052 <2e-16 ***
## power 1.269e-02 8.119e-05 156.312 <2e-16 ***
## kms -2.282e-06 4.769e-08 -47.843 <2e-16 ***
## gear_typeautomatic 8.820e-02 3.574e-02 2.468 0.0136 *
## gear_typemanual -4.412e-02 3.553e-02 -1.242 0.2143
## fuel_typediesel 2.307e-01 2.616e-02 8.818 <2e-16 ***
## fuel_typegasoline 1.870e-02 2.649e-02 0.706 0.4803
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2351 on 19052 degrees of freedom
## Multiple R-squared: 0.8232, Adjusted R-squared: 0.8231
## F-statistic: 1.267e+04 on 7 and 19052 DF, p-value: < 2.2e-16В целом мы видим, что в модели есть 5 статистически значимых предикторов. Регрессия объясняет около 82% зависимой переменной. Также мы отвергаем гипотезы о том, что \(R^2\) был получен случайным образом и модель на константу не хуже модели с предикторами. Также стоит проверить страдает ли модель от мультиколлинеарности или гетероскедастичности (предположение о которой мы высказали ранее).
bptest(model_1)##
## studentized Breusch-Pagan test
##
## data: model_1
## BP = 1198.5, df = 7, p-value < 2.2e-16Как мы видим, нам придется отвергнуть гипотезу о том, что дисперсия остатков гомоскедастична. Поэтому попробуем сделать модель, преобразовав зависимую переменную и предикторы с помощью трансформации Бокса-Кокса, иначе все наши оценки несостоятельны, неустойчивы и неэффективны(
bcprice <- BoxCoxTrans(ndf$price)
bcmonths_old <- BoxCoxTrans(ndf$months_old)
bckms <- BoxCoxTrans(ndf$kms)
bcpower <- BoxCoxTrans(ndf$power)
ndf <- cbind(ndf,
new_price = predict(bcprice, ndf$price),
new_month = predict(bcmonths_old, ndf$months_old),
new_kms = predict(bckms, ndf$kms),
new_power = predict(bcpower, ndf$power))Посмтрим, что выйдет теперь.
model_1bc <- lm(new_price ~ new_month + new_power + new_kms + gear_typeautomatic + gear_typemanual + fuel_typediesel + fuel_typegasoline, data = ndf)
bptest(model_1bc)##
## studentized Breusch-Pagan test
##
## data: model_1bc
## BP = 559.25, df = 7, p-value < 2.2e-16Возможно, поможет обычная стандартизация.
scale_df <- cbind(scale(ndf[c(12, 5,6,11)]), ndf[c('gear_typeautomatic', 'gear_typemanual', 'fuel_typediesel', 'fuel_typegasoline')])
model_1scl <- lm(price ~ months_old + power + kms + gear_typeautomatic + gear_typemanual + fuel_typediesel + fuel_typegasoline, data = scale_df)
bptest(model_1scl)##
## studentized Breusch-Pagan test
##
## data: model_1scl
## BP = 930.01, df = 7, p-value < 2.2e-16Что ж, видимо, я плохой борец с гетероскедастичностью. Но давайте остановим свой выбор на второй модели(с преобразованиями Бокса-Кокса), так как статистика теста Бройша-Пагана оказалась там меньше. Конечно, такая модель вряд ли пригодится в реальной в жизни. Возможно стоит воспользоваться некой другой техникой для предсказания цены машины.
autoplot(model_1bc)
Видимо, я с самого начала забыл про влиятельные наблюдения, посмотрим на них в первой модели.
autoplot(model_1)
Обратим на них более пристальное внимание. Тогда, возможно, мы поправим ситуацию с гетероскедстичностью.
influencePlot(model_1, id.method= 'identify')
Теперь явно видно, что есть ряд точек, от которых страдает наша модель.
model_infl_meas <- influence.measures(model_1)
ndf_influential_obs <- which(apply(model_infl_meas$is.inf, 1, any))
ndf_ninf <- ndf[-ndf_influential_obs,]Вроде всё прошло без ошибок. Но потерялись дамми-переменные. Создадим их снова. Затем попробуем сделать нашу модель заново, исключив влиятельные наблюдения, выбросы и точки высокой напряженности.
ndf_ninf <- ndf_ninf[-c(13:22)] #удалили старые дамми, так как они все стали NA
ndf_ninf$gear_automat <- ifelse(ndf_ninf$gear_type == 'automatic', 1, 0)
ndf_ninf$gear_manual <- ifelse(ndf_ninf$gear_type == 'manual', 1, 0)
ndf_ninf$gear_semi <- ifelse(ndf_ninf$gear_type == 'semi-automatic', 1, 0)
sum(ndf_ninf$gear_semi) # оказывается, мы исключили полуавтомат## [1] 0ndf_ninf <- ndf_ninf[-20]
ndf_ninf$fuel_gasoline <- ifelse(ndf_ninf$fuel_type == 'gasoline', 1, 0)
ndf_ninf$fuel_disel <- ifelse(ndf_ninf$fuel_type == 'diesel', 1, 0)
sum(ndf_ninf$fuel_disel)## [1] 15173Manual и Gasoline пропали, я не знаю, как вернуть их. Поэтому дальше без них.
model_1bc_ninf <- lm(new_price ~ new_month + new_power + new_kms + gear_automat + fuel_disel, data = ndf_ninf)
bptest(model_1bc_ninf)##
## studentized Breusch-Pagan test
##
## data: model_1bc_ninf
## BP = 357.83, df = 5, p-value < 2.2e-16от гетероскедастичности мы все-таки не избавились, но значение статистики Бройша-Пагана заметно уменьшилось по сравнению с тем, что мы имели в самом начале. Посмотрим на результат этой модели.
summary(model_1bc_ninf)##
## Call:
## lm(formula = new_price ~ new_month + new_power + new_kms + gear_automat +
## fuel_disel, data = ndf_ninf)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.1740 -2.6598 0.0426 2.6663 11.9519
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 30.3057185 0.2269221 133.55 <2e-16 ***
## new_month -1.4353231 0.0217306 -66.05 <2e-16 ***
## new_power 2.0456532 0.0134740 151.82 <2e-16 ***
## new_kms -0.0223142 0.0007006 -31.85 <2e-16 ***
## gear_automat 3.2450680 0.0758591 42.78 <2e-16 ***
## fuel_disel 3.8270927 0.0830539 46.08 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.911 on 18130 degrees of freedom
## Multiple R-squared: 0.8211, Adjusted R-squared: 0.821
## F-statistic: 1.664e+04 on 5 and 18130 DF, p-value: < 2.2e-16autoplot(model_1bc_ninf)
Сравнивая с тем, что было изначально, оценки коэффициентов по сути дела не изменились, но теперь у нас появилось больше оснований верить резултатам регрессии.
Что ж, теперь попробуем выбросить один предиктор и сравнить модели.
model_2 <- lm(new_price ~ new_month + new_kms + gear_automat + fuel_disel, data = ndf_ninf) # Удалили power
model_3 <- lm(new_price ~ new_month + new_power + new_kms + gear_automat + fuel_disel + num_owners, data = ndf_ninf) # вдруг зря не взяли num_owners
mtable(model_1bc_ninf, model_2, model_3)##
## Calls:
## model_1bc_ninf: lm(formula = new_price ~ new_month + new_power + new_kms + gear_automat +
## fuel_disel, data = ndf_ninf)
## model_2: lm(formula = new_price ~ new_month + new_kms + gear_automat +
## fuel_disel, data = ndf_ninf)
## model_3: lm(formula = new_price ~ new_month + new_power + new_kms + gear_automat +
## fuel_disel + num_owners, data = ndf_ninf)
##
## ================================================================
## model_1bc_ninf model_2 model_3
## ----------------------------------------------------------------
## (Intercept) 30.306*** 61.211*** 30.991***
## (0.227) (0.151) (0.258)
## new_month -1.435*** -1.418*** -1.428***
## (0.022) (0.033) (0.022)
## new_power 2.046*** 2.047***
## (0.013) (0.013)
## new_kms -0.022*** -0.021*** -0.022***
## (0.001) (0.001) (0.001)
## gear_automat 3.245*** 8.513*** 3.251***
## (0.076) (0.102) (0.076)
## fuel_disel 3.827*** 6.633*** 3.818***
## (0.083) (0.122) (0.083)
## num_owners -0.738***
## (0.132)
## ----------------------------------------------------------------
## R-squared 0.821 0.594 0.821
## adj. R-squared 0.821 0.594 0.821
## sigma 3.911 5.894 3.908
## F 16639.287 6620.419 13894.555
## p 0.000 0.000 0.000
## Log-likelihood -50465.079 -57904.325 -50449.387
## Deviance 277329.060 629917.736 276849.569
## AIC 100944.157 115820.650 100914.773
## BIC 100998.797 115867.484 100977.219
## N 18136 18136 18136
## ================================================================На самом деле ситуация довольно спорная – выбрать первую или третью модель. \(R^2\) в них одинаковый, информационные критерии не дают ощутимой разницы, но все же меньше штрафуют третью модель. При этом F-статистика намного больше у первой модели. И все же первая модель кажется более привлекательной. В ней меньше переменных, а переменная num_owners в третьей модели не обладает какой-то огромной дисперсией, которая могла бы объяснить эту модель, поэтому \(R^2\) в итоге не повышается.
Задание 2.
- Разбейте данные на кластеры, используя метод k-means. Не забудьте показать и рассказать, как вы выбирали количество кластеров. Покажите в датафрейме, какие данные к какому кластеру относятся. Бонус-трек – постройте график с визуализацией кластеров.
- Сделайте иерархическую кластеризацию и постройте дендрограмму. Покажите кластеры на дендрограмме.
Начнем с кластеризации методом k-means. Для начала возьмем переменные, которые было бы интересно интерпретировать.
to_clust <- ndf_ninf %>% select(months_old, power, kms, price, num_owners)Нужны ли для k-means стандартизированные данные или нет, насколько мне известно, это спорный вопрос. Но масштаб у наших данных явно различный, поэтому шкалируем наши переменные.
to_clust <- data.frame(scale(to_clust))Теперь начнем подбор необходимого нам числа кластеров. Начнем с “силуэтов” и ими же закончим (не хватает оперативки на метод локтя:-( ).
fviz_nbclust(to_clust, kmeans, method = "silhouette") +
labs(subtitle = "Silhouette method")
Тут мы можем заметить, что 6 кластеров для нас наиболее оптимальны, так как наиболее удобные силэты получаются именно в данном случае, то есть данные лежат однородно внутри каждого кластера. Теперь проведем саму кластеризацию.
kmclust <- kmeans(to_clust, 6)Добавим метки кластеров в датасет.
to_clust$kmcl <- factor(kmclust$cluster)Попробуем визуализировать нашу кластеризацию
set.seed(123)
darkcols <- brewer.pal(8, "Dark2")
fviz_cluster(kmclust, data = to_clust[1:5],
ellipse.type = "convex",
palette = darkcols,
ggtheme = theme_minimal())
Есть еще вариант, на который у меня не хватило оперативной памяти.
#disttoclust <- dist(to_clust[1:5])
#toclustmds <- cmdscale(disttoclust)
#
#plot(toclustmds, col = kmclust$cluster)Глядя на этот график, трудно говорить хорошим ли вышло деление или нет, но очень сильно выбивается первый кластер, он будто “заполз” под остальные на этом графике.
Интересно посмотреть, как же все кластеризовалось все по месяцам и по цене (можно было бы и дальше, но наверное смысла нет).
ggplot(to_clust, aes(x = months_old, y = price))+
geom_point(col = to_clust$kmcl, alpha = 0.3)
Тоже не особо круто, но в принципе, наверное, интерпретируемо.
Можно прилепить метки к исходному датасету и посмотреть, какие там средние и попробовать интерпретировать.
ndf_ninf$kmcl <- factor(kmclust$cluster)
ndf_ninf %>% group_by(kmcl) %>% summarise_at(vars(5,6,8,11,12), mean)## # A tibble: 6 x 6
## kmcl months_old power num_owners kms price
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 57.6 71.6 1.00 87833. 9743.
## 2 2 58.8 115. 1.00 100653. 17656.
## 3 3 121. 92.5 1.00 154631. 7427.
## 4 4 90.2 97.6 2.12 124428. 12026.
## 5 5 19.7 119. 1.00 21977. 29485.
## 6 6 17.6 80.8 1.00 19625. 15606.Можно заметить, что есть кластеры с явно подержаными и старыми дешевыми машинами, есть довольно новенькие, а есть что-то среднее между ними.
Теперь реализуем иерархическую кластеризацию. для начала найдем матрицу расстояний. Будем пользоваться квадратом Евклидова расстояния, использовать для агрегации метод Варда. Также мне не хватает оперативки, поэтому я сэмплирую 15000 наблюдений.
set.seed(123)
to_clust_sampl = sample(to_clust, 15000)
hc1 <- hclust(dist(to_clust_sampl[1:5]), method = 'ward.D' )построим дендограмму и отметим на ней кластеры.
plot(hc1, main = 'Дендограмма на 15000 наблюдениях')
rect.hclust(hc1, k = 5)
Три или четыре кластера было бы совсем мало – мы бы не смогли разделить различия в левых ветвях.
Задание 3.
Выберете подходящие переменные для эксплораторного факторного анализа. Постройте на них матрицу корреляций. Что можете сказать? Примените метод главных компонент (principal components analysis). И снова, обоснуйте, какое количество компонент кажется оптимальным и почему? Приведите процент объясненной дисперсии первыми двумя компонентами. Визуализируйте результаты PCA. Что говорит вам PCA? Проинтерпретируйте результаты.
Возьмем все количественные переменные, а дальше посмотрим, что окажется для нас важным.
fa_pca <- ndf_ninf[c(5,6,8,11,12)]Построим матрицу корреляций для этих переменных.
cor_m <- cor(fa_pca)
cor_m## months_old power num_owners kms price
## months_old 1.00000000 0.01814416 0.1601872 0.78820986 -0.5459114
## power 0.01814416 1.00000000 0.0308881 0.06886524 0.6521853
## num_owners 0.16018717 0.03088810 1.0000000 0.14454173 -0.0701031
## kms 0.78820986 0.06886524 0.1445417 1.00000000 -0.4902639
## price -0.54591140 0.65218531 -0.0701031 -0.49026390 1.0000000Можно заметить, все переменные коррелируют между собой не сильно или средне, высокая корреляция наблюдается только между пробегом машины и возрастом машины.
Применим метод главных компонент к имеющимся данным.
PC <- prcomp(fa_pca, scale = TRUE )
(PC)## Standard deviations (1, .., p=5):
## [1] 1.5263245 1.1751882 0.9705387 0.4587816 0.3699191
##
## Rotation (n x k) = (5 x 5):
## PC1 PC2 PC3 PC4 PC5
## months_old 0.5652099 -0.2971246 0.152185796 -0.68092323 -0.324711888
## power -0.2326831 -0.7619989 0.162577768 -0.09791146 0.573757862
## num_owners 0.1510618 -0.2461426 -0.957237518 0.01496516 0.008156645
## kms 0.5439464 -0.3453904 0.184535212 0.72310688 -0.167005716
## price -0.5547062 -0.3888367 0.007144449 0.06041013 -0.733080259summary(PC)## Importance of components:
## PC1 PC2 PC3 PC4 PC5
## Standard deviation 1.5263 1.1752 0.9705 0.4588 0.36992
## Proportion of Variance 0.4659 0.2762 0.1884 0.0421 0.02737
## Cumulative Proportion 0.4659 0.7421 0.9305 0.9726 1.00000Нам хватает двух компонент, так как их собственные значения больше единицы, а объясняют они большую часть дисперсии. Первые две компоненты объяснили около 74% дисперсии, что в принципе достаточно много. Собственные значения должны превышать 1 по правилу Кайзера-Харриса.
Теперь визуализируем компоненты.
fviz_pca_var(PC,
col.var = "contrib",
gradient.cols = c("#00AFBB", "#E7B800"),
repel = TRUE
)
Логично было предположить, что переменная num_owners вносит малый вклад в данные, так как ее дисперсия крайне мала(наблюдения принимают значения от 1 до 11). Остальные же переменные вносят довольно большой вклад в данные и их можно использовать в нашем анализе далее. Первая компонента отвечает за “изношенность” (months_old и kms) и цену, а вторая за “лошадок” в машине(power).
Бонус-трек Вообще, тему с PCA я еще до конца так и не просек, но постараюсь визуализировать то, что предлагается.
PC1 <- principal(cor_m, nfactors=2, rotate= 'none') # Возьмем только первые две компоненты с помощью пакета psychset.seed(123)
pred <- predict(PC1, fa_pca) # преобразованные PCA данные
kmpca <- kmeans(pred, 6)
fviz_cluster(kmpca, data = pred,
ellipse.type = "convex",
palette = darkcols,
ggtheme = theme_minimal())
С PCA в универе я еще не знакомился и делаю это задание в рассвет перед дедлайном… НО как я понял фишка в том, чтобы мы могли уловить разброс/дисперсию в данных и спроецировать или найти для него ось так, чтобы он стал для нас более информативен. В итоге я кластерезировал данные предсказанные на основе PCA. Кластеры не приняли какую-то новую форму, а остались такими же, но теперь они спроецировались в двумерное пространство так, что теперь тот “распластавшийся” кластер, который напрягал меня в предыдущем задании очень красиво отобразился на плоскость и не перекрывается другими кластерами.