Practice10
A.Lukyanova
12 мая 2018 г
В работе участвует набор данных Auto. В этом наборе:
-mpg - расход топлива (галлонов на милю)
-cylinders - количество цилиндров (между 4 и 8)
-displacement - объем двигателя в квадратных дюймах
-horsepower - лошадиная сила
-weight - масса автомобиля в кг
-acceleration - время разгона от 0 до 60 миль/час
-year - год выпуска
-origin - страна-производитель
Посчитаем основные характеристики.
data("Auto")
str(Auto)## 'data.frame': 392 obs. of 9 variables:
## $ mpg : num 18 15 18 16 17 15 14 14 14 15 ...
## $ cylinders : num 8 8 8 8 8 8 8 8 8 8 ...
## $ displacement: num 307 350 318 304 302 429 454 440 455 390 ...
## $ horsepower : num 130 165 150 150 140 198 220 215 225 190 ...
## $ weight : num 3504 3693 3436 3433 3449 ...
## $ acceleration: num 12 11.5 11 12 10.5 10 9 8.5 10 8.5 ...
## $ year : num 70 70 70 70 70 70 70 70 70 70 ...
## $ origin : num 1 1 1 1 1 1 1 1 1 1 ...
## $ name : Factor w/ 304 levels "amc ambassador brougham",..: 49 36 231 14 161 141 54 223 241 2 ...Auto1 <- Auto[,-9]
summary(Auto1)## mpg cylinders displacement horsepower
## Min. : 9.00 Min. :3.000 Min. : 68.0 Min. : 46.0
## 1st Qu.:17.00 1st Qu.:4.000 1st Qu.:105.0 1st Qu.: 75.0
## Median :22.75 Median :4.000 Median :151.0 Median : 93.5
## Mean :23.45 Mean :5.472 Mean :194.4 Mean :104.5
## 3rd Qu.:29.00 3rd Qu.:8.000 3rd Qu.:275.8 3rd Qu.:126.0
## Max. :46.60 Max. :8.000 Max. :455.0 Max. :230.0
## weight acceleration year origin
## Min. :1613 Min. : 8.00 Min. :70.00 Min. :1.000
## 1st Qu.:2225 1st Qu.:13.78 1st Qu.:73.00 1st Qu.:1.000
## Median :2804 Median :15.50 Median :76.00 Median :1.000
## Mean :2978 Mean :15.54 Mean :75.98 Mean :1.577
## 3rd Qu.:3615 3rd Qu.:17.02 3rd Qu.:79.00 3rd Qu.:2.000
## Max. :5140 Max. :24.80 Max. :82.00 Max. :3.000Проведем стандартизацию данных и создадим новый массив данных для дальнейшего исследования.
#стандартизация переменных
mpg1 <- (Auto1$mpg-mean(Auto1$mpg))/sd(Auto1$mpg)
cylinders1 <- (Auto1$cylinders-mean(Auto1$cylinders))/sd(Auto1$cylinders)
displacement1 <- (Auto1$displacement-mean(Auto1$displacement))/sd(Auto1$displacement)
horsepower1 <- (Auto1$horsepower-mean(Auto1$horsepower))/sd(Auto1$horsepower)
weight1 <- (Auto1$weight-mean(Auto1$weight))/sd(Auto1$weight)
acceleration1 <- (Auto1$acceleration-mean(Auto1$acceleration))/sd(Auto1$acceleration)
year1 <- (Auto1$year-mean(Auto1$year))/sd(Auto1$year)
origin1 <- (Auto1$origin-mean(Auto1$origin))/sd(Auto1$origin)
#собираем данные в новую матрицу
auto <- data.table(mpg1, cylinders1, displacement1, horsepower1, weight1, acceleration1, year1, origin1)После стандартизации выполним анализ методом главных компонент.
#метод главных компонент
apply(auto, 2, mean)## mpg1 cylinders1 displacement1 horsepower1 weight1
## 1.569283e-16 -9.915610e-17 -5.786680e-17 -1.767977e-16 -8.020330e-18
## acceleration1 year1 origin1
## -2.319512e-16 -1.271467e-15 1.321042e-16apply(auto, 2, var)## mpg1 cylinders1 displacement1 horsepower1 weight1
## 1 1 1 1 1
## acceleration1 year1 origin1
## 1 1 1pr.out=prcomp(auto, scale=TRUE)
names(pr.out)## [1] "sdev" "rotation" "center" "scale" "x"pr.out$center## mpg1 cylinders1 displacement1 horsepower1 weight1
## 1.554702e-16 -9.969350e-17 -5.836105e-17 -1.764373e-16 -8.903734e-18
## acceleration1 year1 origin1
## -2.321299e-16 -1.272013e-15 1.319806e-16pr.out$scale## mpg1 cylinders1 displacement1 horsepower1 weight1
## 1 1 1 1 1
## acceleration1 year1 origin1
## 1 1 1pr.out$rotation## PC1 PC2 PC3 PC4 PC5
## mpg1 0.3858624 0.07663269 -0.2922857863 0.09998251 -0.74036644
## cylinders1 -0.4023885 0.13842878 -0.0722393515 -0.21603551 -0.48261485
## displacement1 -0.4164443 0.12632499 -0.0742362247 -0.13581398 -0.30331627
## horsepower1 -0.4018359 -0.11148007 -0.2360557124 -0.11971643 0.08426839
## weight1 -0.4015758 0.21102000 0.0008939861 -0.32246785 0.13127292
## acceleration1 0.2647309 0.41690206 0.6394351361 -0.49280794 -0.09773197
## year1 0.2138678 0.69046320 -0.5871891991 -0.10601968 0.30134385
## origin1 0.2778682 -0.50150064 -0.3073238207 -0.74328281 0.04739508
## PC6 PC7 PC8
## mpg1 0.387351653 -0.19588516 -0.11513210
## cylinders1 -0.530925481 0.27878265 -0.41774679
## displacement1 -0.006997047 -0.08422855 0.82916553
## horsepower1 0.666709602 0.53504996 -0.13477548
## weight1 0.235859612 -0.72202073 -0.30991105
## acceleration1 0.202933431 0.22891382 0.03518826
## year1 -0.110025917 0.12501506 0.05428840
## origin1 -0.120866629 -0.03452660 0.07951102dim(pr.out$x)## [1] 392 8biplot(pr.out, scale=0)
pr.out$rotation=-pr.out$rotation
pr.out$x=-pr.out$x
biplot(pr.out, scale=0)
pr.out$sdev## [1] 2.3185927 0.9714233 0.9009127 0.6972488 0.4275822 0.3381153 0.2314019
## [8] 0.1788003pr.var=pr.out$sdev^2
pr.var## [1] 5.37587230 0.94366326 0.81164365 0.48615594 0.18282657 0.11432193
## [7] 0.05354682 0.03196954pve=pr.var/sum(pr.var)
pve## [1] 0.671984038 0.117957908 0.101455456 0.060769492 0.022853321 0.014290241
## [7] 0.006693352 0.003996192plot(pve, xlab="Principal Component", ylab="Proportion of Variance Explained", ylim=c(0,1),type='b')
plot(cumsum(pve), xlab="Principal Component", ylab="Cumulative Proportion of Variance Explained", ylim=c(0,1),type='b')
После анализа получаем, что нам лучше оставить 2 первых главных компоненты, так как суммарно они объясняют около 80% всей дисперсии.
Выполним кластерный анализ методом k-средних.
km.out=kmeans(auto,2,nstart=20)
km.out$cluster## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 2 2 2 2 2 2 1 1 1 1 1 2 2 2 1 1 1
## [36] 1 1 1 1 1 1 1 1 1 1 2 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1
## [71] 2 1 1 1 1 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1
## [106] 1 1 2 2 2 2 2 1 2 1 1 2 2 2 2 1 2 2 1 1 1 1 2 2 2 2 1 1 1 1 1 1 1 1 2
## [141] 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 1 2 2 2 2 1 2 1
## [176] 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 2 1 2 2 2 2 1 1 1 1 2 2 2 2 2 1 2 2 1
## [211] 1 1 1 1 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [246] 2 2 1 1 1 1 2 2 2 1 1 1 2 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 1 2 1 2 2 1 2
## [281] 2 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 1 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1
## [316] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [351] 2 2 2 2 2 2 2 2 1 1 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2
## [386] 2 2 2 2 2 2 2plot(auto, col=(km.out$cluster+1), main="K-Means Clustering Results with K=2", pch=20, cex=2)
km.out=kmeans(auto,3,nstart=20)
km.out$cluster## [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 3 3 3 1 1 1 1 1 1 3 2 2 2 2 1 3 1 3 3 3
## [36] 3 3 2 2 2 2 2 2 2 3 3 3 3 3 1 1 1 1 1 1 3 1 3 1 3 3 2 2 2 2 2 2 2 2 2
## [71] 1 2 2 2 2 3 1 3 1 3 1 1 3 1 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 1 2 2 2
## [106] 2 3 1 3 1 1 3 3 1 2 2 1 1 1 3 2 1 1 2 3 3 3 1 3 1 3 3 3 3 2 2 2 2 2 1
## [141] 1 1 1 1 1 1 1 1 1 1 3 3 3 3 2 2 2 2 3 3 3 3 3 2 2 1 3 3 3 1 1 1 3 1 3
## [176] 1 1 3 1 1 1 1 3 3 1 2 2 2 2 3 3 3 3 1 3 1 1 3 3 3 3 1 1 1 3 3 2 3 1 3
## [211] 2 2 2 2 1 1 1 3 1 2 3 2 2 3 3 3 3 2 2 2 2 1 3 1 3 1 1 1 1 1 1 1 1 1 1
## [246] 1 1 3 2 2 3 3 3 3 3 3 3 3 3 3 2 2 2 2 1 1 1 1 1 3 3 1 3 3 3 3 1 1 3 3
## [281] 3 3 3 2 2 2 2 2 2 2 2 1 1 1 3 3 2 1 3 1 1 1 1 3 3 3 1 1 1 1 1 3 3 3 3
## [316] 1 1 1 1 1 1 3 1 1 1 1 1 1 1 1 1 1 1 1 3 3 3 3 1 1 1 1 1 1 1 1 1 1 1 1
## [351] 1 1 1 1 1 1 1 1 3 3 3 3 3 3 1 1 3 3 3 1 1 1 1 1 1 1 1 1 1 1 3 3 3 3 1
## [386] 1 3 3 1 1 3 3plot(auto, col=(km.out$cluster+1), main="K-Means Clustering Results with K=3", pch=20, cex=2)
km.out=kmeans(auto,2,nstart=20)
km.out$tot.withinss## [1] 1553.797km.out=kmeans(auto,3,nstart=20)
km.out$tot.withinss## [1] 1167.322Проанализировав внутриклассовые дисперсии, приходим к выводу, что наиболее оптимальным является разбиение на 3 кластера с ядром равным 20.
Выполним иерархическую кластеризацию.
hc.complete=hclust(dist(auto), method="complete")
hc.average=hclust(dist(auto), method="average")
hc.single=hclust(dist(auto), method="single")
par(mfrow=c(1,3))
plot(hc.complete,main="Complete Linkage", xlab="", sub="", cex=.9)
plot(hc.average, main="Average Linkage", xlab="", sub="", cex=.9)
plot(hc.single, main="Single Linkage", xlab="", sub="", cex=.9)
par(mfrow=c(1,1))
cutree(hc.complete, 2)## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 2 2 2 2 2 2
## [36] 2 2 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1
## [71] 2 1 1 1 1 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 1 1 1
## [106] 1 2 2 2 2 2 2 2 2 1 1 2 2 2 2 1 2 2 1 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 2
## [141] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 2 2 2 2 2 1 1 2 2 2 2 2 2 2 2 2 2
## [176] 2 2 2 2 2 2 2 2 2 2 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2
## [211] 1 1 1 1 2 2 2 2 2 1 2 1 1 2 2 2 2 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [246] 2 2 2 1 1 2 2 2 2 2 2 2 2 2 2 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [281] 2 2 2 2 1 1 2 1 1 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [316] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [351] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [386] 2 2 2 2 2 2 2cutree(hc.average, 2)## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 2 2 2 2 2 2 1 1 1 1 1 2 2 2 1 1 1
## [36] 1 1 1 1 1 1 1 1 1 1 2 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1
## [71] 2 1 1 1 1 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1
## [106] 1 1 2 2 2 2 2 1 2 1 1 2 2 2 2 1 2 1 1 1 1 1 2 2 2 2 1 1 1 1 1 1 1 1 2
## [141] 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 1 2 2 2 2 1 2 1
## [176] 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 2 1 2 1 1
## [211] 1 1 1 1 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 1 2 2 2 2 2
## [246] 2 2 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 1 2 1 2 2 1 1
## [281] 2 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 1 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1
## [316] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [351] 2 2 2 2 2 2 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [386] 2 2 2 2 2 2 2cutree(hc.single, 2)## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [36] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [71] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [106] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [141] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [176] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [211] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [246] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [281] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [316] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [351] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1
## [386] 1 1 1 1 1 1 1cutree(hc.single, 4)## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [36] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [71] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [106] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [141] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [176] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [211] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [246] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [281] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [316] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [351] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 4 1 1 1
## [386] 1 1 1 1 1 1 1xsc=scale(auto)
plot(hclust(dist(xsc), method="complete"), main="Hierarchical Clustering with Scaled Features")
Из всех вариантов остановимся на 2 кластерах, полученных методами дальнего соседа и средней связи.