이 문서 및 군집분석 앱을 작성하는 데에 참조한 문헌들은 다음과 같다.

1. R in action(2nd edition)
2. Practical Data Science with R
3. An introduction to statistical learning
4. clValid, an R package for cluster validation

이 문서에 사용된 R 패키지들은 다음과 같다.

requiredPackages=c("flexclust","NbClust","rattle","cluster","fMultivar","ggplot2")
for(p in requiredPackages){
        if(!require(p,character.only=TRUE)) install.packages(p)
}

웹R에서 만든 cluster analysis 샤이니 앱을 실행하기 위해서는 다음 패키지들이 필요하다.

requiredPackages=c("shiny","d3heatmap","RColorBrewer","DT","fpc","clValid","RankAggreg")
for(p in requiredPackages){
        if(!require(p,character.only=TRUE)) install.packages(p)
}

군집분석

군집분석은 주어진 데이타셋 내에 존재하는 몇 개의 군집을 찾아내는 비지도(unsupervised) 기법이다. 군집(cluster)는 다른 그룹에 속한 다른 관찰치들에 비해 서로 보다 유사한 관찰치들의 그룹이라고 할 수 있는데 그 정의가 정밀하지 않기 때문에 수 많은 군집 방법이 존재하게 된다.

군집분석의 예 : http://shiny.rstudio.com/gallery/kmeans-example.html

군집분석이 유용한 경우

군집분석은 생물학, 행동과학, 마케팅 및 의학분야에서 다양하게 사용된다. 예를 들어 정신과 영역에서는 우울증의 subtype 을 발견하기 위해 우울증 환자의 중상과 인구학적인 특징에 관한 데이타로 군집분석을 할 수 있다. 이를 통해 subtype이 발견된다면 질병에 대한 이해를 넓힐수 있고 보다 효과적인 치료도 가능할 것이다. 마케팅 담당자라면 고객의 인구학적인 특징 및 소비 성향의 유사성에 근거해 몇개의 군집으로 나눔으로써 맞춤형 마케팅을 하는데 활용할 수 있을 것이다. 의학연구자들은 유전체데이타를 이용하여 유사한 표현형 및 공통되는 생물학적 경로를 가지는 유전자와 단백질들을 grouping하는데 군집분석을 사용할 수 있다.

군집분석 방법

  1. Hierarchical agglomerative clustering
  • 모든 관찰치는 자신만의 군집에서 시작하여 유사한 데이타 두 개를 하나의 군집으로 묶는데 이를 모든 데이타가 하나의 군집으로 묶일때까지 반복한다.
  • 사용되는 알고리즘 : single linkage, complete linkage, average linkage, controid, Ward의 방법 등이 있다.
  1. Partitioning clustering
  • 먼저 군집의 갯수 K 를 정한 후 데이타를 무작위로 K개의 군으로 배정한 후 다시 계산하여 군집으로 나눈다.

  • 사용되는 알고리즘 : k-means 및 PAM(partitioning around medoids) 등이 있다.

군집분석 단계

  1. 알맞은 속성 선택(Choose appropriate attributes)

첫번째(이자 아마도 가장 중요한) 단계는 데이타를 군집화하는데 중요하다고 판단되는 속성들을 선택하는 것이다. 예를 들어 우울증에 대한 연구라고 하면 다음과 같은 속성들을 평가할 수 있다. 정신과적 증상, 이학적증상, 발병나이, 우울증의 횟수, 지속기간, 빈도, 입원 횟수, 기능적 상태, 사회력 및 직업력, 현재 나이, 성별, 인종, 사회경제적 상태, 결혼상태, 가족력, 과거 치료에 대한 반응 등. 아무리 복잡하고 철저하게 군집분석을 하더라도 잘못 선택한 속성을 극복할 수 없다.(A sophisticated cluster analysis can’t compaensate for a poor choice of variables.)

  1. 데이타 표준화(Scale the data)

분석에 사용되는 변수들의 범위에 차이가 있는 경우 가장 큰 범위를 갖는 변수가 결과에 가장 큰 영향을 미치게 된다. 이런 결과가 바람직하지 않은 경우 데이타를 표준화 할 수 있다. 가장 많이 사용되는 방법은 각 변수를 평균 0, 표준편차 1로 표준화하는 것이다.

df1 <- apply(mydata, 2, function (x) {(x-mean(x))/sd(x)})

R의 scale() 함수를 사용하면 같은 결과를 얻을 수 있다. 웹R에서도 이 방법을 사용한다. 그 외에 각 변수를 최대값으로 나누는 방법도 있고(df2) 각 값에서 평균을 빼고 median absolute deviation으로 나눌 수도 있다(df3).

df2 <- apply(mydata, 2, function (x) {x/max(x)})
df3 <- apply(mydata, 2, function (x) {(x-mean(x))/mad(x)})
  1. 이상치 선별(Screen for outliers)

많은 군집분석 방법은 이상치에 민감하기 때문에 이상치가 있는 경우 군집분석 결과가 왜곡된다. 단변량 이상치의 경우 outlier 패키지의 함수를 사용할 수 있고 다변량 이상치의 경우 mvoutlier 패키지에 있는 함수들을 사용하여 이상치를 선별하고 제거할 수 있다. 또 다른 방법은 이상치에 대하여 강건한(robust) 군집분석 방법을 쓸 수 있는데 PAM(partitioning around medoids)이 대표적인 방법이다.

  1. 거리의 계산 (Calcuate distance)

두 관찰치 간의 거리를 측정하는 방법은 여러가지가 있는데 “euclidean”, “maximum”, “manhattan”, “canberra”, “binary” 또는 “minkowski” 방법을 사용할 수가 있다. R의 dist()함수를 쓰면 위의 방법들을 이용하여 거리를 계산할 수 있으며 웹R 에서도 이 함수를 사용하여 거리를 계산한다. 디폴트 값은 유클리드 거리(“euclidean”)이다. 유클리드 거리는 다음 절에서 다룬다.

  1. 군집 알고리즘 선택

다음 단계로 군집 방법을 선택하여야 한다. 계층적군집(Hierarchical agglomerative clustering)은 150 관찰치 이하의 적은 데이타에 적합하다. 분할군집은 보다 많은 데이타를 다룰 수 있으나 군집의 갯수를 정해주어야 한다. 계층적 군집/분할 군집을 선택한 후 구체적인 방법을 선택하여야 한다. 각각의 방법은 장단점이 있기 때문에 하나 이상의 방법을 사용해보고 어느 방법이 강건한지 평가해 볼 수 있다.

  1. 하나 이상의 군집분석 결과 얻음

  2. 군집의 갯수 결정

군집분석 최종 결과를 얻기 위해 몇 개의 군집이 있는지 결정해야 한다. NbClust패키지의 NbClust()힘수를 사용할 수 있다. 웹R에서는 clValid패키지를 사용하여 군집 갯수 결정에 도움을 주고자 하였다.

  1. 최종결과 획득

  2. 분석 결과의 시각화

최종 결과를 시각화할 때 계층적 분석은 dendrogram으로 나타내고 분할군집은 이변량 cluster plot으로 시각화한다.

  1. 군집분석 결과의 해석

최종 결과를 얻은 후 그 결과를 해석하고 가능하면 이름도 지어야 한다. 한 군집의 관측치가 갖고 있는 공통점은 무엇인가? 다른 군집과 어떤 점이 다른가? 이 단계는 각 군집의 통계량을 요약함으로써 얻어진다. 연속형 변수의 경우 평균 또는 중앙값을 계산하고 범주형 변수가 있는 경우 범주별로 각 군집의 분포를 보아야 한다.

  1. 결과의 확인(validate the result)

다음과 같은 질문을 가져본다. “이 군집 결과가 실제 존재하는가? 이 데이타셋에서만 또는 이 방법에서만 나타나는 것은 아닌가?” 다른 군집 방법 또는 다른 데이타셋이 사용될 경우 같은 군빕이 얻어질 것인가? 결과의 validation을 위해 fpc, clv 또는 clValid 패키지 등을 사용할 수 있다. 웹R에서는 clValid 패키지를 사용한다.

거리의 계산

유클리드 거리는 다음과 같은 공식으로 계산된다.

\(d_{ij}=\sqrt{\sum^{p}_{p=1}(x_{ip}-x_{jp})^2}\)

여기서 i,j는 관측치이며 P는 변수 번호이다. 에를 들어 nutrient데이타를 살펴보자.

data(nutrient,package="flexclust")
head(nutrient,4)
             energy protein fat calcium iron
BEEF BRAISED    340      20  28       9  2.6
HAMBURGER       245      21  17       9  2.7
BEEF ROAST      420      15  39       7  2.0
BEEF STEAK      375      19  32       9  2.6

첫 두 데이타(beef braised와 hamburger)사이의 유클리드 거리는 다음과 같이 구할 수 있다.

\(d=\sqrt{(340-245)^2+(20-21)^2+(28-17)^2+(9-9)^2+(2.6-2.7)^2}=95.64\)

R의 dist()함수는 데이타프레임 또는 행렬의 모든 행 사이의 거리를 계산하여 행렬 형식으로 결과를 반환헤 준다. 다음과 같이 할 수 있다.

d=dist(nutrient)
as.matrix(d)[1:4,1:4]
             BEEF BRAISED HAMBURGER BEEF ROAST BEEF STEAK
BEEF BRAISED      0.00000   95.6400   80.93429   35.24202
HAMBURGER        95.64000    0.0000  176.49218  130.87784
BEEF ROAST       80.93429  176.4922    0.00000   45.76418
BEEF STEAK       35.24202  130.8778   45.76418    0.00000

관측치 사이의 거리가 크다는 것은 관측치가 유사하지 않다는 것이다. 어떤 관측치와 자신과의 거리는 0이다. beef braised와 hamburger 사이의 거리는 손으로 계산한 값과 같다.

계층적군집 방법

모든 관찰치는 자신만의 군집에서 시작하여 유사한 데이타 두 개를 하나의 군집으로 묶는데 이를 모든 데이타가 하나의 군집으로 묶일때까지 반복한다. 알고리즘은 다음과 같다.

1. 모든 관찰치를 군집으로 정의한다.
2. 모든 군집에 대하여 다른 모든 군집과의 거리를 계산한다.
3. 가장 작은 거리를 갖는 두 군집을 합해 하나의 군집으로 만든다. 따라서 군집의 갯수가 하나 감소한다.
4. 2와3을 반복하여 모든 관찰치가 하나의  군집으로 합쳐질 때까지 반복한다.

2단계에서 군집 사이의 거리를 정의하는 것에 따라 계층적 군집 알고리즘이 달라진다. 가장 많이 쓰이는 다섯가지 방법의 정의는 다음과 같다.

군집방법 두 군집 사이의 거리 정의
single linkage 한 군집의 점과 다른 군집의 점 사이의 가장 짧은 거리(shortest distance)
complete linkage 한 군집의 점과 다른 군집의 점 사이의 가장 긴 거리(longest distance)
average linkage 한 군집의 점과 다른 군집의 점 사이의 평균 거리. UPGMA(unweighted pair group mean averaging)이라고도 한다.
centroid 두 군집의 centroids(변수 평균의 벡터) 사이의 거리.관측치가 하나인 경우 centroid는 변수의 값이 된다
Ward 모든 변수들에 대하여 두 군집의 ANOVA sum of square를 더한 값

single linkage clustering은 긴 시가모양의 군집이 만들어지는 경향이 있으며 이러한 현상을 chaining이라고 한다. chaining은 유사하지 않은 관측치들의 중간 관측치들이 유사하기 때문에 하나의 군집으로 합쳐지는 것를 말한다. complete linkage clustering은 거의 비슷한 직경을 갖는 compact cluster를 만드는 경향이 있으며 이상치에 민감한 것으로 알려져 있다. average linkage clustering은 두 가지 방법의 타협점이다. chaining결향이 덜하고 이상치에도 덜 민감하다. 또한 분산이 적은 군집을 만드는 경향이 있다. Ward의 방법은 적은 관찰치를 갖는 군집을 만드는 경향이 있으며 관찰치의 수와 거의 같은 군집을 만드는 경향이 있다. centroid방법은 단순하고 이해하기 쉬운 거리의 정의를 갖는 매력적인 방법으로 다른 방법들에 비해 이상치에 덜 민감하지만 average linkage나 Ward방법에 비해 수행능력이 떨어진다.

R을 이용해 계층적 분석을 하려면 다음과 같이 한다.

hclust(d, method=)

d는 dist()함수에 의해 만들어지는 거리행렬이고 method로는 “ward.D”, “ward.D2”, “single”, “complete”, “average” (= UPGMA), “mcquitty” (= WPGMA), “median” (= WPGMC) 또는 “centroid” (= UPGMC)를 사용할 수 있다.

R을 이용한 계층적 군집분석

data(nutrient,package="flexclust")
rownames(nutrient)=tolower(rownames(nutrient))
nutrient.scaled=scale(nutrient)

d=dist(nutrient.scaled)
fit.average=hclust(d, method="average")
plot(fit.average,hang=-1,cex=.8,main="Average Linkage Clustering")

몇 개의 군집으로 나누어야 하는가?

library(NbClust)
devAskNewPage(ask=TRUE)
nc=NbClust(nutrient.scaled,distance="euclidean",min.nc=2,max.nc=15,
           method="average")

*** : 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:                                                
* 4 proposed 2 as the best number of clusters 
* 4 proposed 3 as the best number of clusters 
* 2 proposed 4 as the best number of clusters 
* 4 proposed 5 as the best number of clusters 
* 1 proposed 9 as the best number of clusters 
* 1 proposed 10 as the best number of clusters 
* 2 proposed 13 as the best number of clusters 
* 1 proposed 14 as the best number of clusters 
* 4 proposed 15 as the best number of clusters 

                   ***** Conclusion *****                            
 
* According to the majority rule, the best number of clusters is  2 
 
 
******************************************************************* 
devAskNewPage(ask=FALSE)

table(nc$Best.n[1,])

 0  1  2  3  4  5  9 10 13 14 15 
 2  1  4  4  2  4  1  1  2  1  4 
par(mfrow=c(1,1))
barplot(table(nc$Best.n[1,]),xlab="Number of Clusters",ylab="Number of Criteria",
        main="Number of Clusters Chosen by 26 criteria")

최종 결과 획득

clusters<-cutree(fit.average,k=5)
table(clusters)
clusters
 1  2  3  4  5 
 7 16  1  2  1 
aggregate(nutrient,by=list(cluster=clusters),median)
  cluster energy protein fat calcium iron
1       1  340.0      19  29       9 2.50
2       2  170.0      20   8      13 1.45
3       3  160.0      26   5      14 5.90
4       4   57.5       9   1      78 5.70
5       5  180.0      22   9     367 2.50
aggregate(as.data.frame(nutrient.scaled),by=list(cluster=clusters),median)
  cluster     energy    protein        fat    calcium        iron
1       1  1.3101024  0.0000000  1.3785620 -0.4480464  0.08110456
2       2 -0.3696099  0.2352002 -0.4869384 -0.3967868 -0.63743114
3       3 -0.4684165  1.6464016 -0.7534384 -0.3839719  2.40779157
4       4 -1.4811842 -2.3520023 -1.1087718  0.4361807  2.27092763
5       5 -0.2708033  0.7056007 -0.3981050  4.1396825  0.08110456

분석결과 시각화

plot(fit.average,hang=-1,cex=.8,
     main="Average Linkage Clustering\n5 Cluster Solution")
rect.hclust(fit.average,k=5)

분할군집

분할군집에서는 먼저 군집의 갯수 K 를 정한 후 데이타를 무작위로 K개의 군으로 배정한 후 다시 계산하여 군집으로 나눈다. k-means clustering과 PAM을 다룬다.

k-means clustering

k-means 알고리즘

분할 군집에서 가장 많이 사용되는 방법은 k-means clustering이다. 알고리즘은 다음과 같다.

1. K개의 centroids를 선택한다.(K개의 행을 무작위로 선택)
2. 각 데이타를 가장 가까운 centroid에 할당한다.
3. 각 군집에 속한 모든 데이타의 평균으로 centroid를 다시 계산한다.(즉, centroid는 p-개의 길이를 갖는 평균벡터로 p는 변수의 수이다.)
4. 각 데이타를 가장 가까운 centroid에 할당한다.
5. 모든 관측치의 재할당이 일어나지 않거나 최대반복횟수(R에서의 dafault값은 10회)에 도달할 때까지 3과 4를 반복한다.

R에서 가장 가까운 centroid에 할당할때 다음 값을 이용하여 계산한다.

\(ss(k)=\sum^n_{i=1}\sum^p_{j=0}(x_{ij}-\bar{x}_{kj})^2\)

여기서 k는 군집이고 \(x_{ij}\)는 i번째 관측치의 j번째 변수이고 \(\bar{x}_{kj}\)는 k번째 군집의 j번째 변수의 평균이고 p는 변수의 갯수이다.

k-means의 장단점

k-means는 계층적 군집분석에 비해 큰 데이타셋에서 사용할 수 있으며 관측치가 군집에 영구히 할당되는 것이 아니라 최종결과를 개선시키는 방향으로 이동한다. 하지만 평균을 사용하기 때문에 연속형변수에만 적용될 수 있으며 이상치에 심하게 영향을 받는다. 또한 non-convex형태의(예를 들어 U모양) 군집이 있는 경우 잘 수행되지 않는다.

R을 이용한 k-means 군집분석

k-means 군집 분석을 할 때 무작위로 K개의 행을 선택하므로 실행할 때마다 결과가 달라진다. set.seed()함수를 쓰면 재현 가능한 결과를 얻을 수 있다. K값을 결정하기 위해 계층적 군집분석에서 사용했던 NbClust()함수를 이용할 수 있다.

data(wine,package="rattle")
head(wine)
  Type Alcohol Malic  Ash Alcalinity Magnesium Phenols Flavanoids
1    1   14.23  1.71 2.43       15.6       127    2.80       3.06
2    1   13.20  1.78 2.14       11.2       100    2.65       2.76
3    1   13.16  2.36 2.67       18.6       101    2.80       3.24
4    1   14.37  1.95 2.50       16.8       113    3.85       3.49
5    1   13.24  2.59 2.87       21.0       118    2.80       2.69
6    1   14.20  1.76 2.45       15.2       112    3.27       3.39
  Nonflavanoids Proanthocyanins Color  Hue Dilution Proline
1          0.28            2.29  5.64 1.04     3.92    1065
2          0.26            1.28  4.38 1.05     3.40    1050
3          0.30            2.81  5.68 1.03     3.17    1185
4          0.24            2.18  7.80 0.86     3.45    1480
5          0.39            1.82  4.32 1.04     2.93     735
6          0.34            1.97  6.75 1.05     2.85    1450
#str(wine)
#summary(wine)
df=scale(wine[-1])
require(NbClust)
set.seed(1234)
nc=NbClust(df,min.nc=2,max.nc=15,method="kmeans")

*** : 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:                                                
* 4 proposed 2 as the best number of clusters 
* 15 proposed 3 as the best number of clusters 
* 1 proposed 10 as the best number of clusters 
* 1 proposed 12 as the best number of clusters 
* 1 proposed 14 as the best number of clusters 
* 1 proposed 15 as the best number of clusters 

                   ***** Conclusion *****                            
 
* According to the majority rule, the best number of clusters is  3 
 
 
******************************************************************* 
table(nc$Best.n[1,])

 0  1  2  3 10 12 14 15 
 2  1  4 15  1  1  1  1 
par(mfrow=c(1,1))
barplot(table(nc$Best.n[1,]),xlab="Number of Clusters",ylab="Number of Criteria",
        main="Number of Clusters Chosen by 26 criteria")

또한 다음과 같은 코드로 군집의 갯수에 따른 total within-groups sums of squares를 계산해 plot을 그려보면 도움이 된다.

wssplot=function(data,nc=15,seed=1234,plot=TRUE){
        wss <- (nrow(data)-1)*sum(apply(data,2,var))
        
        for( i in 2:nc){
                set.seed(seed)
                wss[i]<-sum(kmeans(data,centers=i)$withinss)
        }
        if(plot) plot(1:nc,wss,type="b",xlab="Number pf Clusters",ylab="Within group sum of squares")
        wss
}

wssplot(df)

 [1] 2301.0000 1649.4400 1270.7491 1168.6143 1098.7390 1039.2957  977.5410
 [8]  952.5328  920.4558  883.7607  846.7963  806.6972  744.7018  729.1297
[15]  702.2454

wssplot에서 bend가 있는 곳이 적절한 군집의 갯수를 시사해준다. 적절한 군집의 갯수가 3개로 판단되므로 이를 이용해 k-means clustering을 시행한다.

fit.km=kmeans(df,3,nstart=25)
fit.km$cluster
  [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 3 3 3 3 3 3 3 3 3 2 2 1 2 2 2 2 2 2 2 2
 [71] 2 2 2 3 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 3 2 2 2 2 2 2 2 2 2
[106] 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 3 2 2 2 2 2 2 2 2 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
fit.km$size
[1] 51 65 62
fit.km$centers
     Alcohol      Malic        Ash Alcalinity   Magnesium     Phenols
1  0.1644436  0.8690954  0.1863726  0.5228924 -0.07526047 -0.97657548
2 -0.9234669 -0.3929331 -0.4931257  0.1701220 -0.49032869 -0.07576891
3  0.8328826 -0.3029551  0.3636801 -0.6084749  0.57596208  0.88274724
   Flavanoids Nonflavanoids Proanthocyanins      Color        Hue
1 -1.21182921    0.72402116     -0.77751312  0.9388902 -1.1615122
2  0.02075402   -0.03343924      0.05810161 -0.8993770  0.4605046
3  0.97506900   -0.56050853      0.57865427  0.1705823  0.4726504
    Dilution    Proline
1 -1.2887761 -0.4059428
2  0.2700025 -0.7517257
3  0.7770551  1.1220202
aggregate(wine[-1],by=list(clusters=fit.km$cluster),mean)
  clusters  Alcohol    Malic      Ash Alcalinity Magnesium  Phenols
1        1 13.13412 3.307255 2.417647   21.24118  98.66667 1.683922
2        2 12.25092 1.897385 2.231231   20.06308  92.73846 2.247692
3        3 13.67677 1.997903 2.466290   17.46290 107.96774 2.847581
  Flavanoids Nonflavanoids Proanthocyanins    Color       Hue Dilution
1  0.8188235     0.4519608        1.145882 7.234706 0.6919608 1.696667
2  2.0500000     0.3576923        1.624154 2.973077 1.0627077 2.803385
3  3.0032258     0.2920968        1.922097 5.453548 1.0654839 3.163387
    Proline
1  619.0588
2  510.1692
3 1100.2258

k-means clustering이 Type변수에 저장되어 있는 Type과 어는 정도 일치하는지 평가할 수 있다.

ct.km<-table(wine$Type,fit.km$cluster)
ct.km
   
     1  2  3
  1  0  0 59
  2  3 65  3
  3 48  0  0

Adjusted Rand Index는 Type과 cluster의 일치 정도를 정량화하는데 사용할 수 있다. -1은 전혀 일치하지 않는 것이고 1은 완벽하게 일치하는 것이다.

library(flexclust)
randIndex(ct.km)
     ARI 
0.897495 

Partitioning around medoids(PAM)

k-means clustering 은 평균을 이용하기 떄문에 이상치에 민감한 단점이 있다. 보다 강건한 방법은 partitioning around medoids(PAM) 방법이다. k-means clustering에서 각 군집을 centroid(변수들의 평균 벡터)로 나타내는 것과 달리 각 군집은 하나의 관찰치(medoid라고 부른다)로 대표된다. k-mean에서 유클리드 거리를 사용하는 것과 달리 PAM에서는 다른 거리 측정법도 사용할 수 있기 때문에 연속형 변수들 뿐만 아니라 mixed data type에도 적합시킬 수 있다.

PAM 알고리즘

  1. K개의 관찰치(medoid)를 무작위로 선택한다.
  2. 모든 관찰치에서 각medoid까지의 거리를 계산한다.
  3. 각 관찰치를 가장 가까운 medoid에 할당한다.
  4. 각 관찰치와 해당하는 medoid사이의 거리의 총합(총비용,total cost)을 계산한다.
  5. medoid가 아닌 점 하나를 선택하여 그 점에 할당된 medoid와 바꾼다.
  6. 모든 관찰치들을 가장 가까운 medoid에 할당한다.
  7. 총비용을 다시 계산한다.
  8. 다시계산한 총비용이 더 작다면 새 점들을 medoid로 유지한다.
  9. medoid가 바뀌지 않을 때까지 5-8단계를 반복한다.

PAM방법에서 사용하는 수학적인 방법에 대한 예는 https://en.wikipedia.org/wiki/K-medoids 를 참조한다.

R을 이용한 PAM

cluster패키지의 pam()함수를 이용하여 PAM 을 시행할 수 있다. 다음과 같은 형식으로 사용한다.

pam(x,k,metric="eucladean", stand=FALSE)

여기서 x는 데이터 행렬 또는 데이터프레임이고 k는 군집의 갯수, metric은 거리 측정 방법이고 stand는 거리를 측정하기 전에 변수들을 표준화할 것인지를 나타내는 논리값이다. wine 데이터에 PAM을 적용해보면 다음과 같다.

library(cluster)
set.seed(1234)
fit.pam <- pam(wine[-1],k=3,stand=TRUE)
fit.pam$medoids
     Alcohol Malic  Ash Alcalinity Magnesium Phenols Flavanoids
[1,]   13.48  1.81 2.41       20.5       100    2.70       2.98
[2,]   12.25  1.73 2.12       19.0        80    1.65       2.03
[3,]   13.40  3.91 2.48       23.0       102    1.80       0.75
     Nonflavanoids Proanthocyanins Color  Hue Dilution Proline
[1,]          0.26            1.86   5.1 1.04     3.47     920
[2,]          0.37            1.63   3.4 1.00     3.17     510
[3,]          0.43            1.41   7.3 0.70     1.56     750

PAM에서 사용되는 medoid는 wine 데이터에 포함되어 있는 실제 관측치인데 이 경우 36, 107, 175번쨰 관측치로 이 세 관측치가 세개의 군집을 대표한다.

clusplot(fit.pam, main="Bivariate Cluster Plot")

Bivariate plot은 각 관측치들을 두개의 주성분을 좌표로 하여 산점도로 나타낸 것이다. 각 군집은 각 군집에 속한 모든 관측치를 포함하는 가장 작은 타원으로 표시되어 있다.

이 경우 PAM의 수행능력은 k-means에 비해 떨어진다.

ct.pam=table(wine$Type, fit.pam$clustering)

randIndex(ct.pam)
      ARI 
0.6994957 

adjusted Rand index가 k-mean의 0.9에서 0.7로 감소되었다.

군집이 과연 존재하는가?

군집분석은 주어진 데이터 내에 존재하는 보다 가까운 subgoup을 발견해내는 방법이다. 이 방법은 매우 잘 동작한다. 사실 너무 잘 되기 때문에 존재하지도 않는 군집을 찾아내기도 한다. 다음의 코드를 보자.

library(fMultivar)
set.seed(1234)
df <- rnorm2d(1000,rho=.5)
df <- as.data.frame(df)

fMultivar 패키지의 rnorm2d 함수로 0.5의 상관을 갖는 이변량 정규분포를 갖는 1000개의 관측치를 무작위 추출하였다. 이를 산점도로 보면 데이터 내에 군집이 없다는 것을 알 수 있다.

plot(df, main="Bivariate normal Distribution with rho=0.5")

wssplot() 과 NbClust() 함수를 써서 군집의 갯수를 결정해본다.

wssplot(df)

 [1] 2010.8260 1018.6313  756.8010  601.0979  502.5761  425.6289  382.3161
 [8]  328.2312  287.1210  274.9018  252.4364  231.1382  211.3068  198.0347
[15]  183.4492
library(NbClust)
nc<-NbClust(df, min.nc=2, max.nc=15, method="kmeans")

*** : 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:                                                
* 8 proposed 2 as the best number of clusters 
* 5 proposed 3 as the best number of clusters 
* 1 proposed 4 as the best number of clusters 
* 1 proposed 5 as the best number of clusters 
* 1 proposed 8 as the best number of clusters 
* 4 proposed 10 as the best number of clusters 
* 1 proposed 12 as the best number of clusters 
* 2 proposed 13 as the best number of clusters 

                   ***** Conclusion *****                            
 
* According to the majority rule, the best number of clusters is  2 
 
 
******************************************************************* 
dev.new()
barplot(table(nc$Best.n[1,]),
        xlab="Number of Clusters",ylab="Number of Criteria",
        main="Numer of Clusters Chosen by 26 Criteria")

wssplot()은 세 개의 군집이 있다는 것을 시사해주고 NbClust()함수에 사용되는 많은 기준들은 두개 또는 세개의 군집이 있다는 것을 시사해준다. 이 결과에 근거하여 PAM을 이용해 두개의 군집으로 분석해보면 다음과 같다.

library(ggplot2)
library(cluster)
fit <-pam(df,k=2)
df$clustering <- factor(fit$clustering)
ggplot(data=df, aes(x=V1,y=V2,color=clustering, shape=clustering))+
        geom_point() + ggtitle("Clustering of Bivariate Normal Data")

이와 같은 군집은 분명히 인위적이다. 실제로는 군집이 존재하지 않는다. 이와 같은 실수를 피할 수 있는 방법이 있을까? 절대 안전한 것은 아니지만 NbClust패키지에서 제공하는 Cubic Cluster Criteria가 이러한 경우 도움이 될 수 있다.

plot(nc$All.index[,4], type="o", ylab="CCC", xlab="Number of clusters", col="blue")

CCC의 값이 모두 음수이고 두 개 이상의 군집에서 감소한다면 분포가 단봉형(unimodal)이라고 할 수 있다. 군집분석이 실제 존재하지 않는 군집을 만들어내는 능력이 있기 때문에 검증 단계가 필요하다. 군집분석에 의해 발견된 군집이 실제로 존재한다고 확인하기 위해서는 이러한 결과가 강건하고 반복적이어야 한다. 다른 군집분석 방법을 사용하고 다른 표본들에서 재현 가능한지 시도해보고 같은 결과가 나온다면 결과 대해 보다 확신을 가질 수 있을 것이다.

웹 R의 군집분석

http://r-ggplot2.com:3838/cluster