2 Cluster Analysis based on k-means

k-means의 장점

안정적인 결과를 찾을 때까지 각 반복마다 군집 간에 case 이동이 가능하다.
예측 변수가 많은 경우, 다른 군집분석 알고리듬에 비해 수행하는 시간이 적게 걸린다.
알고리듬 구현이 매우 간단하다.

k-means의 단점

범주형 예측 변수는 거리 계산을 할 수 없으므로 분석에 사용할 수 없다.
군집의 개수를 미리 결정해야 한다.
예측 변수들의 scale을 동일하게 설정해야 한다.
이상치에 민감하다.
초기 중심값 선택에 매우 민감하다.
비구형 군집에는 잘 작동하지 않는다.

실습 자료 : “USArrests” 데이터셋은 1973년 미국의 50개 주별로 살인(Murder), 폭행(Assault) 그리고 강간(Rape) 범죄의 10만명당 체포 건수와 도시 인구 비율을 포함하고 있다.

2.1 데이터 불러오기

pacman::p_load("data.table", 
               "tidyverse", 
               "dplyr",
               "caret",
               "GGally",                       # For ggpairs
               "mlr",
               "clue",
               "parallelMap",                  # For parallelStartSocket
               "parallel")                     # For detectCores

data("USArrests")                              # 데이터 불러오기

USArrests %>%
  as_tibble

# A tibble: 50 × 4
   Murder Assault UrbanPop  Rape
    <dbl>   <int>    <int> <dbl>
 1   13.2     236       58  21.2
 2   10       263       48  44.5
 3    8.1     294       80  31  
 4    8.8     190       50  19.5
 5    9       276       91  40.6
 6    7.9     204       78  38.7
 7    3.3     110       77  11.1
 8    5.9     238       72  15.8
 9   15.4     335       80  31.9
10   17.4     211       60  25.8
# ℹ 40 more rows

2.2 데이터 탐색

ggpairs(USArrests,
        upper = list(continuous = "density"),
        lower = list(continuous = wrap("points", size = 0.5)),
        diag = list(continuous = "densityDiag")) +
  theme_bw()

# 상관계수 그래프
ggcorr(USArrests,               # 데이터
       label = TRUE,            # 라벨 명시 여부
       label_round = 3,         # 상관계수 소숫점 이하 자릿수
       label_size = 3,          # 상관계수 글자 크기
       low = "steelblue",       # 상관계수가 음수일 때 색깔
       mid = "white",           # 상관계수가 0에 가까울 때 색깔
       high = "darkred")        # 상관계수가 양수일 때 색깔

2.3 데이터 분할

# Partition (Training Dataset : Test Dataset = 8:2)
set.seed(200)
ind <- sample(1:nrow(USArrests), 0.8*nrow(USArrests))       # Index를 이용하여 8:2로 분할

USArrests.trd <- USArrests[ind,]                            # Training Dataset
USArrests.ted <- USArrests[-ind,]                           # Test Dataset

2.4 데이터 전처리

# Standardization
preProcValues <- preProcess(USArrests.trd,
                            method = c("center", "scale"))  # Standardization 정의 -> Training Dataset에 대한 평균과 표준편차 계산 

USArrests.trd <- predict(preProcValues, USArrests.trd)      # Standardization for Training Dataset
USArrests.ted <- predict(preProcValues, USArrests.ted)      # Standardization for Test Dataset

glimpse(USArrests.trd)                                      # 데이터 구조 확인

Rows: 40
Columns: 4
$ Murder   <dbl> -0.28335915, -0.82987994, 1.64134450, -1.13878299, 1.11858549, -1.01997413, -0.35464447, 1.87896224, -0.75859462, 1.35620322, 2.04529465, 0.90472953, -0.87740349, -0.37840624, -0.21…
$ Assault  <dbl> -0.66262772, -0.20320584, 1.37532065, -1.06314937, 1.05725935, -0.49770705, -0.62728758, 1.02191920, -0.70974792, 0.30333625, 1.13971969, 1.62270167, -0.89822869, 0.89233867, -0.132…
$ UrbanPop <dbl> 0.56553882, 0.63579209, -1.12053964, 0.14401920, 1.19781824, 1.12756497, -0.76927330, 0.14401920, -0.13699387, -0.34775368, -1.40155272, 0.21427247, -1.33129945, 0.56553882, 0.28452…
$ Rape     <dbl> -0.62276898, 0.61022086, 0.20649852, -0.62276898, 2.77068093, 0.25014418, -0.45909777, 0.17376428, -0.44818635, 0.68660076, -0.38271786, 0.78480349, -0.85190869, -0.52456625, -0.066…

glimpse(USArrests.ted)                                      # 데이터 구조 확인

Rows: 10
Columns: 4
$ Murder   <dbl> 0.5958265, 1.8789622, 2.3541977, -1.1625448, 0.6908736, -0.3546445, -0.7348328, 0.3582087, 0.8572060, -0.9724506
$ Assault  <dbl> 1.1868399, 2.0350034, 0.5742774, -0.4977070, 1.0219192, -0.5566073, -0.1560856, 0.1855358, 1.0808194, 0.1384156
$ UrbanPop <dbl> -1.1205396, 1.1275650, -0.2775004, -0.6990200, 1.3383248, 0.1440192, 1.4788313, 0.4250323, 1.5490846, 1.6193379
$ Rape     <dbl> 2.6070097, 1.2321715, 0.5665752, -0.6991489, 0.3701697, -0.2845151, -0.4700092, 0.8284491, 0.5993094, -1.3429223

2.5 모형 훈련

2.5.1 Resampling을 이용한 최적의 초모수 조합 찾기

2.5.1.1 Define Task

문제 유형에 따라 Task를 정의하는 데 사용하는 함수는 다음과 같다.

문제 유형	함수
회귀 문제	`makeRegrTask()`
이진 또는 다중 클래스 분류 문제	`makeClassifTask()`
생존분석	`makeSurvTask()`
군집분석	`makeClusterTask()`
다중 라벨 분류 문제	`makeMultilabelTask()`
비용 민감 분류 문제	`makeCostSensTask()`

## k-means : 군집분석
USArrests.Task <- makeClusterTask(data = USArrests.trd)   # Training Dataset
USArrests.Task

Unsupervised task: USArrests.trd
Type: cluster
Observations: 40
Features:
   numerics     factors     ordered functionals 
          4           0           0           0 
Missings: FALSE
Has weights: FALSE
Has blocking: FALSE
Has coordinates: FALSE

2.5.1.2 Define Learner

학습자(Learner)는 함수 makeLearner()를 이용하여 정의한다.
함수 makeLearner()의 첫 번째 인자 cl에는 사용하고자 하는 머신러닝 알고리듬을 문제 유형.알고리듬과 관련된 R 함수 이름 형태로 입력한다.
- 예를 들어,
  - 회귀 문제 : "regr.알고리듬과 관련된 R 함수 이름"
  - 클래스 분류 문제 : "classif.알고리듬과 관련된 R 함수 이름"
  - 생존분석 : "surv.알고리듬과 관련된 R 함수 이름"
  - 군집분석 : "cluster.알고리듬과 관련된 R 함수 이름"
  - 다중 라벨 문제 : "multilabel.알고리듬과 관련된 R 함수 이름"
- Package "mlr"에서 사용할 수 있는 알고리듬은 여기를 통해서 확인할 수 있다.

# 초모수 집합
getParamSet("cluster.kmeans")  # cluster.kmeans : Clustering based on "kmeans" function

              Type len           Def                             Constr Req Tunable Trafo
centers    untyped   -             -                                  -   -    TRUE     -
iter.max   integer   -            10                           1 to Inf   -    TRUE     -
nstart     integer   -             1                           1 to Inf   -    TRUE     -
algorithm discrete   - Hartigan-Wong Hartigan-Wong,Lloyd,Forgy,MacQueen   -    TRUE     -
trace      logical   -             -                                  -   -   FALSE     -

Caution! 특정 머신러닝 알고리듬이 가지고 있는 초모수는 함수 getParamSet()를 이용하여 확인할 수 있다.

USArrests.Learner <- makeLearner(cl = "cluster.kmeans",                # 함수 kmeans를 이용하여 군집분석 수행
                                 par.vals = list(iter.max = 100,       # 최대 반복 수
                                                 nstart = 25))         # 수행 횟수
USArrests.Learner

Learner cluster.kmeans from package stats,clue
Type: cluster
Name: K-Means; Short name: kmeans
Class: cluster.kmeans
Properties: numerics,prob
Predict-Type: response
Hyperparameters: centers=2,iter.max=100,nstart=25

Caution! 함수 makeLearner()의 인자 par.vals에 함수 kmeans()의 옵션을 입력할 수 있다. iter.max에는 최대 반복 수를, nstart에는 k-means 수행 횟수를 지정할 수 있다. k-means는 초기 중심값을 랜덤하게 선택하기 때문에 이 과정에서 다양한 결과가 나타날 수 있다. 그래서 옵션 nstart를 이용하여 수행 횟수를 늘려 최대한 다양한 초기 중심값에 대해 k-means를 수행하고 최적의 결과를 찾을 수 있다.

2.5.1.3 Define Search Space

함수 makeDiscreteParam() 또는 makeNumericParam()을 이용하여 초모수의 검색 범위를 정의한다.
- 함수 makeDiscreteParam() : 검색 범위를 특정값으로 정의하는 경우 사용
- 함수 makeNumericParam() : 검색 범위를 구간으로 정의하는 경우 사용
그러고나서, 함수 makeParamSet()를 이용하여 정의한 검색 범위을 ParamSet 객체로 만든다.

# 초모수 "centers" (군집 수)와 알고리듬의 검색 범위 정의 
tune.hyper <- makeParamSet( 
  makeDiscreteParam("centers",                                         # 군집 수 
                    values = 3:7),                                     # 군집 수에 대한 검색 범위
  makeDiscreteParam("algorithm",                                       # k-means 알고리듬
                    values = c("Lloyd", "MacQueen", "Hartigan-Wong"))) # 알고리듬에 대한 검색 범위
tune.hyper

              Type len Def                       Constr Req Tunable Trafo
centers   discrete   -   -                    3,4,5,6,7   -    TRUE     -
algorithm discrete   -   - Lloyd,MacQueen,Hartigan-Wong   -    TRUE     -

Caution! k-means를 수행하기 위해 “Lloyd”, “MacQueen” 그리고 “Hartigan-Wong” 알고리듬을 검색 범위로 정의하였다. 각 알고리듬의 수행 절차는 다음과 같다.

2.5.1.4 Define Tuning Method

정의한 검색 범위 내에서 후보 초모수 조합을 찾는다.
가장 많이 사용하는 방법은 그리드 검색(Grid Search)과 랜덤 검색(Random Search)이며, 각각 함수 makeTuneControlGrid()와 makeTuneControlRandom()을 이용하여 정의할 수 있다.

gridSearch <- makeTuneControlGrid()                                    # 그리드 검색               
gridSearch

Tune control: TuneControlGrid
Same resampling instance: TRUE
Imputation value: <worst>
Start: <NULL>

Tune threshold: FALSE
Further arguments: resolution=10

Result! 함수 makeTuneControlGrid()를 이용하여 위에서 정의한 초모수 “centers”와 알고리듬의 검색 범위에 해당하는 모든 조합을 후보 초모수 조합으로 설정한다.

2.5.1.5 Define Resampling Strategy

kFold <- makeResampleDesc("CV", iters = 5)                             # 5-Fold Cross Validation
kFold

Resample description: cross-validation with 5 iterations.
Predict: test
Stratification: FALSE

2.5.1.6 Perform Tuning

최적의 초모수 조합을 찾기 위해 함수 tuneParams()를 이용한다.

parallelStartSocket(cpus = detectCores())                              # 병렬 처리

set.seed(100)
tunedK <- tuneParams(task = USArrests.Task,                            # Defined Task in 5-1-1
                     learner = USArrests.Learner,                      # Defined Learner in 5-1-2
                     par.set = tune.hyper,                             # Defined Search Space in 5-1-3
                     control = gridSearch,                             # Defined Tuning Method in 5-1-4
                     resampling = kFold,                               # Defined Resampling Strategy in 5-1-5
                     measures = list(db))                              # Davies-Bouldin Index (군집내 유사성이 높고 군집간 유사성이 낮을수록 값이 낮음)

tunedK

Tune result:
Op. pars: centers=7; algorithm=Hartigan-Wong
db.test.mean=0.5285825

Result! 최적의 초모수 조합을 찾기 위해 “Silhouette Index”를 사용하였다. “Silhouette Index”는 다른 군집과 비교하여 case가 현재 속한 군집과 얼마나 유사한지를 나타내는 척도로 값이 높을수록 군집화가 잘 되었다는 것을 의미한다.

tunedK$x$centers                                                       # 최적의 군집 수

[1] 7

tunedK$x$algorithm                                                     # 최적의 알고리듬

[1] "Hartigan-Wong"

# 튜닝 과정 시각화
kMeansTuningData <- generateHyperParsEffectData(tunedK)                # Extract Hyperparameter Effect from Tuning Result
kMeansTuningData$data

   centers     algorithm db.test.mean iteration exec.time
1        3         Lloyd    0.7914754         1      0.58
2        4         Lloyd    0.7204947         2      0.72
3        5         Lloyd    0.7140210         3      0.70
4        6         Lloyd    0.6604494         4      0.61
5        7         Lloyd    0.5892118         5      0.62
6        3      MacQueen    0.7662724         6      0.66
7        4      MacQueen    0.7349811         7      0.06
8        5      MacQueen    0.7286601         8      0.08
9        6      MacQueen    0.7524750         9      0.09
10       7      MacQueen    0.5848960        10      0.08
11       3 Hartigan-Wong    0.7662724        11      0.09
12       4 Hartigan-Wong    0.7349811        12      0.08
13       5 Hartigan-Wong    0.7267509        13      0.11
14       6 Hartigan-Wong    0.7990965        14      0.12
15       7 Hartigan-Wong    0.5285825        15      0.12

# 데이터 구조 변환
TuningData <- pivot_longer(kMeansTuningData$data,
                           cols = -c(centers, iteration, algorithm, exec.time),
                           names_to = "Metric",
                           values_to = "Value")

TuningData %>%
  as_tibble

# A tibble: 15 × 6
   centers algorithm     iteration exec.time Metric       Value
     <int> <chr>             <int>     <dbl> <chr>        <dbl>
 1       3 Lloyd                 1    0.580  db.test.mean 0.791
 2       4 Lloyd                 2    0.72   db.test.mean 0.720
 3       5 Lloyd                 3    0.7    db.test.mean 0.714
 4       6 Lloyd                 4    0.61   db.test.mean 0.660
 5       7 Lloyd                 5    0.620  db.test.mean 0.589
 6       3 MacQueen              6    0.66   db.test.mean 0.766
 7       4 MacQueen              7    0.0600 db.test.mean 0.735
 8       5 MacQueen              8    0.0800 db.test.mean 0.729
 9       6 MacQueen              9    0.0900 db.test.mean 0.752
10       7 MacQueen             10    0.0800 db.test.mean 0.585
11       3 Hartigan-Wong        11    0.0900 db.test.mean 0.766
12       4 Hartigan-Wong        12    0.0800 db.test.mean 0.735
13       5 Hartigan-Wong        13    0.110  db.test.mean 0.727
14       6 Hartigan-Wong        14    0.120  db.test.mean 0.799
15       7 Hartigan-Wong        15    0.120  db.test.mean 0.529

ggplot(TuningData, aes(centers, Value, col = algorithm)) +
  facet_wrap(~ Metric, scales = "free_y") +
  geom_line() +
  geom_point() +
  theme_bw()

Result! 초모수 “centers = 4”이고 알고리듬이 “Lloyd”일 때 “Silhouette Index”값이 가장 높다는 것을 알 수 있다.

2.5.2 최적의 초모수 조합과 함께 모형 훈련

2.5.2.1 Redefine Learner

# Redefine Learner with 최적의 초모수 조합
tunedKMeans <- setHyperPars(USArrests.Learner,                               # Defined Learner in 5-1-2
                            par.vals = list(centers = tunedK$x$centers,      # 최적의 군집 수
                                            algorithm = tunedK$x$algorithm)) # 최적의 알고리듬
tunedKMeans

Learner cluster.kmeans from package stats,clue
Type: cluster
Name: K-Means; Short name: kmeans
Class: cluster.kmeans
Properties: numerics,prob
Predict-Type: response
Hyperparameters: centers=7,iter.max=100,nstart=25,algorithm=Hartigan-Wong

2.5.2.2 Train Model

tunedKMeansModel <- train(tunedKMeans,                                       # Defined Learner in 5-2-1
                          USArrests.Task)                                    # Defined Task in 5-1-1
tunedKMeansModel

Model for learner.id=cluster.kmeans; learner.class=cluster.kmeans
Trained on: task.id = USArrests.trd; obs = 40; features = 4
Hyperparameters: centers=7,iter.max=100,nstart=25,algorithm=Hartigan-Wong

kMeanModel <- getLearnerModel(tunedKMeansModel)                              # 예측 모형 추출
kMeanModel

K-means clustering with 7 clusters of sizes 5, 6, 6, 7, 7, 6, 3

Cluster means:
       Murder    Assault   UrbanPop       Rape
1  0.56255999  1.1067356  1.1837676  1.9283197
2  1.27699731  0.9532023  0.1908547  0.5829423
3  0.08098804 -0.2660328 -0.4765513 -0.2572366
4 -1.13878299 -1.1439268 -1.2409738 -1.1698985
5 -0.44629702 -0.0736253  0.7662624  0.1581765
6 -0.81007846 -0.9512389  0.4952855 -0.6409547
7  1.66510628  1.5245346 -1.2844639 -0.2226838

Clustering vector:
  Pennsylvania     Washington South Carolina      Minnesota         Nevada           Utah        Montana      Louisiana       Nebraska      Tennessee    Mississippi       Maryland   South Dakota 
             6              5              7              6              1              5              3              2              6              2              7              2              4 
      Delaware       Oklahoma     New Mexico       Colorado     New Jersey      Wisconsin       Virginia   North Dakota        Arizona         Hawaii  New Hampshire     California        Alabama 
             5              5              2              1              5              6              3              4              1              6              4              1              2 
      Michigan           Ohio         Oregon    Connecticut          Texas       Arkansas North Carolina          Maine        Indiana        Wyoming       Kentucky           Iowa        Vermont 
             1              5              5              6              2              3              7              4              3              3              3              4              4 
 West Virginia 
             4 

Within cluster sum of squares by cluster:
[1] 3.970393 4.544595 2.689821 3.834536 5.417810 3.459916 1.052706
 (between_SS / total_SS =  84.0 %)

Available components:

[1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss" "betweenss"    "size"         "iter"         "ifault"

# Training Dataset의 각 군집별 특징 시각화
USArrests.trd.clus <- mutate(USArrests.trd,
                             kMeansCluster = as.factor(kMeanModel$cluster))  # 예측 모형에 의한 군집 결과
USArrests.trd.clus %>%
  as_tibble

# A tibble: 40 × 5
   Murder Assault UrbanPop   Rape kMeansCluster
    <dbl>   <dbl>    <dbl>  <dbl> <fct>        
 1 -0.283  -0.663    0.566 -0.623 6            
 2 -0.830  -0.203    0.636  0.610 5            
 3  1.64    1.38    -1.12   0.206 7            
 4 -1.14   -1.06     0.144 -0.623 6            
 5  1.12    1.06     1.20   2.77  1            
 6 -1.02   -0.498    1.13   0.250 5            
 7 -0.355  -0.627   -0.769 -0.459 3            
 8  1.88    1.02     0.144  0.174 2            
 9 -0.759  -0.710   -0.137 -0.448 6            
10  1.36    0.303   -0.348  0.687 2            
# ℹ 30 more rows

ggpairs(USArrests.trd.clus, aes(col = kMeansCluster),
        upper = list(continuous = "density")) +
  theme_bw()

2.6 예측

# 예측 군집 생성
Pred <- predict(tunedKMeansModel, 
                newdata = USArrests.ted)         # predict(Trained Model, Test Dataset)
Pred %>%
  as_tibble

# A tibble: 10 × 1
   response
      <int>
 1        1
 2        2
 3        2
 4        4
 5        2
 6        6
 7        5
 8        5
 9        1
10        6