군집분석 Cluster Analysis
1. 계층적 군집분석
1. 개요
각 객체의 유사성을 측정하여 유사성이 높은 대상 집단을 분류하고, 군집 간의 상이성을 규명하는 분석 방법.
특성에 따라 관측치들을 여러 개의 배타적인 집단으로 나눔.
범주(그룹)에 대한 사전 정보가 없다.
* 그룹의 갯수나 특성에 대한 사전 정보가 있는 경우 분류분석(Classification) 사용
* 요인분석은 유사한 “변수”를 함께 묶어 주고, 군집분석은 “데이터”를 묶어 주는 차이점
군집의 갯수나 구조에 대한 가정 없이 각 데이터 간의 거리를 기준으로 군집화 유도
2. 유사성 측도
1) 연속형 변수
* 유클리드 거리 (Euclidean distance) : 주로 표준화된 자료에 사용
* 맨하탄 거리 (Manhattan)
* 민코우스키 거리 (Minkowski)
* 마할라노비스 거리 (Mahalanobis) : 공분산 행렬을 고려
2) 범주형 변수
* 자카드 거리(Jaccard)
3. 계층적 군집분석 (Hierachical, Agglomerative Clustering)
데이터 간의 유사성을 계산해 가장 가까운 객체들부터 차례로 군집화
한 번 군집에 속하면 이동 불가능
dendrogram을 사용해 군집 형성 과정 파악 가능
군집 간의 거리 차이에 큰 변화를 보이는 경우를 고려해 군집 갯수 결정
1) 최단 연결법 (Single linkage)
distance matrix 에서 가장 작은 값을 취하여 군집화
군집과 군집/데이터 간의 거리 중 최단거리(min) 값을 거리로 산정
길게 늘어진 형태의 군집이 형성될 수 있음
2) 최장 연결법 (Complete linkage)
군집과 군집/데이터 간의 거리 중 최장거리(max) 값을 거리로 산정
둥그런 형태의 군집 형성
단점 : 이상치에 민감
3) 평균 연결법 (Average linkage)
군집과 군집/데이터 간의 거리의 평균거리(mean) 값을 거리로 산정
4) Ward 연결법 (Ward’s method)
군집 간 정보의 손실을 최소화하는 군집화
군집 내 편차들의 제곱합을 고려하여 군집 내 거리(within cluster distance, ESS)를 최소화.
비슷한 크기의 군집을 생성하는 경향
# USArrests data
# Arrests per 100,000 residents for assault, murder, and rape in each of the 50 US states in 1973.
# percent of the population living in urban areas.
df = USArrests
dist(df, method = "manhattan") # "euclidean", "maximum", "manhattan", "canberra", "binary" or "minkowski"
## Alabama Alaska Arizona Arkansas California Colorado
## Alaska 63.5
## Arizona 94.9 78.4
## Arkansas 60.1 101.2 146.2
## California 96.6 60.9 39.5 148.3
## Colorado 74.8 96.9 99.9 62.1 88.0
## Connecticut 165.0 222.1 211.7 120.9 215.2 127.2
## Delaware 28.7 81.8 81.4 76.6 84.9 64.9
## Florida 133.9 122.0 49.2 194.0 85.1 147.3
## Georgia 35.8 90.1 117.5 45.9 119.2 47.4
## Hawaii 223.9 281.0 264.6 181.2 262.1 184.1
## Idaho 137.6 186.7 222.3 85.5 225.8 137.8
## Illinois 43.6 69.9 57.3 98.1 53.0 67.2
## Indiana 136.2 193.3 206.9 95.1 210.4 122.4
## Iowa 201.9 257.0 286.6 155.8 290.1 202.1
## Kansas 139.4 196.5 208.1 95.3 211.6 123.6
## Kentucky 141.4 186.5 229.3 87.1 231.0 145.2
## Louisiana 24.2 59.7 75.1 84.3 76.8 81.0
## Maine 184.5 227.6 269.2 126.4 272.7 184.7
## Maryland 81.5 74.0 25.4 137.8 63.1 121.3
## Massachusetts 127.7 184.8 168.4 83.6 161.9 87.9
## Michigan 50.0 45.5 53.1 107.9 46.6 62.8
## Minnesota 188.8 245.9 257.5 144.7 261.0 173.0
## Mississippi 44.0 41.5 92.9 84.7 94.6 118.8
## Missouri 81.2 124.3 129.7 40.9 131.4 45.6
## Montana 144.0 191.1 228.7 89.9 232.2 144.2
## Nebraska 151.6 208.7 228.3 107.5 231.8 143.8
## Nevada 64.8 47.7 62.1 122.9 42.6 62.6
## New Hampshire 203.8 256.9 288.5 155.7 292.0 204.0
## New Jersey 116.2 173.3 156.9 72.1 142.4 76.4
## New Mexico 73.7 57.8 23.4 130.2 40.9 99.1
## New York 53.0 66.5 53.9 108.9 43.6 73.8
## North Carolina 119.3 108.4 97.8 159.6 135.5 193.7
## North Dakota 231.3 268.4 316.0 171.2 319.5 231.5
## Ohio 139.1 195.8 189.4 98.4 192.9 104.9
## Oklahoma 102.8 159.9 167.5 59.7 171.0 83.0
## Oregon 102.4 143.3 152.9 61.7 156.4 68.4
## Pennsylvania 157.2 214.3 213.9 113.1 217.4 129.4
## Rhode Island 113.7 170.8 154.4 69.6 143.9 73.9
## South Carolina 55.5 42.4 61.8 99.6 69.5 127.7
## South Dakota 180.8 217.9 265.5 120.7 269.0 181.0
## Tennessee 54.7 106.8 136.2 22.8 137.9 52.1
## Texas 61.8 115.7 103.1 50.9 104.8 23.0
## Utah 149.7 203.4 187.0 109.0 190.5 106.5
## Vermont 235.0 272.1 319.7 174.9 323.2 235.2
## Virginia 90.2 147.3 165.7 48.5 168.4 81.6
## Washington 120.2 167.3 164.9 79.5 168.4 80.4
## West Virginia 193.4 230.5 278.1 133.3 281.6 193.6
## Wisconsin 212.0 269.1 280.7 167.9 284.2 196.2
## Wyoming 89.0 146.1 169.7 44.9 173.2 85.2
mahalanobis(df, colMeans(df), cov(df)) # mahalanobis(x, center : 각 변수의 평균, cov : 공분산행렬)
## Alabama Alaska Arizona Arkansas California
## 2.3361665 15.1680638 5.7240749 1.4744001 6.5198834
## Colorado Connecticut Delaware Florida Georgia
## 5.1675431 3.1199978 5.9206962 4.5555037 9.5556141
## Hawaii Idaho Illinois Indiana Iowa
## 7.3838682 2.6996606 2.5572838 1.2856909 2.0936916
## Kansas Kentucky Louisiana Maine Maryland
## 0.5659241 3.6753352 4.5743258 3.0391048 3.2438436
## Massachusetts Michigan Minnesota Mississippi Missouri
## 3.4493654 2.2467698 1.6179094 7.8368676 0.9402315
## Montana Nebraska Nevada New Hampshire New Jersey
## 1.0947097 0.7563206 8.1390645 2.2557657 4.0342537
## New Mexico New York North Carolina North Dakota Ohio
## 2.3134509 2.9272068 12.6102433 4.5070960 1.8374468
## Oklahoma Oregon Pennsylvania Rhode Island South Carolina
## 0.1216993 3.0387808 1.8044534 9.7843187 4.7363035
## South Dakota Tennessee Texas Utah Vermont
## 2.7168038 3.6333816 3.8972486 2.5391464 7.1116744
## Virginia Washington West Virginia Wisconsin Wyoming
## 0.2960040 2.2992245 3.8954199 2.3234779 0.5746895
# linkage methods
x = data.frame(c(1,3,6,12,20))
rownames(x) = c("a","b","c","d","e")
d = dist(x)
par(mfcol=c(2,2))
hc = hclust(d, method = "single") # single linkage (최단연결법)
plot(hc)
hc2 = hclust(d, method = "complete") # complete linkage (최장연결법)
plot(hc2)
hc3 = hclust(d, method = "average") # average linkage (평균연결법)
plot(hc3)
hc4 = hclust(d, method = "ward.D") # ward method
plot(hc4)
par(mfcol=c(1,1))
hc$height
## [1] 2 3 6 8
# USArrests
distUS = dist(scale(df))
hc = hclust(distUS, method = "single")
plot(hc)
hc2 = hclust(distUS, method = "complete")
plot(hc2)
hc3 = hclust(distUS, method = "average")
plot(hc3)
hc4 = hclust(distUS, method = "ward.D")
plot(hc4)
# Cluster Plot
# 위 결과중 '평균연결법' 사용.
h3result = cutree(hc3, k=5) # k : 그룹의 갯수
plot(df$Murder, df$Assault, col=h3result, pch=h3result)
text(df$Murder, df$Assault, labels = rownames(df), col=h3result, pos = 1)
plot(df$Murder, df$UrbanPop, col=h3result, pch=h3result)
text(df$Murder, df$UrbanPop, labels = rownames(df), col=h3result, pos = 1)
df$cluster = h3result
par(mfcol=c(2,2))
for (i in 1:4) {
boxplot(df[,i] ~ h3result, main = names(df)[i])
}
par(mfcol=c(1,1))
library(psych)
describe(df)
## vars n mean sd median trimmed mad min max range skew
## Murder 1 50 7.79 4.36 7.25 7.53 5.41 0.8 17.4 16.6 0.37
## Assault 2 50 170.76 83.34 159.00 168.47 110.45 45.0 337.0 292.0 0.22
## UrbanPop 3 50 65.54 14.47 66.00 65.88 17.79 32.0 91.0 59.0 -0.21
## Rape 4 50 21.23 9.37 20.10 20.36 8.60 7.3 46.0 38.7 0.75
## cluster 5 50 3.44 1.20 4.00 3.55 1.48 1.0 5.0 4.0 -0.87
## kurtosis se
## Murder -0.95 0.62
## Assault -1.15 11.79
## UrbanPop -0.87 2.05
## Rape 0.08 1.32
## cluster -0.12 0.17
describeBy(df, group = h3result)
## $`1`
## vars n mean sd median trimmed mad min max range skew
## Murder 1 7 14.67 1.69 14.4 14.67 1.78 13.0 17.4 4.4 0.37
## Assault 2 7 251.29 48.38 249.0 251.29 44.48 188.0 337.0 149.0 0.42
## UrbanPop 3 7 54.29 8.54 58.0 54.29 11.86 44.0 66.0 22.0 -0.04
## Rape 4 7 21.69 4.03 22.2 21.69 5.34 16.1 26.9 10.8 -0.13
## cluster 5 7 1.00 0.00 1.0 1.00 0.00 1.0 1.0 0.0 NaN
## kurtosis se
## Murder -1.66 0.64
## Assault -1.09 18.28
## UrbanPop -1.90 3.23
## Rape -1.68 1.52
## cluster NaN 0.00
##
## $`2`
## vars n mean sd median trimmed mad min max range skew
## Murder 1 1 10.0 NA 10.0 10.0 0 10.0 10.0 0 NA
## Assault 2 1 263.0 NA 263.0 263.0 0 263.0 263.0 0 NA
## UrbanPop 3 1 48.0 NA 48.0 48.0 0 48.0 48.0 0 NA
## Rape 4 1 44.5 NA 44.5 44.5 0 44.5 44.5 0 NA
## cluster 5 1 2.0 NA 2.0 2.0 0 2.0 2.0 0 NA
## kurtosis se
## Murder NA NA
## Assault NA NA
## UrbanPop NA NA
## Rape NA NA
## cluster NA NA
##
## $`3`
## vars n mean sd median trimmed mad min max range skew
## Murder 1 12 10.88 2.16 11.20 10.73 1.85 7.9 15.4 7.5 0.34
## Assault 2 12 256.92 45.46 254.50 257.00 51.89 178.0 335.0 157.0 -0.16
## UrbanPop 3 12 78.33 7.02 80.00 78.20 6.67 67.0 91.0 24.0 -0.02
## Rape 4 12 32.25 6.73 31.45 31.70 6.67 24.0 46.0 22.0 0.61
## cluster 5 12 3.00 0.00 3.00 3.00 0.00 3.0 3.0 0.0 NaN
## kurtosis se
## Murder -0.72 0.62
## Assault -1.06 13.12
## UrbanPop -1.07 2.03
## Rape -0.92 1.94
## cluster NaN 0.00
##
## $`4`
## vars n mean sd median trimmed mad min max range skew
## Murder 1 23 5.53 2.09 5.9 5.45 2.22 2.6 9.7 7.1 0.19
## Assault 2 23 129.43 43.35 120.0 128.95 43.00 46.0 238.0 192.0 0.23
## UrbanPop 3 23 68.91 11.23 67.0 68.79 10.38 50.0 89.0 39.0 0.09
## Rape 4 23 17.79 4.80 16.5 17.61 4.45 8.3 29.3 21.0 0.30
## cluster 5 23 4.00 0.00 4.0 4.00 0.00 4.0 4.0 0.0 NaN
## kurtosis se
## Murder -1.12 0.44
## Assault 0.14 9.04
## UrbanPop -1.00 2.34
## Rape 0.00 1.00
## cluster NaN 0.00
##
## $`5`
## vars n mean sd median trimmed mad min max range skew
## Murder 1 7 2.70 1.58 2.2 2.70 0.15 0.8 5.7 4.9 0.75
## Assault 2 7 65.14 17.58 57.0 65.14 17.79 45.0 86.0 41.0 0.11
## UrbanPop 3 7 46.29 9.09 45.0 46.29 8.90 32.0 57.0 25.0 -0.20
## Rape 4 7 9.89 1.99 9.5 9.89 2.52 7.3 12.8 5.5 0.06
## cluster 5 7 5.00 0.00 5.0 5.00 0.00 5.0 5.0 0.0 NaN
## kurtosis se
## Murder -0.82 0.60
## Assault -2.08 6.65
## UrbanPop -1.59 3.43
## Rape -1.70 0.75
## cluster NaN 0.00
##
## attr(,"call")
## by.data.frame(data = x, INDICES = group, FUN = describe, type = type)
2. 비계층적 군집분석 : K-means clustering
n개의 객체를 g개의 군집으로 나눌 수 있는 모든 가능한 방법을 점검해 최적화된 군집을 형성
K-means clustering
원하는 군집의 갯수와 초기값 seed를 정해 seed 중심으로 군집을 형성한다.
각 객체는 거리가 가장 가까운 seed가 있는 군집으로 분류된다.
각 군집의 seed 값을 다시 계산한다.
모든 객체를 갱신된 seed와 비교하여 필요시 적절한 군집으로 이동시킨다.
더이상 객체의 이동이 없을 때까지 위 과정을 반복한다.
1) 장점
주어진 데이터의 내부 구조에 대한 사전정보 없이 의미있는 자료 구조를 찾을 수 있다.
다양한 형태의 데이터에 적용 가능
분석 방법 적용도 쉽다.
2) 단점
가중치와 거리에 대한 정의가 어렵다.
초기 군집수 결정이 어렵다.
사전에 주어진 목적이 없기 때문에 결과 해석이 어려울 수 있다.
kdata = iris
kdata$Species = NULL
head(kdata)
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 5.1 3.5 1.4 0.2
## 2 4.9 3.0 1.4 0.2
## 3 4.7 3.2 1.3 0.2
## 4 4.6 3.1 1.5 0.2
## 5 5.0 3.6 1.4 0.2
## 6 5.4 3.9 1.7 0.4
# 3개의 군집으로 나누는 경우
kmc = kmeans(kdata, 3)
kmc
## K-means clustering with 3 clusters of sizes 38, 62, 50
##
## Cluster means:
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 6.850000 3.073684 5.742105 2.071053
## 2 5.901613 2.748387 4.393548 1.433871
## 3 5.006000 3.428000 1.462000 0.246000
##
## Clustering vector:
## [1] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## [36] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [71] 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 2 2 2 1 2 1 1 1
## [106] 1 2 1 1 1 1 1 1 2 2 1 1 1 1 2 1 2 1 2 1 1 2 2 1 1 1 1 1 2 1 1 1 1 2 1
## [141] 1 1 2 1 1 1 2 1 1 2
##
## Within cluster sum of squares by cluster:
## [1] 23.87947 39.82097 15.15100
## (between_SS / total_SS = 88.4 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss"
## [5] "tot.withinss" "betweenss" "size" "iter"
## [9] "ifault"
table(kmc$cluster, iris$Species)
##
## setosa versicolor virginica
## 1 0 2 36
## 2 0 48 14
## 3 50 0 0
plot(kdata$Sepal.Length, kdata$Sepal.Width, col = kmc$cluster, xlab = "Sepal.Length", ylab = "Sepal.Width")
points(kmc$centers[, c(1:2)], col = "blue", pch = 4, cex = 2)
plot(kdata$Petal.Length, kdata$Petal.Width, col = kmc$cluster, xlab = "Petal.Length", ylab = "Petal.Width")
points(kmc$centers[, c(3:4)], col = "blue", pch = 4, cex = 2)
# USArrests
kmUSA = kmeans(df, 3)
kmUSA
## K-means clustering with 3 clusters of sizes 20, 14, 16
##
## Cluster means:
## Murder Assault UrbanPop Rape cluster
## 1 4.270000 87.5500 59.75000 14.39000 4.350000
## 2 8.214286 173.2857 70.64286 22.84286 3.357143
## 3 11.812500 272.5625 68.31250 28.37500 2.375000
##
## Clustering vector:
## Alabama Alaska Arizona Arkansas California
## 3 3 3 2 3
## Colorado Connecticut Delaware Florida Georgia
## 2 1 3 3 2
## Hawaii Idaho Illinois Indiana Iowa
## 1 1 3 1 1
## Kansas Kentucky Louisiana Maine Maryland
## 1 1 3 1 3
## Massachusetts Michigan Minnesota Mississippi Missouri
## 2 3 1 3 2
## Montana Nebraska Nevada New Hampshire New Jersey
## 1 1 3 1 2
## New Mexico New York North Carolina North Dakota Ohio
## 3 3 3 1 1
## Oklahoma Oregon Pennsylvania Rhode Island South Carolina
## 2 2 1 2 3
## South Dakota Tennessee Texas Utah Vermont
## 1 2 2 1 1
## Virginia Washington West Virginia Wisconsin Wyoming
## 2 2 1 1 2
##
## Within cluster sum of squares by cluster:
## [1] 19268.310 9151.857 19579.613
## (between_SS / total_SS = 86.5 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss"
## [5] "tot.withinss" "betweenss" "size" "iter"
## [9] "ifault"
plot(df$Murder, df$Assault, type = "n", xlab = "Murder", ylab = "Assault", main = "3 Cluster (K-means)")
text(df$Murder, df$Assault, labels = rownames(df), cex = 0.8, col = kmUSA$cluster)
points(kmUSA$centers[, c(1,2)], col = "blue", pch = 4, cex = 2)
# 클러스터 갯수 결정시 지표
# within_SS : Within cluster sum of squares by cluster: 군집내 거리의 제곱합. 작을수록 좋다.
# between_SS 군집간 거리의 제곱합. 클수록 좋다.
kmUSA$withinss
## [1] 19268.310 9151.857 19579.613
wss = c()
for (k in 1:20) {
km = kmeans(df, k)
wss[k] = sum(km$withinss)
}
wss
## [1] 355878.142 96428.739 47999.780 37688.090 28266.149 22682.508
## [7] 16588.922 13283.929 12845.085 11378.673 11630.426 13328.633
## [13] 7565.675 6796.296 6560.102 6909.992 5560.443 5908.651
## [19] 4692.008 4165.023
plot(1:20, wss, type = 'l')
points(wss)
3. 모형 기반 군집분석 (Model-based clustering)
확률분포에 대한 정보가 있을 경우 이를 활용하여 군집분석
k번째 군집에 속한 관측치 x 가 다변량 정규분포인 확률밀도함수 f 를 가진다고 가정
각 객체각 각 군집에 속할 사후 확률을 계산하여 가장 확률이 높은 군집으로 할당
모형선택방법 중 BIC (Bayesian information criterion)를 사용하여 이 값이 최대가 되는 모형을 선택
lamda * I : I = identical matrix. 각 dimension의 분산이 같고 각 변수가 서로 독립적이다.
모든 모형에 대해 계산을 한 후 BIC를 기준으로 모형을 선택한다.
library(mclust)
mc = Mclust(USArrests)
summary(mc)
## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust VEI (diagonal, equal shape) model with 3 components:
##
## log.likelihood n df BIC ICL
## -757.5611 50 20 -1593.363 -1597.743
##
## Clustering table:
## 1 2 3
## 20 20 10
# 최종 선택 모형 : Mclust VEI (diagonal, equal shape) model with 3 components:
plot(mc)
# BIC plot을 보고 적당한 클러스터 갯수를 비교하여 선택한다.
# classification plot : 타원 = 공분산의 모양.
# uncertainry : posterior probability (사후 확률) 기준으로 어느 그룹에 속할지 확실하지 않은 관측치.
mc$classification
## Alabama Alaska Arizona Arkansas California
## 1 1 1 2 1
## Colorado Connecticut Delaware Florida Georgia
## 1 2 2 1 1
## Hawaii Idaho Illinois Indiana Iowa
## 2 3 1 2 3
## Kansas Kentucky Louisiana Maine Maryland
## 2 2 1 3 1
## Massachusetts Michigan Minnesota Mississippi Missouri
## 2 1 3 1 1
## Montana Nebraska Nevada New Hampshire New Jersey
## 2 2 1 3 2
## New Mexico New York North Carolina North Dakota Ohio
## 1 1 1 3 2
## Oklahoma Oregon Pennsylvania Rhode Island South Carolina
## 2 2 2 2 1
## South Dakota Tennessee Texas Utah Vermont
## 3 1 1 2 3
## Virginia Washington West Virginia Wisconsin Wyoming
## 2 2 3 3 2
par(mfcol=c(2,2))
for (i in 1:4) {
boxplot(USArrests[,1] ~ mc$classification, main = names(USArrests)[i])
}
par(mfcol=c(1,1))
# compare groups 1
USArrests_s = as.data.frame(scale(USArrests))
par(mfcol=c(1,3))
for (i in 1:3) {
boxplot(USArrests_s[mc$classification == i,], ylim = c(-3,3), main = i)
}
par(mfcol=c(1,1))
# compare groups 2
library(psych)
result = describeBy(USArrests[,1], group = mc$classification, mat = T) # murder
result
## item group1 vars n mean sd median trimmed mad min max
## X11 1 1 1 20 12.165 2.684904 12.15 12.11250 2.89107 7.9 17.4
## X12 2 2 1 20 5.965 1.877646 6.00 5.89375 2.00151 3.2 9.7
## X13 3 3 1 10 2.680 1.293402 2.40 2.53750 0.44478 0.8 5.7
## range skew kurtosis se
## X11 9.5 0.1552937 -1.0023747 0.6003628
## X12 6.5 0.1872682 -1.0161795 0.4198543
## X13 4.9 1.0025078 0.4867888 0.4090096
library(ggplot2)
ggplot(result, aes(group1, mean, fill = group1)) +
geom_bar(position = position_dodge(), stat = "identity") +
geom_errorbar(aes(ymin = mean-2*se, ymax = mean+2*se), width = 0.1)
# confidence inteval이 서로 겹치지 않기 때문에 그룹간 차이가 명확하게 있다.
# 세 그룹간 비교 --> anova
model = aov(USArrests$Murder ~ factor(mc$classification))
summary(model)
## Df Sum Sq Mean Sq F value Pr(>F)
## factor(mc$classification) 2 710.5 355.3 76.24 1.76e-15 ***
## Residuals 47 219.0 4.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 사후 검정 : Tukey Test
TukeyHSD(model)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = USArrests$Murder ~ factor(mc$classification))
##
## $`factor(mc$classification)`
## diff lwr upr p adj
## 2-1 -6.200 -7.852027 -4.547973 0.0000000
## 3-1 -9.485 -11.508312 -7.461688 0.0000000
## 3-2 -3.285 -5.308312 -1.261688 0.0007996
# 각 그룹별 비교시 p-value가 0.05보다 작기 때문에 그룹간 확실한 차이가 있다.
# cluster 갯수 4개로 지정.
mc4 = Mclust(USArrests, 4)
summary(mc4)
## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust EEI (diagonal, equal volume and shape) model with 4 components:
##
## log.likelihood n df BIC ICL
## -755.3597 50 23 -1600.696 -1606.279
##
## Clustering table:
## 1 2 3 4
## 12 7 20 11
plot(mc4)
# Practice 4
# Jet
# 22개 미국 전투기
# - FFD: 처음 비행 날짜
# - CAR: 비행기가 항공모함에 착륙 가능여부
# - SPR: 단위무게 당 출력에 비례하는 특정한 출력
# - RGF: 비행범위 요인
# - PLF: 비행기의 총 무게의 일부분으로서의 탑재량
# - SLF: 일관된 무게 요인
# 1. 계층적군집분석
# A. FFD와 CAR를 제외한 변수를 표준화 한 후 계층적 군집화를 시행하고 덴드로그램을 그리시오.
jet = read.csv("data/jet.csv", header = T)
rownames(jet) = jet$X
jet = jet[, -c(1,2,7)] # FFD, CAR 제외.
jet
## SPR RGF PLF SLF
## FH-1 1.468 3.30 0.166 0.10
## FJ-1 1.605 3.64 0.154 0.10
## F-86A 2.168 4.87 0.177 2.90
## F9F-2 2.054 4.72 0.275 1.10
## F-94A 2.467 4.11 0.298 1.00
## F3D-1 1.294 3.75 0.150 0.90
## F-89A 2.183 3.97 0.000 2.40
## XF10F-1 2.426 4.65 0.117 1.80
## F9F-6 2.607 3.84 0.155 2.30
## F-100A 4.567 4.92 0.138 3.20
## F4D-1 4.588 3.82 0.249 3.50
## F11F-1 3.618 4.32 0.143 2.80
## F-101A 5.855 4.53 0.172 2.50
## F3H-2 2.898 4.48 0.178 3.00
## F-102A 3.880 5.39 0.101 3.00
## F-8A 0.455 4.99 0.008 2.64
## F-104B 8.088 4.50 0.251 2.70
## F-105B 6.502 5.20 0.366 2.90
## YF-107A 6.081 5.65 0.106 2.90
## F-106A 7.105 5.40 0.089 3.20
## F-4B 8.548 4.20 0.222 2.90
## F-111A 6.321 6.45 0.187 2.00
jet_s = scale(jet)
distJet = dist(jet_s)
par(mfrow=c(2,2))
hc1 = hclust(distJet, method = "single")
plot(hc1, main = "Single")
hc2 = hclust(distJet, method = "complete")
plot(hc2, main = "Complete")
hc3 = hclust(distJet, method = "average")
plot(hc3, main = "Average")
hc4 = hclust(distJet, method = "ward.D")
plot(hc4, main = "Ward.D")
par(mfcol=c(1,1))
# B. A의 결과를 사용해 두 개의 집단으로 관측치를 분류하고 각 집단의 특징을 원변수 관점에서 비교하시오.
# 2개의 군집으로 나누기 위해서는 complete linkage 또는 ward method를 선택하는 것이 좋다.
# (1) complete linkage 사용하는 경우.
result_hc = cutree(hc2, k=2)
par(mfrow=c(2,3))
plot(jet$SPR, jet$RGF, col=result_hc, pch=result_hc, xlab = "SPR", ylab = "RGF")
text(jet$SPR, jet$RGF, labels = rownames(jet), col=result_hc, pos = 1)
abline(v = 3.75)
plot(jet$SPR, jet$PLF, col=result_hc, pch=result_hc, xlab = "SPR", ylab = "PLF")
text(jet$SPR, jet$PLF, labels = rownames(jet), col=result_hc, pos = 1)
abline(v = 3.75)
plot(jet$SPR, jet$SLF, col=result_hc, pch=result_hc, xlab = "SPR", ylab = "SLF")
text(jet$SPR, jet$SLF, labels = rownames(jet), col=result_hc, pos = 1)
abline(v = 3.75)
plot(jet$RGF, jet$PLF, col=result_hc, pch=result_hc, xlab = "RGF", ylab = "PLF")
text(jet$RGF, jet$PLF, labels = rownames(jet), col=result_hc, pos = 1)
plot(jet$RGF, jet$SLF, col=result_hc, pch=result_hc, xlab = "RGF", ylab = "SLF")
text(jet$RGF, jet$SLF, labels = rownames(jet), col=result_hc, pos = 1)
plot(jet$PLF, jet$SLF, col=result_hc, pch=result_hc, xlab = "PLF", ylab = "SLF")
text(jet$PLF, jet$SLF, labels = rownames(jet), col=result_hc, pos = 1)
par(mfcol=c(1,1))
# complete linkage 방법을 사용하여 두 개의 군집으로 나눌 경우에는 SPR 변수의 영향을 가장 크게 받았음을 알 수 있다.
# SPR 값이 약 3.7 보다 작은 전투기는 cluster 1, 큰 전투기는 cluster 2로 군집이 형성되었다.
table(result_hc)
## result_hc
## 1 2
## 12 10
# (2) ward method 사용하는 경우
result_hc2 = cutree(hc4, k=2)
par(mfrow=c(2,3))
plot(jet$SPR, jet$RGF, col=result_hc2, pch=result_hc2, xlab = "SPR", ylab = "RGF")
text(jet$SPR, jet$RGF, labels = rownames(jet), col=result_hc2, pos = 1)
plot(jet$SPR, jet$PLF, col=result_hc2, pch=result_hc2, xlab = "SPR", ylab = "PLF")
text(jet$SPR, jet$PLF, labels = rownames(jet), col=result_hc2, pos = 1)
plot(jet$RGF, jet$PLF, col=result_hc2, pch=result_hc2, xlab = "RGF", ylab = "PLF")
text(jet$RGF, jet$PLF, labels = rownames(jet), col=result_hc2, pos = 1)
plot(jet$SLF, jet$SPR, col=result_hc2, pch=result_hc2, xlab = "SLF", ylab = "SPR")
text(jet$SLF, jet$SPR, labels = rownames(jet), col=result_hc2, pos = 1)
abline(v = 1.5)
plot(jet$SLF, jet$RGF, col=result_hc2, pch=result_hc2, xlab = "SLF", ylab = "RGF")
text(jet$SLF, jet$RGF, labels = rownames(jet), col=result_hc2, pos = 1)
abline(v = 1.5)
plot(jet$SLF, jet$PLF, col=result_hc2, pch=result_hc2, xlab = "SLF", ylab = "PLF")
text(jet$SLF, jet$PLF, labels = rownames(jet), col=result_hc2, pos = 1)
abline(v = 1.5)
par(mfcol=c(1,1))
# ward method를 사용하여 두 개의 군집으로 나눌 경우에는 SLF 변수의 분포가 가장 큰 영향을 끼쳤음을 알 수 있다.
# SLF 값이 약 1.5 보다 작은 전투기는 cluster 1, 큰 전투기는 cluster 2로 군집이 형성되었다.
table(result_hc2)
## result_hc2
## 1 2
## 5 17
# C. 두 집단을 주성분을 이용해 2차원 산점도로 표현하시오.
# 즉, 제1 주성분과 제2 주성분을 사용한 산점도에서 두 개의 집단을 서로 다른 마크와 색으로 표현하시오.
# 주성분분석
pca <- prcomp(jet, scale = T)
summary(pca)
## Importance of components:
## PC1 PC2 PC3 PC4
## Standard deviation 1.4064 1.0779 0.7422 0.55609
## Proportion of Variance 0.4945 0.2905 0.1377 0.07731
## Cumulative Proportion 0.4945 0.7850 0.9227 1.00000
# 2개의 주성분으로 전체 변동의 78%를 설명할 수 있다.
library(ggfortify)
jet_c = jet
jet_c$cluster1 = factor(result_hc) # single linkage
jet_c$cluster2 = factor(result_hc2) # ward method
jet_c
## SPR RGF PLF SLF cluster1 cluster2
## FH-1 1.468 3.30 0.166 0.10 1 1
## FJ-1 1.605 3.64 0.154 0.10 1 1
## F-86A 2.168 4.87 0.177 2.90 1 2
## F9F-2 2.054 4.72 0.275 1.10 1 1
## F-94A 2.467 4.11 0.298 1.00 1 1
## F3D-1 1.294 3.75 0.150 0.90 1 1
## F-89A 2.183 3.97 0.000 2.40 1 2
## XF10F-1 2.426 4.65 0.117 1.80 1 2
## F9F-6 2.607 3.84 0.155 2.30 1 2
## F-100A 4.567 4.92 0.138 3.20 2 2
## F4D-1 4.588 3.82 0.249 3.50 2 2
## F11F-1 3.618 4.32 0.143 2.80 1 2
## F-101A 5.855 4.53 0.172 2.50 2 2
## F3H-2 2.898 4.48 0.178 3.00 1 2
## F-102A 3.880 5.39 0.101 3.00 2 2
## F-8A 0.455 4.99 0.008 2.64 1 2
## F-104B 8.088 4.50 0.251 2.70 2 2
## F-105B 6.502 5.20 0.366 2.90 2 2
## YF-107A 6.081 5.65 0.106 2.90 2 2
## F-106A 7.105 5.40 0.089 3.20 2 2
## F-4B 8.548 4.20 0.222 2.90 2 2
## F-111A 6.321 6.45 0.187 2.00 2 2
# single linkage
autoplot(pca, data = jet_c, colour = 'cluster1', shape = F, label = T, loadings = T,
label.size = 7, loadings.label.size = 5, main = "PCA (single linkage cluster)",
loadings.label = T, loadings.colour = "blue", loadings.label.colour = "blue")
# ward method
autoplot(pca, data = jet_c, colour = 'cluster2', shape = F, label = T, loadings = T,
label.size = 7, loadings.label.size = 5, main = "PCA (ward method cluster)",
loadings.label = T, loadings.colour = "blue", loadings.label.colour = "blue")
# 2. 비계층적군집분석
# A. 군집 개수 1~5까지를 사용해 k-means clustering을 시행하고
# 얻은 within-group sum of squares를 저장하고 그래프로 표현하여 적절한 군집 개수를 판단하시오.
wss = c()
for (k in 1:5) {
km = kmeans(jet, k)
wss[k] = sum(km$withinss)
}
wss
## [1] 150.88282 50.79957 31.32728 25.00716 18.71070
plot(1:5, wss, type = 'l')
points(wss)
# 2에서 팔꿈치 생김
# 군집 3~5개 일 경우 서로 wss 차이가 많지 않음
# ---> 2개의 군집으로 결정
# B. K-means clustering을 이용해 2개의 집단으로 군집화하고 그 결과를 1번의 B, C와 같이 탐색하시오.
result_km = kmeans(jet, 2)
result_km
## K-means clustering with 2 clusters of sizes 13, 9
##
## Cluster means:
## SPR RGF PLF SLF
## 1 2.240231 4.310000 0.1478462 1.849231
## 2 6.406111 4.963333 0.1977778 2.866667
##
## Clustering vector:
## FH-1 FJ-1 F-86A F9F-2 F-94A F3D-1 F-89A XF10F-1 F9F-6
## 1 1 1 1 1 1 1 1 1
## F-100A F4D-1 F11F-1 F-101A F3H-2 F-102A F-8A F-104B F-105B
## 2 2 1 2 1 1 1 2 2
## YF-107A F-106A F-4B F-111A
## 2 2 2 2
##
## Within cluster sum of squares by cluster:
## [1] 28.95941 21.84015
## (between_SS / total_SS = 66.3 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss"
## [5] "tot.withinss" "betweenss" "size" "iter"
## [9] "ifault"
par(mfcol=c(1,2))
plot(jet$SPR, jet$RGF, type = "n", xlab = "SPR", ylab = "RGF", main = "2 Cluster (K-means)")
text(jet$SPR, jet$RGF, labels = rownames(jet), cex = 0.8, col = result_km$cluster)
points(result_km$centers[, c(1,2)], col = "blue", pch = 4, cex = 2)
plot(jet$SLF, jet$SPR, type = "n", xlab = "SLF", ylab = "SPR", main = "2 Cluster (K-means)")
text(jet$SLF, jet$SPR, labels = rownames(jet), cex = 0.8, col = result_km$cluster)
points(result_km$centers[, c(4,1)], col = "blue", pch = 4, cex = 2)
par(mfcol=c(1,1))
# 3. 모형기반 군집화를 통해 최적의 군집 개수를 찾고 그 결과를 1번의 B, C와 같이 탐색하시오.
library(mclust)
result_mc = Mclust(jet)
summary(result_mc)
## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust VVI (diagonal, varying volume and shape) model with 3 components:
##
## log.likelihood n df BIC ICL
## -43.76189 22 26 -167.8909 -168.7453
##
## Clustering table:
## 1 2 3
## 3 5 14
# 최종 선택 모형 : 군집 3개. VVI model
plot(result_mc)
result_mc$classification
## FH-1 FJ-1 F-86A F9F-2 F-94A F3D-1 F-89A XF10F-1 F9F-6
## 1 1 3 2 2 1 2 2 2
## F-100A F4D-1 F11F-1 F-101A F3H-2 F-102A F-8A F-104B F-105B
## 3 3 3 3 3 3 3 3 3
## YF-107A F-106A F-4B F-111A
## 3 3 3 3
rownames(jet[result_mc$classification == 1, ])
## [1] "FH-1" "FJ-1" "F3D-1"
rownames(jet[result_mc$classification == 2, ])
## [1] "F9F-2" "F-94A" "F-89A" "XF10F-1" "F9F-6"
rownames(jet[result_mc$classification == 3, ])
## [1] "F-86A" "F-100A" "F4D-1" "F11F-1" "F-101A" "F3H-2" "F-102A"
## [8] "F-8A" "F-104B" "F-105B" "YF-107A" "F-106A" "F-4B" "F-111A"
jet[order(result_mc$uncertainty, decreasing = T), ][1, ]
## SPR RGF PLF SLF
## F-86A 2.168 4.87 0.177 2.9
result_km3 = kmeans(jet, 3)
plot(jet$SPR, jet$RGF, type = "n", xlab = "SPR", ylab = "RGF", main = "3 Cluster (K-means)")
text(jet$SPR, jet$RGF, labels = rownames(jet), cex = 0.8, col = result_km3$cluster)
points(result_km3$centers[, c(1,2)], col = "blue", pch = 4, cex = 2)
