クラスター分析(基礎編)
基礎編
使用パッケージ
＜ポイント＞
・{cluster} ：クラスタリングアルゴリズムを計算
・{factoextra} ：クラスタリング結果のggplot2ベースのエレガントな可視化データ基準化
概要
＜ポイント＞
・データセットの行は観測値(個体)であり、列は変数であること
・データに欠損値があれば除去または推定する必要がある
・変数を比較可能にするためには、標準化する必要がある
--- 平均値がゼロで標準偏差が1になるように変数を変換データロード
# データロード
data("USArrests")
# データ格納
df <- USArrests %>% as_tibble()
df## # A tibble: 50 x 4
## Murder Assault UrbanPop Rape
## * <dbl> <int> <int> <dbl>
## 1 13.2 236 58 21.2
## 2 10.0 263 48 44.5
## 3 8.10 294 80 31.0
## 4 8.80 190 50 19.5
## 5 9.00 276 91 40.6
## 6 7.90 204 78 38.7
## 7 3.30 110 77 11.1
## 8 5.90 238 72 15.8
## 9 15.4 335 80 31.9
## 10 17.4 211 60 25.8
## # ... with 40 more rows欠損値除去
df2 <- df %>% na.omit()
df2## # A tibble: 50 x 4
## Murder Assault UrbanPop Rape
## <dbl> <int> <int> <dbl>
## 1 13.2 236 58 21.2
## 2 10.0 263 48 44.5
## 3 8.10 294 80 31.0
## 4 8.80 190 50 19.5
## 5 9.00 276 91 40.6
## 6 7.90 204 78 38.7
## 7 3.30 110 77 11.1
## 8 5.90 238 72 15.8
## 9 15.4 335 80 31.9
## 10 17.4 211 60 25.8
## # ... with 40 more rows基準化
# 基準化
X <- df2 %>% scale()
X## Murder Assault UrbanPop Rape
## [1,] 1.24256408 0.78283935 -0.52090661 -0.003416473
## [2,] 0.50786248 1.10682252 -1.21176419 2.484202941
## [3,] 0.07163341 1.47880321 0.99898006 1.042878388
## [4,] 0.23234938 0.23086801 -1.07359268 -0.184916602
## [5,] 0.27826823 1.26281442 1.75892340 2.067820292
## [6,] 0.02571456 0.39885929 0.86080854 1.864967207
## [7,] -1.03041900 -0.72908214 0.79172279 -1.081740768
## [8,] -0.43347395 0.80683810 0.44629400 -0.579946294
## [9,] 1.74767144 1.97077766 0.99898006 1.138966691
## [10,] 2.20685994 0.48285493 -0.38273510 0.487701523
## [11,] -0.57123050 -1.49704226 1.20623733 -0.110181255
## [12,] -1.19113497 -0.60908837 -0.79724965 -0.750769945
## [13,] 0.59970018 0.93883125 1.20623733 0.295524916
## [14,] -0.13500142 -0.69308401 -0.03730631 -0.024769429
## [15,] -1.28297267 -1.37704849 -0.58999237 -1.060387812
## [16,] -0.41051452 -0.66908525 0.03177945 -0.345063775
## [17,] 0.43898421 -0.74108152 -0.93542116 -0.526563903
## [18,] 1.74767144 0.93883125 0.03177945 0.103348309
## [19,] -1.30593210 -1.05306531 -1.00450692 -1.434064548
## [20,] 0.80633501 1.55079947 0.10086521 0.701231086
## [21,] -0.77786532 -0.26110644 1.34440885 -0.526563903
## [22,] 0.99001041 1.01082751 0.58446551 1.480613993
## [23,] -1.16817555 -1.18505846 0.03177945 -0.676034598
## [24,] 1.90838741 1.05882502 -1.48810723 -0.441152078
## [25,] 0.27826823 0.08687549 0.30812248 0.743936999
## [26,] -0.41051452 -0.74108152 -0.86633540 -0.515887425
## [27,] -0.80082475 -0.82507715 -0.24456358 -0.505210947
## [28,] 1.01296983 0.97482938 1.06806582 2.644350114
## [29,] -1.30593210 -1.36504911 -0.65907813 -1.252564419
## [30,] -0.08908257 -0.14111267 1.62075188 -0.259651949
## [31,] 0.82929443 1.37080881 0.30812248 1.160319648
## [32,] 0.76041616 0.99882813 1.41349461 0.519730957
## [33,] 1.19664523 1.99477641 -1.41902147 -0.547916860
## [34,] -1.60440462 -1.50904164 -1.48810723 -1.487446939
## [35,] -0.11204199 -0.60908837 0.65355127 0.017936483
## [36,] -0.27275797 -0.23710769 0.16995096 -0.131534211
## [37,] -0.66306820 -0.14111267 0.10086521 0.861378259
## [38,] -0.34163624 -0.77707965 0.44629400 -0.676034598
## [39,] -1.00745957 0.03887798 1.48258036 -1.380682157
## [40,] 1.51807718 1.29881255 -1.21176419 0.135377743
## [41,] -0.91562187 -1.01706718 -1.41902147 -0.900240639
## [42,] 1.24256408 0.20686926 -0.45182086 0.605142783
## [43,] 1.12776696 0.36286116 0.99898006 0.455672088
## [44,] -1.05337842 -0.60908837 0.99898006 0.178083656
## [45,] -1.28297267 -1.47304350 -2.31713632 -1.071064290
## [46,] 0.16347111 -0.17711080 -0.17547783 -0.056798864
## [47,] -0.86970302 -0.30910395 0.51537975 0.530407436
## [48,] -0.47939280 -1.07706407 -1.83353601 -1.273917376
## [49,] -1.19113497 -1.41304662 0.03177945 -1.113770203
## [50,] -0.22683912 -0.11711392 -0.38273510 -0.601299251
## attr(,"scaled:center")
## Murder Assault UrbanPop Rape
## 7.788 170.760 65.540 21.232
## attr(,"scaled:scale")
## Murder Assault UrbanPop Rape
## 4.355510 83.337661 14.474763 9.366385# 確認：平均
X %>% apply(2, mean) %>% round(4)## Murder Assault UrbanPop Rape
## 0 0 0 0# 確認：標準偏差
X %>% apply(2, sd) %>% round(4)## Murder Assault UrbanPop Rape
## 1 1 1 1一連の処理
USArrests %>%
drop_na() %>%
scale()## Murder Assault UrbanPop Rape
## Alabama 1.24256408 0.78283935 -0.52090661 -0.003416473
## Alaska 0.50786248 1.10682252 -1.21176419 2.484202941
## Arizona 0.07163341 1.47880321 0.99898006 1.042878388
## Arkansas 0.23234938 0.23086801 -1.07359268 -0.184916602
## California 0.27826823 1.26281442 1.75892340 2.067820292
## Colorado 0.02571456 0.39885929 0.86080854 1.864967207
## Connecticut -1.03041900 -0.72908214 0.79172279 -1.081740768
## Delaware -0.43347395 0.80683810 0.44629400 -0.579946294
## Florida 1.74767144 1.97077766 0.99898006 1.138966691
## Georgia 2.20685994 0.48285493 -0.38273510 0.487701523
## Hawaii -0.57123050 -1.49704226 1.20623733 -0.110181255
## Idaho -1.19113497 -0.60908837 -0.79724965 -0.750769945
## Illinois 0.59970018 0.93883125 1.20623733 0.295524916
## Indiana -0.13500142 -0.69308401 -0.03730631 -0.024769429
## Iowa -1.28297267 -1.37704849 -0.58999237 -1.060387812
## Kansas -0.41051452 -0.66908525 0.03177945 -0.345063775
## Kentucky 0.43898421 -0.74108152 -0.93542116 -0.526563903
## Louisiana 1.74767144 0.93883125 0.03177945 0.103348309
## Maine -1.30593210 -1.05306531 -1.00450692 -1.434064548
## Maryland 0.80633501 1.55079947 0.10086521 0.701231086
## Massachusetts -0.77786532 -0.26110644 1.34440885 -0.526563903
## Michigan 0.99001041 1.01082751 0.58446551 1.480613993
## Minnesota -1.16817555 -1.18505846 0.03177945 -0.676034598
## Mississippi 1.90838741 1.05882502 -1.48810723 -0.441152078
## Missouri 0.27826823 0.08687549 0.30812248 0.743936999
## Montana -0.41051452 -0.74108152 -0.86633540 -0.515887425
## Nebraska -0.80082475 -0.82507715 -0.24456358 -0.505210947
## Nevada 1.01296983 0.97482938 1.06806582 2.644350114
## New Hampshire -1.30593210 -1.36504911 -0.65907813 -1.252564419
## New Jersey -0.08908257 -0.14111267 1.62075188 -0.259651949
## New Mexico 0.82929443 1.37080881 0.30812248 1.160319648
## New York 0.76041616 0.99882813 1.41349461 0.519730957
## North Carolina 1.19664523 1.99477641 -1.41902147 -0.547916860
## North Dakota -1.60440462 -1.50904164 -1.48810723 -1.487446939
## Ohio -0.11204199 -0.60908837 0.65355127 0.017936483
## Oklahoma -0.27275797 -0.23710769 0.16995096 -0.131534211
## Oregon -0.66306820 -0.14111267 0.10086521 0.861378259
## Pennsylvania -0.34163624 -0.77707965 0.44629400 -0.676034598
## Rhode Island -1.00745957 0.03887798 1.48258036 -1.380682157
## South Carolina 1.51807718 1.29881255 -1.21176419 0.135377743
## South Dakota -0.91562187 -1.01706718 -1.41902147 -0.900240639
## Tennessee 1.24256408 0.20686926 -0.45182086 0.605142783
## Texas 1.12776696 0.36286116 0.99898006 0.455672088
## Utah -1.05337842 -0.60908837 0.99898006 0.178083656
## Vermont -1.28297267 -1.47304350 -2.31713632 -1.071064290
## Virginia 0.16347111 -0.17711080 -0.17547783 -0.056798864
## Washington -0.86970302 -0.30910395 0.51537975 0.530407436
## West Virginia -0.47939280 -1.07706407 -1.83353601 -1.273917376
## Wisconsin -1.19113497 -1.41304662 0.03177945 -1.113770203
## Wyoming -0.22683912 -0.11711392 -0.38273510 -0.601299251
## attr(,"scaled:center")
## Murder Assault UrbanPop Rape
## 7.788 170.760 65.540 21.232
## attr(,"scaled:scale")
## Murder Assault UrbanPop Rape
## 4.355510 83.337661 14.474763 9.366385距離測定
距離の測定方法
＜ポイント＞
・距離尺度の選択はクラスタリングの重要なステップ
--- 2つの要素(x, y)の類似度がどのように計算されるかを定義
・距離測定の前にはデータ標準化が必須
＜古典的な距離測定の方法＞
・一般的なクラスタリング関数では、ユークリッド距離がデフォルト設定が多い
--- ユークリッド距離 : 直線距離 (ピタゴラス距離)
--- マンハッタン距離 : 市街地距離
＜相関に基づく距離＞
・ユークリッド距離で離れていても、相関距離は相関があれば類似しているとみなす
--- ピアソン相関距離 : 2つのプロファイル間の線形関係の度合いを測定
--- コサイン距離 : ピアソン相関距離の変形
--- スピアマン相関距離 : 順位情報に基づく相関
--- ケンドール相関距離 : 順位の対応関係に基づく相関距離測定の関数
＜ポイント＞
1. stats::dist()
・数値データのみを入力として受け入れて距離行列を作成
・基本的な距離尺度をカバー
--- ユークリッド/最大化/マンハッタン/キャンベラ/バイナリ/ミンコフスキー
2. factoextra::get_dist()
・数値データのみを入力として受け入れて距離行列を作成
・相関ベースの距離測定をサポート
3. cluster::daisy()
・数値以外の変数型を扱うことができるデータ準備
＜ポイント＞
・USArrestsデータセットを使用する
・50行から15行をランダムに選択してデータセットを縮小# 乱数設定
set.seed(123)
# 行をランダムに選択
ss <- 1:50 %>% sample(15)
df <- USArrests[ss, ]
df %>% head()## Murder Assault UrbanPop Rape
## Iowa 2.2 56 57 11.3
## Rhode Island 3.4 174 87 8.3
## Maryland 11.3 300 67 27.8
## Tennessee 13.2 188 59 26.9
## Utah 3.2 120 80 22.9
## Arizona 8.1 294 80 31.0# データ標準化
df.scaled <- df %>% scale()
df.scaled %>% head()## Murder Assault UrbanPop Rape
## Iowa -1.2508408 -1.49648631 -0.6270029 -1.0725081
## Rhode Island -1.0079591 -0.06238736 1.6857126 -1.4332171
## Maryland 0.5910122 1.46893865 0.1439023 0.9113914
## Tennessee 0.9755749 0.10775998 -0.4728218 0.8031787
## Utah -1.0484394 -0.71866993 1.1460790 0.3222334
## Arizona -0.0566724 1.39601836 1.1460790 1.2961477ユークリッド距離
＜ポイント＞
・ユークリッド距離などの基本的な距離尺度はdist()で算出できる
・dist()関数で生成された距離情報はあまり見やすくない
--- as.matrix()関数を使用して距離ベクトルを行列に再フォーマットすることができる
・個別要素ごとの距離を測定している
--- 変数間のペアワイズ距離を計算したい場合はt()で転置して計算# 元データのサイズ
df.scaled %>% dim()## [1] 15 4# ユークリッド距離に基づく距離行列
# --- 下三角行列で示される
# --- print()しても見やすくない
dist.eucl <- df.scaled %>% dist(method = "euclidean")
dist.eucl %>% class()## [1] "dist"# matrixクラスに変換
X <- dist.eucl %>% as.matrix()
X %>% class()## [1] "matrix"# 表示
X[1:5, 1:5] %>% round(1)## Iowa Rhode Island Maryland Tennessee Utah
## Iowa 0.0 2.8 4.1 3.3 2.4
## Rhode Island 2.8 0.0 3.6 3.7 2.0
## Maryland 4.1 3.6 0.0 1.5 3.0
## Tennessee 3.3 3.7 1.5 0.0 2.8
## Utah 2.4 2.0 3.0 2.8 0.0# 変数間のペアワイズ距離
df.scaled %>% t() %>% dist(method = "euclidean") %>% as.matrix()## Murder Assault UrbanPop Rape
## Murder 0.000000 2.479688 5.628043 3.228057
## Assault 2.479688 0.000000 4.681875 2.757250
## UrbanPop 5.628043 4.681875 0.000000 4.355974
## Rape 3.228057 2.757250 4.355974 0.000000相関に基づく距離
＜ポイント＞
・get_dist()を使用するとdist()の距離尺度以外に、相関に基づく距離が計算できる
--- "pearson", "spearman" or "kendall"# 距離行列の作成
dist.cor <- df.scaled %>% get_dist(method = "pearson")
# Display a subset
dist.cor %>% as.matrix() %>% .[1:5, 1:5] %>% round(1)## Iowa Rhode Island Maryland Tennessee Utah
## Iowa 0.0 0.4 1.9 1.5 0.1
## Rhode Island 0.4 0.0 1.5 2.0 0.4
## Maryland 1.9 1.5 0.0 0.7 1.6
## Tennessee 1.5 2.0 0.7 0.0 1.6
## Utah 0.1 0.4 1.6 1.6 0.0混合データの距離
＜ポイント＞
・cluster::daisy()を使用すると数値データ以外についても距離行列を作成できる# データロード
data(flower)
# データ確認
flower %>% head(3)## V1 V2 V3 V4 V5 V6 V7 V8
## 1 0 1 1 4 3 15 25 15
## 2 1 0 0 2 1 3 150 50
## 3 0 1 0 3 3 1 150 50# データ構造
flower %>% glimpse()## Observations: 18
## Variables: 8
## $ V1 <fct> 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0
## $ V2 <fct> 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0
## $ V3 <fct> 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1
## $ V4 <fct> 4, 2, 3, 4, 5, 4, 4, 2, 3, 5, 5, 1, 1, 4, 3, 4, 2, 2
## $ V5 <ord> 3, 1, 3, 2, 2, 3, 3, 2, 1, 2, 3, 2, 2, 2, 2, 2, 2, 1
## $ V6 <ord> 15, 3, 1, 16, 2, 12, 13, 7, 4, 14, 8, 9, 6, 11, 10, 18, 17, 5
## $ V7 <dbl> 25, 150, 150, 125, 20, 50, 40, 100, 25, 100, 45, 90, 20, 80...
## $ V8 <dbl> 15, 50, 50, 50, 15, 40, 20, 15, 15, 60, 10, 25, 10, 30, 20,...# 距離行列の作成
dd <- flower %>% daisy()距離行列の可視化
＜ポイント＞
・fviz_dist()を使うと距離行列をヒートマップで可視化することができる
--- データセットはdistクラス# 距離行列をヒートマップ化
df %>%
scale() %>%
dist(method = "euclidean") %>%
fviz_dist()
分析４ステップ
はじめに
＜ポイント＞
・stats::k-means()を使用してパーティション化クラスタリングを実装する
・k-means()の適用前後で行うべき準備や検証を確認する
・factoextra::eclust()を使用すると以下のワークフローの簡素化が可能
＜分析ステップ＞
・1：クラスタ性の評価
・2：クラスター数を見積り
・3：meansの計算
・4：クラスタ検証統計データ準備
＜ポイント＞
・USArrestsデータセットを使用# データロード
data(USArrests)
# データサイズ
USArrests %>% dim()## [1] 50 4# 基準化
df <- USArrests %>% scale()
df %>% head()## Murder Assault UrbanPop Rape
## Alabama 1.24256408 0.7828393 -0.5209066 -0.003416473
## Alaska 0.50786248 1.1068225 -1.2117642 2.484202941
## Arizona 0.07163341 1.4788032 0.9989801 1.042878388
## Arkansas 0.23234938 0.2308680 -1.0735927 -0.184916602
## California 0.27826823 1.2628144 1.7589234 2.067820292
## Colorado 0.02571456 0.3988593 0.8608085 1.8649672071：クラスタ性の評価
＜ポイント＞
・ホプキンス統計量はデータセットにクラスタ性が存在するかを示す
--- 有意水準の0.5より小さいと、クラスタが存在すると判断
・順序付き類似度画像(DOI)も併せて表示される# クラスタ性の評価
res <- df %>% get_clust_tendency(40, graph = TRUE)
# ホプキンス統計量
# --- 有意水準は0.5
res$hopkins_stat## [1] 0.3440875# プロット
res$plot
2：クラスター数を見積り
＜ポイント＞
・k-meansクラスタリングなどでは、生成クラスタ数を指定するためにが必要がある
・ギャップ統計量でクラスタの最適数を見積もるための計算をする
--- ブートストラップ法で計算
・今回はクラスタ数は４が適切であることを示唆# 乱数設定
set.seed(123)
# Compute the gap statistic
gap_stat <-
df %>%
clusGap(FUN = kmeans, nstart = 25, K.max = 10, B = 500) %T>%
print()## Clustering Gap statistic ["clusGap"] from call:
## clusGap(x = ., FUNcluster = kmeans, K.max = 10, B = 500, nstart = 25)
## B=500 simulated reference sets, k = 1..10; spaceH0="scaledPCA"
## --> Number of clusters (method 'firstSEmax', SE.factor=1): 4
## logW E.logW gap SE.sim
## [1,] 3.458369 3.639717 0.1813481 0.03833784
## [2,] 3.135112 3.369717 0.2346052 0.03511187
## [3,] 2.977727 3.231689 0.2539628 0.03559512
## [4,] 2.826221 3.116905 0.2906841 0.03618572
## [5,] 2.738868 3.019569 0.2807013 0.03781053
## [6,] 2.669860 2.933315 0.2634555 0.03776678
## [7,] 2.598748 2.854657 0.2559083 0.03844657
## [8,] 2.531626 2.782849 0.2512235 0.03936957
## [9,] 2.468162 2.715364 0.2472015 0.04045979
## [10,] 2.394884 2.651371 0.2564875 0.04184419# ギャップ統計量の可視化
gap_stat %>% fviz_gap_stat()
3：k-meansの計算
＜ポイント＞
・Gap統計量に基づいてk-meansを実行
--- 乱数シードを検証時と一致させている点に注意# 乱数設定
set.seed(123)
# k平均法による句ｒスタリング
km.res <- df %>% kmeans(4, nstart = 25)
km.res %>% print()## K-means clustering with 4 clusters of sizes 13, 16, 13, 8
##
## Cluster means:
## Murder Assault UrbanPop Rape
## 1 -0.9615407 -1.1066010 -0.9301069 -0.96676331
## 2 -0.4894375 -0.3826001 0.5758298 -0.26165379
## 3 0.6950701 1.0394414 0.7226370 1.27693964
## 4 1.4118898 0.8743346 -0.8145211 0.01927104
##
## Clustering vector:
## Alabama Alaska Arizona Arkansas California
## 4 3 3 4 3
## Colorado Connecticut Delaware Florida Georgia
## 3 2 2 3 4
## Hawaii Idaho Illinois Indiana Iowa
## 2 1 3 2 1
## Kansas Kentucky Louisiana Maine Maryland
## 2 1 4 1 3
## Massachusetts Michigan Minnesota Mississippi Missouri
## 2 3 1 4 3
## Montana Nebraska Nevada New Hampshire New Jersey
## 1 1 3 1 2
## New Mexico New York North Carolina North Dakota Ohio
## 3 3 4 1 2
## Oklahoma Oregon Pennsylvania Rhode Island South Carolina
## 2 2 2 2 4
## South Dakota Tennessee Texas Utah Vermont
## 1 4 3 2 1
## Virginia Washington West Virginia Wisconsin Wyoming
## 2 2 1 1 2
##
## Within cluster sum of squares by cluster:
## [1] 11.952463 16.212213 19.922437 8.316061
## (between_SS / total_SS = 71.2 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss"
## [5] "tot.withinss" "betweenss" "size" "iter"
## [9] "ifault"# 所属クラスタの確認
km.res$cluster %>% head(20)## Alabama Alaska Arizona Arkansas California Colorado
## 4 3 3 4 3 3
## Connecticut Delaware Florida Georgia Hawaii Idaho
## 2 2 3 4 2 1
## Illinois Indiana Iowa Kansas Kentucky Louisiana
## 3 2 1 2 1 4
## Maine Maryland
## 1 3# 可視化
km.res %>% fviz_cluster(USArrests)
4：クラスタ検証統計
＜ポイント＞
・シルエットはクラスタ内の他のオブジェクトと近隣クラスターの類似度を測定する
--- 値は-1から+1の間で表示される
・+1に近い場合は、要素が適切にクラスタ化されていることを示す
--- その要素はグループ内の他のオブジェクトと類似している
・-1に近い場合は、要素のクラスタへの割り当てが不適切であることを示唆する
--- 割当先を変えることで改善する可能性あり
・今回はクラスタ３で負のシルエット値を有するサンプルがあることが分かる# クラスタ情報に基づいてシルエット分析
sil <- km.res$cluster %>% silhouette(dist(df))
# 結果表示
sil %>%
set_rownames(rownames(USArrests)) %>%
.[,1:3] %>%
head()## cluster neighbor sil_width
## Alabama 4 3 0.48577530
## Alaska 3 4 0.05825209
## Arizona 3 2 0.41548326
## Arkansas 4 2 0.11870947
## California 3 2 0.43555885
## Colorado 3 2 0.32654235# 可視化
sil %>% fviz_silhouette()## cluster size ave.sil.width
## 1 1 13 0.37
## 2 2 16 0.34
## 3 3 13 0.27
## 4 4 8 0.39
クラスタ分析の強化
はじめに
＜ポイント＞
・factoextra::eclust()は標準パッケージの関数と比較していくつかのメリットを提供
・メリットは分析ステップの2-4を一括で行うことを意味する
＜メリット＞
・クラスタ分析のワークフローを簡素化
・階層クラスタリングとパーティション化クラスタリングの両方で使用できる
・ギャップ統計を自動的に計算して適切な数のクラスタを推定
・シルエット情報を自動的に提供
・ggplot2を使って美しいグラフを描く
＜分析ステップ＞
・1：クラスタ性の評価
・2：クラスター数を見積り
・3：meansの計算
・4：クラスタ検証統計データ準備
＜ポイント＞
・USArrestsデータセットを使用# データロード
data(USArrests)
# データサイズ
USArrests %>% dim()## [1] 50 4# 基準化
df <- USArrests %>% scale()
df %>% head()## Murder Assault UrbanPop Rape
## Alabama 1.24256408 0.7828393 -0.5209066 -0.003416473
## Alaska 0.50786248 1.1068225 -1.2117642 2.484202941
## Arizona 0.07163341 1.4788032 0.9989801 1.042878388
## Arkansas 0.23234938 0.2308680 -1.0735927 -0.184916602
## California 0.27826823 1.2628144 1.7589234 2.067820292
## Colorado 0.02571456 0.3988593 0.8608085 1.864967207k平均法
# k平均法
# --- プロット表示あり
res.km <- df %>% eclust(FUNcluster = "kmeans", graph = FALSE)
res.km %>% print()## K-means clustering with 3 clusters of sizes 13, 29, 8
##
## Cluster means:
## Murder Assault UrbanPop Rape
## 1 0.6950701 1.0394414 0.72263703 1.27693964
## 2 -0.7010700 -0.7071522 -0.09924526 -0.57773737
## 3 1.4118898 0.8743346 -0.81452109 0.01927104
##
## Clustering vector:
## Alabama Alaska Arizona Arkansas California
## 3 1 1 3 1
## Colorado Connecticut Delaware Florida Georgia
## 1 2 2 1 3
## Hawaii Idaho Illinois Indiana Iowa
## 2 2 1 2 2
## Kansas Kentucky Louisiana Maine Maryland
## 2 2 3 2 1
## Massachusetts Michigan Minnesota Mississippi Missouri
## 2 1 2 3 1
## Montana Nebraska Nevada New Hampshire New Jersey
## 2 2 1 2 2
## New Mexico New York North Carolina North Dakota Ohio
## 1 1 3 2 2
## Oklahoma Oregon Pennsylvania Rhode Island South Carolina
## 2 2 2 2 3
## South Dakota Tennessee Texas Utah Vermont
## 2 3 1 2 2
## Virginia Washington West Virginia Wisconsin Wyoming
## 2 2 2 2 2
##
## Within cluster sum of squares by cluster:
## [1] 19.922437 53.354791 8.316061
## (between_SS / total_SS = 58.4 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss"
## [5] "tot.withinss" "betweenss" "size" "iter"
## [9] "ifault" "silinfo" "nbclust" "data"
## [13] "gap_stat"# 最適クラスタ数
res.km$nbclust## [1] 3# Step2：ギャップ統計量
res.km$gap_stat %>% fviz_gap_stat()
# Step3：クラスタ結果の表示
res.km %>% fviz_cluster()
# Step4：シルエット分析
res.km %>% fviz_silhouette()## cluster size ave.sil.width
## 1 1 13 0.31
## 2 2 29 0.38
## 3 3 8 0.39
階層クラスタリング
# 階層クラスター分析
# --- プロット表示あり
res.hc <- df %>% eclust("hclust")
res.hc %>% print()##
## Call:
## stats::hclust(d = x, method = hc_method)
##
## Cluster method : ward.D2
## Distance : euclidean
## Number of objects: 50# Step2：ギャップ統計量
res.hc$gap_stat %>% fviz_gap_stat()
# Step3：クラスタ結果の表示
# --- デンドログラム
res.hc %>% fviz_dend(rect = TRUE)
# Step3：クラスタ結果の表示
# --- クラスタープロット
res.hc %>% fviz_cluster()
# Step4：シルエット分析
res.hc %>% fviz_silhouette()## cluster size ave.sil.width
## 1 1 19 0.26
## 2 2 19 0.28
## 3 3 12 0.43
クイックスタート
準備
はじめに
＜ポイント＞
・多様なクラスタリング手法を概観する
＜クラスタリング手法＞
・パーティショニング型クラスタリング
・階層的クラスタリング
・ファジィクラスタリング
・密度ベースのクラスタリング
・モデルベースのクラスタリングデータ準備
# データロード
data("USArrests")
# 欠損値除去後に基準化
my_data <- USArrests %>% na.omit() %>%scale()
# データ確認
my_data %>% head()## Murder Assault UrbanPop Rape
## Alabama 1.24256408 0.7828393 -0.5209066 -0.003416473
## Alaska 0.50786248 1.1068225 -1.2117642 2.484202941
## Arizona 0.07163341 1.4788032 0.9989801 1.042878388
## Arkansas 0.23234938 0.2308680 -1.0735927 -0.184916602
## California 0.27826823 1.2628144 1.7589234 2.067820292
## Colorado 0.02571456 0.3988593 0.8608085 1.864967207距離測定
＜ポイント＞
get_dist()：データ行列の行間の距離行列を計算
--- dist()と比較して、相関ベースの距離測定をサポート
fviz_dist()：距離行列を視覚化# 距離行列の作成
res.dist <- USArrests %>% get_dist(stand = TRUE, method = "pearson")
res.dist %>% as.matrix() %>% .[1:5, 1:5] %>% round(2)## Alabama Alaska Arizona Arkansas California
## Alabama 0.00 0.71 1.45 0.09 1.87
## Alaska 0.71 0.00 0.83 0.37 0.81
## Arizona 1.45 0.83 0.00 1.18 0.29
## Arkansas 0.09 0.37 1.18 0.00 1.59
## California 1.87 0.81 0.29 1.59 0.00# データ属性
res.dist %>% attributes()## $Labels
## [1] "Alabama" "Alaska" "Arizona" "Arkansas"
## [5] "California" "Colorado" "Connecticut" "Delaware"
## [9] "Florida" "Georgia" "Hawaii" "Idaho"
## [13] "Illinois" "Indiana" "Iowa" "Kansas"
## [17] "Kentucky" "Louisiana" "Maine" "Maryland"
## [21] "Massachusetts" "Michigan" "Minnesota" "Mississippi"
## [25] "Missouri" "Montana" "Nebraska" "Nevada"
## [29] "New Hampshire" "New Jersey" "New Mexico" "New York"
## [33] "North Carolina" "North Dakota" "Ohio" "Oklahoma"
## [37] "Oregon" "Pennsylvania" "Rhode Island" "South Carolina"
## [41] "South Dakota" "Tennessee" "Texas" "Utah"
## [45] "Vermont" "Virginia" "Washington" "West Virginia"
## [49] "Wisconsin" "Wyoming"
##
## $Size
## [1] 50
##
## $call
## as.dist.default(m = 1 - res.cor)
##
## $class
## [1] "dist"
##
## $Diag
## [1] FALSE
##
## $Upper
## [1] FALSEres.dist %>%
fviz_dist(gradient = list(low = "#00AFBB", mid = "white", high = "#FC4E07"))
分割型クラスタリング
クラスタ数の決定
＜ポイント＞
・データセットをkグループのセットに細分するクラスタリングテクニック
・kの数はは分析者があらかじめ指定しなければならない
・GAP統計量がクラスタ数を決定する目安となる
--- 今回の場合は推奨クラスタ数は３# ギャップ統計量を可視化
my_data %>% fviz_nbclust(kmeans, method = "gap_stat")
k-meansクラスタリング
＜ポイント＞
・GAP統計量に基づいてクラスタ数を３に設定
・PC1とPC2に基づいてマップ化# 乱数設定
set.seed(123)
# クラスタリング
km.res <- my_data %>% kmeans(3, nstart = 25)
km.res %>% attributes()## $names
## [1] "cluster" "centers" "totss" "withinss"
## [5] "tot.withinss" "betweenss" "size" "iter"
## [9] "ifault"
##
## $class
## [1] "kmeans"# 可視化
km.res %>%
fviz_cluster(data = my_data,
ellipse.type = "convex",
palette = "jco",
ggtheme = theme_minimal())
pamクラスタリング
＜ポイント＞
・GAP統計量に基づいてクラスタ数を３に設定
・PC1とPC2に基づいてマップ化# クラスタリング
pam.res <- my_data %>% pam(3)
pam.res %>% attributes()## $names
## [1] "medoids" "id.med" "clustering" "objective" "isolation"
## [6] "clusinfo" "silinfo" "diss" "call" "data"
##
## $class
## [1] "pam" "partition"# 可視化
pam.res %>% fviz_cluster()
階層型クラスタリング
クラスタリング
＜ポイント＞
・階層的クラスタリングは生成するクラスタ数をあらかじめ指定する必要がない# 階層的クラスタリング
res.hc <-
USArrests %>%
scale() %>%
dist(method = "euclidean") %>%
hclust(method = "ward.D2") %T>%
print()##
## Call:
## hclust(d = ., method = "ward.D2")
##
## Cluster method : ward.D2
## Distance : euclidean
## Number of objects: 50# 属性確認
res.hc %>% attributes()## $names
## [1] "merge" "height" "order" "labels" "method"
## [6] "call" "dist.method"
##
## $class
## [1] "hclust"デンドログラム
＜ポイント＞
・結果はツリーベースの表現でデンドログラムとよばれる# デンドログラムの作成
res.hc %>%
fviz_dend(k = 4,
cex = 0.5,
k_colors = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"),
color_labels_by_k = TRUE,
rect = TRUE
)
クラスタリングの検証と評価
クラスタリング傾向の評価
＜ポイント＞
・ホプキンスの統計および視覚的アプローチによりクラスタリング傾向を評価できる
＜ホプキンス統計＞
・統計値が1にづくほど大幅にクラスタ化可能であることを示す
＜視覚的アプローチ＞
・順序付けされた相違度画像における対角線に沿った四角形の有色ブロックの数を数えることによってクラスタリング傾向を検出する。#
iris[, -5] %>%
scale() %>%
get_clust_tendency(n = 50,
gradient = list(low = "steelblue",
high = "white"))## $hopkins_stat
## [1] 0.2002686
##
## $plot
最適クラスタ数の決定
＜ポイント＞
・NbClustでは、最良のクラスタ数を決定するための30個のインデックスを提供される# 乱数設定
set.seed(123)
# 最適クラスタ数の決定
res.nbclust <-
USArrests %>%
scale() %>%
NbClust::NbClust(distance = "euclidean",
min.nc = 2,
max.nc = 10,
method = "complete",
index ="all") 
## *** : 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:
## * 9 proposed 2 as the best number of clusters
## * 4 proposed 3 as the best number of clusters
## * 6 proposed 4 as the best number of clusters
## * 2 proposed 5 as the best number of clusters
## * 1 proposed 8 as the best number of clusters
## * 1 proposed 10 as the best number of clusters
##
## ***** Conclusion *****
##
## * According to the majority rule, the best number of clusters is 2
##
##
## *******************************************************************クラスタリング検証統計
# 乱数設定
set.seed(123)
# 階層クラスタリング(k=3でカット)
res.hc <-
iris[, -5] %>%
scale() %>%
eclust("hclust", k = 3, graph = FALSE)
# デンドログラムによる可視化
res.hc %>%
fviz_dend(palette = "jco",
rect = TRUE,
show_labels = FALSE)
# シルエット検査
res.hc %>% fviz_silhouette()## cluster size ave.sil.width
## 1 1 49 0.63
## 2 2 30 0.44
## 3 3 71 0.32
# 負のシルエットを持つサンプル
res.hc$silinfo$widths[, 1:3] %>%
as.data.frame() %>%
arrange(sil_width) %>%
dplyr::filter(sil_width < 0)## cluster neighbor sil_width
## 1 3 2 -0.23036371
## 2 3 2 -0.16107155
## 3 3 2 -0.14761137
## 4 3 2 -0.10091884
## 5 3 2 -0.05302402
## 6 3 2 -0.04756835
## 7 3 2 -0.01789603
## 8 3 2 -0.01269799高度なクラスタリング手法
ハイブリッド
＜ポイント＞
・階層的K平均クラスタリング：k平均結果を改善するハイブリッドアプローチ
・HCPC：主成分の階層的クラスタリングファジィ
＜ポイント＞
・ファジィクラスタリングとは、ソフトメソッドとも呼ばれる
・ファジークラスタリング、アイテムが複数のクラスタのメンバーであることができます
・ファジィc-平均の方法は、最も人気のあるファジークラスタリングアルゴリズムモデルベース
密度ベース
分割型クラスタリング
k平均法
概要
＜ポイント＞
・各観測値はk個のクラスタのどれか1つに所属する
・クラスタの個数kはならかの方法に基づいて予め指定する必要がある
・計算にはユーグリッド距離が使われている
--- 2乗したユーグリッド距離を最小化
＜アルゴリズム＞
1. ランダムにk個のセントロイド(平均値を多次元に拡張したもの)をランダムに決定
2. 各観測値とセントロイドのユーグリッド距離に基づいてクラスタ分類
3. 分類後の各クラスタから新たなセントロイドを計算して再分類
4. 再分類がなくなるまで2と3を繰り返す
＜メリット＞
・アルゴリズムが単純なので高速処理(大規模処理)が可能
＜デメリット＞
・クラスタ数を分析者が指定する必要がある
・異常値に敏感
・結果がランダムに決まる初期値や、データの並び順に依存する
＜代替手法＞
・PAMクラスタリング(メドイド法)
--- 異常値に強いという点でk-meansを代替データ
# データロード
data("USArrests")
# データ基準化
df <- USArrests %>% scale()
# 確認
df %>% head()## Murder Assault UrbanPop Rape
## Alabama 1.24256408 0.7828393 -0.5209066 -0.003416473
## Alaska 0.50786248 1.1068225 -1.2117642 2.484202941
## Arizona 0.07163341 1.4788032 0.9989801 1.042878388
## Arkansas 0.23234938 0.2308680 -1.0735927 -0.184916602
## California 0.27826823 1.2628144 1.7589234 2.067820292
## Colorado 0.02571456 0.3988593 0.8608085 1.864967207関数
＜構文＞
stats::kmeans(x, centers, iter.max = 10, nstart = 1)
＜引数＞
・x ：数値形式の行列又はデータフレーム
・centers ：クラスター数
・iter.max：許容される反復の最大回数(デフォルト値は10)。
・nstart ：ランダム開始パーティションの数(デフォルトは1)
クラスタ数の見積もり
＜ポイント＞
・k-meansでは分析者がクラスタ数を決定しなければならない
・ここではクラスタ内平方和(分散)を小さくするという観点で見積もる(wss)
--- fviz_nbclust()は silhouette / wss / gap_stat が選択可能
--- 平方和は４以降は効率的に逓減しなくなる
＜参考＞
・fcp::pamk()のようにクラスタ数を関数内で自動決定するものもある
--- シルエット分析などに基づいて決定
df %>%
fviz_nbclust(kmeans, method = "wss") +
geom_vline(xintercept = 4, linetype = 2)
クラスタリングの計算
＜ポイント＞
・k-meansはk個のランダムに選択された重心から始まる
--- 乱数シードの設定することで再現性が担保できる
・k-meansはランダムな開始割り当てに敏感であるため、nstart=25と指定
--- 25のランダムな開始割当てを試行、クラスタの変動が最も小さいものを選択# 乱数設定
set.seed(123)
# クラスタリング
# --- 見積もりで決定した４を使用
km.res <- df %>% kmeans(centers = 4, nstart = 25)
km.res %>% print()## K-means clustering with 4 clusters of sizes 13, 16, 13, 8
##
## Cluster means:
## Murder Assault UrbanPop Rape
## 1 -0.9615407 -1.1066010 -0.9301069 -0.96676331
## 2 -0.4894375 -0.3826001 0.5758298 -0.26165379
## 3 0.6950701 1.0394414 0.7226370 1.27693964
## 4 1.4118898 0.8743346 -0.8145211 0.01927104
##
## Clustering vector:
## Alabama Alaska Arizona Arkansas California
## 4 3 3 4 3
## Colorado Connecticut Delaware Florida Georgia
## 3 2 2 3 4
## Hawaii Idaho Illinois Indiana Iowa
## 2 1 3 2 1
## Kansas Kentucky Louisiana Maine Maryland
## 2 1 4 1 3
## Massachusetts Michigan Minnesota Mississippi Missouri
## 2 3 1 4 3
## Montana Nebraska Nevada New Hampshire New Jersey
## 1 1 3 1 2
## New Mexico New York North Carolina North Dakota Ohio
## 3 3 4 1 2
## Oklahoma Oregon Pennsylvania Rhode Island South Carolina
## 2 2 2 2 4
## South Dakota Tennessee Texas Utah Vermont
## 1 4 3 2 1
## Virginia Washington West Virginia Wisconsin Wyoming
## 2 2 1 1 2
##
## Within cluster sum of squares by cluster:
## [1] 11.952463 16.212213 19.922437 8.316061
## (between_SS / total_SS = 71.2 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss"
## [5] "tot.withinss" "betweenss" "size" "iter"
## [9] "ifault"関数結果
＜出力項目＞
・cluster : 各要素が割り当てられるクラスタ
・centers : クラスターセンターのマトリックス
・totss :
・withinss : クラスター内合計のベクトル(クラスターごとに1つのコンポーネント)
・tot.withinss : withinssの平方和
・betweenss :
・size : 各クラスタの観測数# クラスタの属性情報
km.res %>% attributes()## $names
## [1] "cluster" "centers" "totss" "withinss"
## [5] "tot.withinss" "betweenss" "size" "iter"
## [9] "ifault"
##
## $class
## [1] "kmeans"# 所属クラスタ
km.res$cluster %>% head()## Alabama Alaska Arizona Arkansas California Colorado
## 4 3 3 4 3 3# クラスタごとの個数
km.res$size## [1] 13 16 13 8# クラスタ平均
km.res$centers## Murder Assault UrbanPop Rape
## 1 -0.9615407 -1.1066010 -0.9301069 -0.96676331
## 2 -0.4894375 -0.3826001 0.5758298 -0.26165379
## 3 0.6950701 1.0394414 0.7226370 1.27693964
## 4 1.4118898 0.8743346 -0.8145211 0.01927104クラスタ集計
＜ポイント＞
・k-meansによるクラスタ情報は.$clusterに格納される
・クラスタごとに４つの変数の平均値を確認(4次元)# 元のデータセットにクラスタ情報を追加
dd <-
USArrests %>%
bind_cols(cluster = km.res$cluster) %>%
group_by(cluster) %T>%
print()## # A tibble: 50 x 5
## # Groups: cluster [4]
## Murder Assault UrbanPop Rape cluster
## <dbl> <int> <int> <dbl> <int>
## 1 13.2 236 58 21.2 4
## 2 10.0 263 48 44.5 3
## 3 8.10 294 80 31.0 3
## 4 8.80 190 50 19.5 4
## 5 9.00 276 91 40.6 3
## 6 7.90 204 78 38.7 3
## 7 3.30 110 77 11.1 2
## 8 5.90 238 72 15.8 2
## 9 15.4 335 80 31.9 3
## 10 17.4 211 60 25.8 4
## # ... with 40 more rows# 各クラスタの変数の平均
dd %>% summarise_all(funs(mean))## # A tibble: 4 x 5
## cluster Murder Assault UrbanPop Rape
## <int> <dbl> <dbl> <dbl> <dbl>
## 1 1 3.60 78.5 52.1 12.2
## 2 2 5.66 139 73.9 18.8
## 3 3 10.8 257 76.0 33.2
## 4 4 13.9 244 53.8 21.4可視化
＜ポイント＞
・クラスタ分析はプロットで確認することで非常に理解が深まる
--- 通常だと変数分の次元数(今回は4次元)となる
・解決策としてPCAによる次元削減により2次元に変換してプロットするkm.res %>%
fviz_cluster(data = df,
palette = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"),
ellipse.type = "euclid",
star.plot = TRUE,
repel = TRUE,
ggtheme = theme_minimal()
)
kメドイド法
概要
＜ポイント＞
・非類似度行列を用いて非類似度自体を最小化する
--- 異常値に対してk平均法よりも頑健
--- 非類似度を最小化する点をメドイドと呼ぶ
＜アルゴリズム＞
1. 非類似度行列に基づいてランダムにk個のメドイドに決定
2. 各観測値とメドイドとの非類似度に基づいてクラスタ分類
3. 分類後の各クラスタから新たなメドイドを計算して再分類
4. 再分類がなくなるまで2と3を繰り返す
＜メリット＞
・異常値に対してk-meansよりも頑健
＜デメリット＞
・大規模データの処理には不向き
--- 非類似度行列は、観測点数の2乗で拡大していくデータ
# データロード
data("USArrests")
# データ基準化
df <- USArrests %>% scale()
# 確認
df %>% head()## Murder Assault UrbanPop Rape
## Alabama 1.24256408 0.7828393 -0.5209066 -0.003416473
## Alaska 0.50786248 1.1068225 -1.2117642 2.484202941
## Arizona 0.07163341 1.4788032 0.9989801 1.042878388
## Arkansas 0.23234938 0.2308680 -1.0735927 -0.184916602
## California 0.27826823 1.2628144 1.7589234 2.067820292
## Colorado 0.02571456 0.3988593 0.8608085 1.864967207関数
＜構文＞
cluster::pam(x, k, metric = "euclidean", stand = FALSE)
＜引数＞
・x ：数値形式の行列又はデータフレーム
：daisy()又はdist()から出力した非類似度行列
・k ：クラスター数
・metric ：使用される距離尺度(euclidean / manhattan)
・stand ：TRUEの場合は非類似度を計算する前に基準化を実施
クラスタ数の見積もり
＜ポイント＞
・k-meansと同様に分析者がクラスタ数を決定しなければならない
・ここでは平均シルエット法で見積もるdf %>%
fviz_nbclust(kmeans, method = "silhouette") +
theme_classic()
クラスタリングの計算
＜ポイント＞
・k-meansはk個のランダムに選択された重心から始まる
--- 乱数シードの設定することで再現性が担保できる
・k-meansはランダムな開始割り当てに敏感であるため、nstart=25と指定
--- 25のランダムな開始割当てを試行、クラスタの変動が最も小さいものを選択# クラスタリング
# --- 見積もりで決定した２を使用
pam.res <- df %>% pam(k = 2)
pam.res %>% print()## Medoids:
## ID Murder Assault UrbanPop Rape
## New Mexico 31 0.8292944 1.3708088 0.3081225 1.1603196
## Nebraska 27 -0.8008247 -0.8250772 -0.2445636 -0.5052109
## Clustering vector:
## 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 2 1 2 2
## Kansas Kentucky Louisiana Maine Maryland
## 2 2 1 2 1
## Massachusetts Michigan Minnesota Mississippi Missouri
## 2 1 2 1 1
## Montana Nebraska Nevada New Hampshire New Jersey
## 2 2 1 2 2
## New Mexico New York North Carolina North Dakota Ohio
## 1 1 1 2 2
## Oklahoma Oregon Pennsylvania Rhode Island South Carolina
## 2 2 2 2 1
## South Dakota Tennessee Texas Utah Vermont
## 2 1 1 2 2
## Virginia Washington West Virginia Wisconsin Wyoming
## 2 2 2 2 2
## Objective function:
## build swap
## 1.441358 1.368969
##
## Available components:
## [1] "medoids" "id.med" "clustering" "objective" "isolation"
## [6] "clusinfo" "silinfo" "diss" "call" "data"関数結果
＜出力項目＞
・cluster : 各要素が割り当てられるクラスタ
・centers : クラスターセンターのマトリックス
・totss :
・withinss : クラスター内合計のベクトル(クラスターごとに1つのコンポーネント)
・tot.withinss : withinssの平方和
・betweenss :
・size : 各クラスタの観測数# クラスタの属性情報
pam.res %>% attributes()## $names
## [1] "medoids" "id.med" "clustering" "objective" "isolation"
## [6] "clusinfo" "silinfo" "diss" "call" "data"
##
## $class
## [1] "pam" "partition"# 所属クラスタ
pam.res$clustering %>% head()## Alabama Alaska Arizona Arkansas California Colorado
## 1 1 1 2 1 1# クラスタごとの個数
pam.res$id.med## [1] 31 27# メドイド
pam.res$medoids## Murder Assault UrbanPop Rape
## New Mexico 0.8292944 1.3708088 0.3081225 1.1603196
## Nebraska -0.8008247 -0.8250772 -0.2445636 -0.5052109クラスタ集計
＜ポイント＞
・k-meansによるクラスタ情報は.$clusterに格納される
・クラスタごとに４つの変数の平均値を確認(4次元)# 元のデータセットにクラスタ情報を追加
dd <-
USArrests %>%
bind_cols(cluster = pam.res$cluster) %>%
group_by(cluster) %T>%
print()## # A tibble: 50 x 5
## # Groups: cluster [2]
## Murder Assault UrbanPop Rape cluster
## <dbl> <int> <int> <dbl> <int>
## 1 13.2 236 58 21.2 1
## 2 10.0 263 48 44.5 1
## 3 8.10 294 80 31.0 1
## 4 8.80 190 50 19.5 2
## 5 9.00 276 91 40.6 1
## 6 7.90 204 78 38.7 1
## 7 3.30 110 77 11.1 2
## 8 5.90 238 72 15.8 2
## 9 15.4 335 80 31.9 1
## 10 17.4 211 60 25.8 1
## # ... with 40 more rows# 各クラスタの変数の平均
dd %>% summarise_all(funs(mean))## # A tibble: 2 x 5
## cluster Murder Assault UrbanPop Rape
## <int> <dbl> <dbl> <dbl> <dbl>
## 1 1 12.2 255 68.4 29.2
## 2 2 4.87 114 63.6 15.9可視化
＜ポイント＞
・クラスタ分析はプロットで確認することで非常に理解が深まる
--- 通常だと変数分の次元数(今回は4次元)となる
・解決策としてPCAによる次元削減により2次元に変換してプロットするpam.res %>%
fviz_cluster(palette = c("#00AFBB", "#FC4E07"),
ellipse.type = "t",
repel = TRUE,
ggtheme = theme_classic()
)
大規模クラスタリング
概要
＜ポイント＞
・数千件以上の観測値をもつデータに対してメドイド法を適用
--- cluster::pam()が持つ計算時間と大規模行列のメモリ問題を軽減
・データセット全体ではなく、データを分割してサンプルごとの最適メドイドを生成
・サンプリングバイアスを最小限に抑えるために、予め指定された回数繰り返します
＜アルゴリズム＞
1. データセットをランダムに分割
＜メリット＞
・大規模データに対してメドイド法を適用できる
＜デメリット＞
・データ
＜ポイント＞
・今回は乱数を用いてデータセットを作成
・２列のデータセットを２種類作成して列結合# 乱数設定
set.seed(1234)
# データ基準化
df <-
rbind(cbind(rnorm(200, 0, 8),
rnorm(200, 0, 8)),
cbind(rnorm(300, 50, 8),
rnorm(300, 50, 8))) %>%
set_colnames(c("x", "y")) %>%
set_rownames(paste0("S", 1:nrow(.)))
# サイズ
df %>% dim()## [1] 500 2# 確認
df %>% head()## x y
## S1 -9.656526 3.881815
## S2 2.219434 5.574150
## S3 8.675529 1.484111
## S4 -18.765582 5.605868
## S5 3.432998 2.493448
## S6 4.048447 6.083699clara関数
＜構文＞
cluster::clara(x, k, metric = "euclidean", stand = FALSE,
samples = 5, pamLike = FALSE)
＜引数＞
・x ：数値データ行列またはデータフレーム(欠損値は許容)
・k ：クラスタ数
・metric ：距離行列の測定方法(euclidean / mahattan)
・stand ：TRUEの場合、非類似度を計算する前に基準化を実行
・samples ：データセットから抽出されるサンプルの数
・pamLike ：クラスタ数の見積もり
＜ポイント＞
・k-meansと同様に分析者がクラスタ数を決定しなければならない
・ここでは平均シルエット法で見積もる
・プロットから推奨クラスター数は２と分かるdf %>%
fviz_nbclust(clara, method = "silhouette") +
theme_classic()
クラスタリングの計算
＜表示項目＞
・medoids ：クラスタを表すオブジェクト
・clustering ：各オブジェクトのクラスタ番号を含むベクトル
・sample ：ベストサンプルの観測値のラベル# クラスタリング
# --- 見積もりで決定した２を使用
clara.res <- df %>% clara(2, samples = 50, pamLike = TRUE) %T>% print()## Call: clara(x = ., k = 2, samples = 50, pamLike = TRUE)
## Medoids:
## x y
## S121 -1.531137 1.145057
## S455 48.357304 50.233499
## Objective function: 9.87862
## Clustering vector: Named int [1:500] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...
## - attr(*, "names")= chr [1:500] "S1" "S2" "S3" "S4" "S5" "S6" "S7" ...
## Cluster sizes: 200 300
## Best sample:
## [1] S37 S49 S54 S63 S68 S71 S76 S80 S82 S101 S103 S108 S109 S118
## [15] S121 S128 S132 S138 S144 S162 S203 S210 S216 S231 S234 S249 S260 S261
## [29] S286 S299 S304 S305 S312 S315 S322 S350 S403 S450 S454 S455 S456 S465
## [43] S488 S497
##
## Available components:
## [1] "sample" "medoids" "i.med" "clustering" "objective"
## [6] "clusinfo" "diss" "call" "silinfo" "data"関数結果
＜出力項目＞
・cluster : 各要素が割り当てられるクラスタ
・centers : クラスターセンターのマトリックス
・totss :
・withinss : クラスター内合計のベクトル(クラスターごとに1つのコンポーネント)
・tot.withinss : withinssの平方和
・betweenss :
・size : 各クラスタの観測数# クラスタの属性情報
clara.res %>% attributes()## $names
## [1] "sample" "medoids" "i.med" "clustering" "objective"
## [6] "clusinfo" "diss" "call" "silinfo" "data"
##
## $class
## [1] "clara" "partition"# 所属クラスタ
clara.res$clustering %>% head()## S1 S2 S3 S4 S5 S6
## 1 1 1 1 1 1# メドイド
clara.res$medoids## x y
## S121 -1.531137 1.145057
## S455 48.357304 50.233499クラスタ集計
＜ポイント＞
・元のデータセットに作成したクラスター情報を結合# 元のデータセットにクラスタ情報を追加
dd <- df %>% cbind(cluster = clara.res$cluster)
dd %>% head()## x y cluster
## S1 -9.656526 3.881815 1
## S2 2.219434 5.574150 1
## S3 8.675529 1.484111 1
## S4 -18.765582 5.605868 1
## S5 3.432998 2.493448 1
## S6 4.048447 6.083699 1# クラスタごとに集計
dd %>% as_tibble() %>% group_by(cluster) %>% tally()## # A tibble: 2 x 2
## cluster n
## <dbl> <int>
## 1 1.00 200
## 2 2.00 300可視化
＜ポイント＞
・クラスタ分析はプロットで確認することで非常に理解が深まる
--- 通常だと変数分の次元数(今回は4次元)となる
・解決策としてPCAによる次元削減により2次元に変換してプロットするclara.res %>%
fviz_cluster(palette = c("#00AFBB", "#FC4E07"),
ellipse.type = "t",
geom = "point",
pointsize = 1,
ggtheme = theme_classic()
)
階層型クラスタリング
凝集型階層クラスタ
概要
＜ポイント＞
・ボトムアップ形式で階層型クラスタリングを作成する
・階層型クラスタリングの最も一般的なタイプ
＜アルゴリズム＞
1. 個別の観測点を1つのクラスタ(リーフ)とみなす
2. 反復ステップにおいて最も非類似度が小さいクラスタを統合
3. 全ての観測点が１つのクラスタに入るまで繰り返すデータ
# データロード
data("USArrests")
# データ基準化
df <- USArrests %>% scale()
# データサイズ
df %>% dim()## [1] 50 4# データ確認
df %>% head()## Murder Assault UrbanPop Rape
## Alabama 1.24256408 0.7828393 -0.5209066 -0.003416473
## Alaska 0.50786248 1.1068225 -1.2117642 2.484202941
## Arizona 0.07163341 1.4788032 0.9989801 1.042878388
## Arkansas 0.23234938 0.2308680 -1.0735927 -0.184916602
## California 0.27826823 1.2628144 1.7589234 2.067820292
## Colorado 0.02571456 0.3988593 0.8608085 1.864967207類似性指標
＜ポイント＞
・類似性指標は距離行列や非類似度などとも呼ばれる
・以下のような関数で測定するのが一般的
--- stats::dist()
--- factoextra::get_dist()
--- cluster::daisy()# 距離行列の作成
res.dist <- df %>% dist(method = "euclidean")
# 通常の行列に変換(見やすくするため)
res.dist %>% as.matrix() %>% .[1:6, 1:6]## Alabama Alaska Arizona Arkansas California Colorado
## Alabama 0.000000 2.703754 2.293520 1.289810 3.263110 2.651067
## Alaska 2.703754 0.000000 2.700643 2.826039 3.012541 2.326519
## Arizona 2.293520 2.700643 0.000000 2.717758 1.310484 1.365031
## Arkansas 1.289810 2.826039 2.717758 0.000000 3.763641 2.831051
## California 3.263110 3.012541 1.310484 3.763641 0.000000 1.287619
## Colorado 2.651067 2.326519 1.365031 2.831051 1.287619 0.000000計算
＜ポイント＞
・dist()などで作成した距離行列を受けて観測点からクラスタを作成していく
＜結合方法＞
・ward.D ：ウォード法
・ward.D2 ：ウォード法
・single ：
・complete ：
・average ：
・mcquitty ：
・median ：
・centroid ：# ststs::hclust()
res.hc <- res.dist %>% hclust(method = "ward.D2")
res.hc %>% print()##
## Call:
## hclust(d = ., method = "ward.D2")
##
## Cluster method : ward.D2
## Distance : euclidean
## Number of objects: 50デンドログラム
＜ポイント＞
・デンドログラム(樹形図)は階層型クラスタリングをグラフィカルに表示する方法
・plot()でも表示できるが、fviz_dend()を使用するとより美しいグラフィックが可能
・高さ(Y軸)はオブジェクト同士の距離を示すres.hc %>% fviz_dend(cex = 0.5)
クラスタツリーの確認
＜ポイント＞
・ツリー内の距離(高さ)が元の距離を正確に反映しているかどうかを評価する
・共役距離とdist()関数によって生成された元の距離データとの間の相関を計算# コーフェン行列の算出
res.coph <- res.hc %>% cophenetic()
res.coph %>% as.matrix() %>% dim()## [1] 50 50# 距離行列
res.dist %>% as.matrix() %>% dim()## [1] 50 50# 距離行列とコーフェン行列の相関
cor(res.dist, res.coph)## [1] 0.6975266# 平均リンケージ方法を使用して再度クラスタを作成
res.hc2 <- res.dist %>% hclust(method = "average")
cor(res.dist, cophenetic(res.hc2))## [1] 0.7180382樹形図のカット
＜ポイント＞
・データをクラスタに分割するため階層ツリーを所定の高さで切断できる
・カット結果はデンドログラムや散布図に反映することも可能# ツリーが４グループになるようにカット
grp <- res.hc %>% cutree(k = 4)
grp %>% head(n = 4)## Alabama Alaska Arizona Arkansas
## 1 2 2 3# テーブル化
grp %>% table()## .
## 1 2 3 4
## 7 12 19 12# クラスタのメンバーの取得
df %>% rownames() %>% .[grp == 1]## [1] "Alabama" "Georgia" "Louisiana" "Mississippi"
## [5] "North Carolina" "South Carolina" "Tennessee"# カット結果をデンドログラムに反映
res.hc %>%
fviz_dend(k = 4,
cex = 0.5,
k_colors = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"),
color_labels_by_k = TRUE,
rect = TRUE
)
# カット結果を散布図に反映
list(data = df, cluster = grp) %>%
fviz_cluster(palette = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"),
ellipse.type = "convex",
repel = TRUE,
show.clust.cent = FALSE,
ggtheme = theme_minimal())
参考：agnes関数
＜ポイント＞
・凝集型クラスタリングはhclust()が有名だが、{cluster}でも関数が提供されている
・cluster::agnes()を使用するとクラスタ分析の以下の操作を一括して行える
--- scale() + dist() + hclust()# {cluster}による凝集型クラスタリング
res.agnes <-
USArrests %>%
cluster::agnes(stand = TRUE, metric = "euclidean", method = "ward") %T>%
print()## Call: cluster::agnes(x = ., metric = "euclidean", stand = TRUE, method = "ward")
## Agglomerative coefficient: 0.934098
## Order of objects:
## [1] Alabama Louisiana Georgia Tennessee
## [5] Mississippi South Carolina North Carolina Alaska
## [9] California Nevada Colorado Arizona
## [13] Maryland New Mexico Michigan Florida
## [17] Illinois New York Texas Arkansas
## [21] Kentucky Virginia Wyoming Indiana
## [25] Kansas Oklahoma Ohio Pennsylvania
## [29] Missouri Oregon Washington Connecticut
## [33] Rhode Island Delaware Hawaii Utah
## [37] Massachusetts New Jersey Idaho Montana
## [41] Nebraska Iowa New Hampshire Maine
## [45] Minnesota Wisconsin North Dakota Vermont
## [49] South Dakota West Virginia
## Height (summary):
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.2620 0.9627 1.3330 2.2030 2.2250 16.3200
##
## Available components:
## [1] "order" "height" "ac" "merge" "diss" "call"
## [7] "method" "order.lab" "data"# デンドログラムの作成
res.agnes %>% fviz_dend(cex = 0.6, k = 4)
区分型階層クラスタ
概要
＜ポイント＞
・トップダウン形式で階層型クラスタリングを作成する
＜アルゴリズム＞
1. 全ての観測点を1つの大きなクラスタに含める
2. 反復ステップにおいて最も非類似度が大きいクラスタを２つに分割
3. 全ての観測点が独自のクラスタに入るまで繰り返す計算
＜ポイント＞
・区分型階層クラスタリングはcluster::diana()を使用するres.diana <- USArrests %>% diana(stand = TRUE, metric = "euclidean")
res.diana %>% print()## Merge:
## [,1] [,2]
## [1,] -15 -29
## [2,] -13 -32
## [3,] -23 -49
## [4,] -14 -36
## [5,] -20 -31
## [6,] -16 -38
## [7,] -37 -47
## [8,] 1 -19
## [9,] 4 -35
## [10,] -41 -48
## [11,] -46 -50
## [12,] -12 -27
## [13,] -21 -30
## [14,] -24 -40
## [15,] 2 -43
## [16,] -17 -26
## [17,] -1 -42
## [18,] -10 -18
## [19,] 9 6
## [20,] -4 11
## [21,] -11 -44
## [22,] -7 -39
## [23,] 8 -34
## [24,] 10 -45
## [25,] 5 -22
## [26,] 17 18
## [27,] 14 -33
## [28,] 12 3
## [29,] -5 -28
## [30,] -3 25
## [31,] 29 -6
## [32,] 15 -25
## [33,] 30 32
## [34,] 22 13
## [35,] 23 24
## [36,] 19 7
## [37,] 20 -8
## [38,] 34 21
## [39,] 28 16
## [40,] 37 36
## [41,] 26 27
## [42,] 39 35
## [43,] 33 -9
## [44,] 43 31
## [45,] 40 38
## [46,] -2 44
## [47,] 45 42
## [48,] 41 46
## [49,] 48 47
## Order of objects:
## [1] Alabama Tennessee Georgia Louisiana
## [5] Mississippi South Carolina North Carolina Alaska
## [9] Arizona Maryland New Mexico Michigan
## [13] Illinois New York Texas Missouri
## [17] Florida California Nevada Colorado
## [21] Arkansas Virginia Wyoming Delaware
## [25] Indiana Oklahoma Ohio Kansas
## [29] Pennsylvania Oregon Washington Connecticut
## [33] Rhode Island Massachusetts New Jersey Hawaii
## [37] Utah Idaho Nebraska Minnesota
## [41] Wisconsin Kentucky Montana Iowa
## [45] New Hampshire Maine North Dakota South Dakota
## [49] West Virginia Vermont
## Height:
## [1] 1.0340801 1.3687929 1.0408349 2.8205587 0.9771986 1.3960934 5.4623592
## [8] 4.0837483 1.4720236 0.6743186 1.3385319 1.9655812 0.4336501 1.0044328
## [15] 1.8355779 2.9848751 3.0534495 1.4621581 1.7044355 7.4792599 1.2084105
## [22] 0.8788445 2.1852059 2.5927716 0.6281577 0.8500187 1.0553617 0.6765523
## [29] 2.1808837 0.7311305 3.7810250 1.2903008 1.9665793 0.9626545 2.5032675
## [36] 1.2698293 5.3463933 0.9132211 1.4571402 0.6215418 2.5227753 1.0181020
## [43] 2.9549994 0.2619577 0.7930730 1.2918266 2.1082699 0.8721679 1.3129132
## Divisive coefficient:
## [1] 0.8530481
##
## Available components:
## [1] "order" "height" "dc" "merge" "diss" "call"
## [7] "order.lab" "data"可視化
res.diana %>%
fviz_dend(cex = 0.5,
k = 4,
palette = "jco"
)
樹形図の比較
概要
＜ポイント＞
・階層型クラスタリングで作成した２つの樹形図を比較する
・{dendextend}によりサポートされているデータ準備
＜ポイント＞
・USArrestsデータセットを50行から10行に減らしたものを使う# データロード
data("USArrests")
# 乱数設定
set.seed(123)
# 10レコードのデータセットを作成
ss <- sample(1:50, 10)
df <- USArrests %>% scale() %>% .[ss, ]
df## Murder Assault UrbanPop Rape
## Iowa -1.28297267 -1.37704849 -0.58999237 -1.06038781
## Rhode Island -1.00745957 0.03887798 1.48258036 -1.38068216
## Maryland 0.80633501 1.55079947 0.10086521 0.70123109
## Tennessee 1.24256408 0.20686926 -0.45182086 0.60514278
## Utah -1.05337842 -0.60908837 0.99898006 0.17808366
## Arizona 0.07163341 1.47880321 0.99898006 1.04287839
## Mississippi 1.90838741 1.05882502 -1.48810723 -0.44115208
## Wisconsin -1.19113497 -1.41304662 0.03177945 -1.11377020
## Virginia 0.16347111 -0.17711080 -0.17547783 -0.05679886
## Maine -1.30593210 -1.05306531 -1.00450692 -1.43406455計算
＜ポイント＞
・距離測定はユーグリッド距離を使用
・階層クラスタの形成を平均法とウォード法の２種類を使用
・2種類のクラスタから各デンドログラムを作成してリスト化# 距離行列の計算
res.dist <- dist(df, method = "euclidean")
# クラスタ1：平均法
hc1 <- res.dist %>% hclust(method = "average")
hc1##
## Call:
## hclust(d = ., method = "average")
##
## Cluster method : average
## Distance : euclidean
## Number of objects: 10# クラスタ2：ウォード法
hc2 <- res.dist %>% hclust(method = "ward.D2")
hc2##
## Call:
## hclust(d = ., method = "ward.D2")
##
## Cluster method : ward.D2
## Distance : euclidean
## Number of objects: 10# デンドログラム作成
dend1 <- hc1 %>% as.dendrogram()
dend2 <- hc2 %>% as.dendrogram()
# 2つのデンドログラムをリスト化
dend_list <- dendlist(dend1, dend2)
dend_list## [[1]]
## 'dendrogram' with 2 branches and 10 members total, at height 3.415145
##
## [[2]]
## 'dendrogram' with 2 branches and 10 members total, at height 6.657022
##
## attr(,"class")
## [1] "dendlist"樹形図の比較
＜ポイント＞
・タングレグラムは2つの樹形図を横に並びにしてラベルを線で結ぶ# タングレグラム
tanglegram(dend1, dend2)
# 各種オプションを設定
tanglegram(dend1, dend2,
highlight_distinct_edges = FALSE,
common_subtrees_color_lines = FALSE,
common_subtrees_color_branches = TRUE,
main = paste("entanglement =", round(entanglement(dend_list), 2))
)
樹形図の相関比較
# コーフェン行列
dend_list %>% cor.dendlist(method = "cophenetic")## [,1] [,2]
## [1,] 1.0000000 0.9646883
## [2,] 0.9646883 1.0000000# ベーカー行列
dend_list %>% cor.dendlist(method = "baker")## [,1] [,2]
## [1,] 1.0000000 0.9622885
## [2,] 0.9622885 1.00000002つのツリー間の相関
# コーフェン行列
cor_cophenetic(dend1, dend2)## [1] 0.9646883# ベーカー行列
cor_bakers_gamma(dend1, dend2)## [1] 0.9622885複数の樹形図の同時比較
# 複数のデンドログラムを作成
dend1 <- df %>% dist %>% hclust("complete") %>% as.dendrogram
dend2 <- df %>% dist %>% hclust("single") %>% as.dendrogram
dend3 <- df %>% dist %>% hclust("average") %>% as.dendrogram
dend4 <- df %>% dist %>% hclust("centroid") %>% as.dendrogram
# リスト作成
dend_list <- dendlist("Complete" = dend1,
"Single" = dend2,
"Average" = dend3,
"Centroid" = dend4)
# 各デンドログラムの相関算出
cors <- dend_list %>% cor.dendlist()
cors %>% round(2)## Complete Single Average Centroid
## Complete 1.00 0.76 0.99 0.75
## Single 0.76 1.00 0.80 0.84
## Average 0.99 0.80 1.00 0.74
## Centroid 0.75 0.84 0.74 1.00# 可視化
cors %>% corrplot("pie", "lower")
樹形図の可視化
準備
＜ポイント＞
・USArrestsデータセットを使用
・階層クラスタリングを実行(ユーグリッド距離 + ウォード法)# データロード
data("USArrests")
# データ確認
USArrests %>% as_tibble() %>% head()## # A tibble: 6 x 4
## Murder Assault UrbanPop Rape
## <dbl> <int> <int> <dbl>
## 1 13.2 236 58 21.2
## 2 10.0 263 48 44.5
## 3 8.10 294 80 31.0
## 4 8.80 190 50 19.5
## 5 9.00 276 91 40.6
## 6 7.90 204 78 38.7# クラスタリング
# --- 距離測定 ：ユーグリッド距離
# --- リンケージ：ウォード法
hc <-
USArrests %>%
scale() %>%
dist(method = "euclidean") %>%
hclust(method = "ward.D2")樹形図(基礎)
＜ポイント＞
・fviz_dend()を使うことでggplot2ベースのデンドログラムが作成できる
・k引数でk個のクラスタにカットして色分けできる
・ggtheme引数でggplot2の各種テーマを指定可能# 基本的なデンドログラム
hc %>%
fviz_dend(cex = 0.5,
main = "Dendrogram - ward.D2",
xlab = "Objects",
ylab = "Distance",
sub = "")
# 指定カット数のクラスタに色分け
# --- # Cut in four groups
hc %>%
fviz_dend(k = 4,
cex = 0.5,
k_colors = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"),
color_labels_by_k = TRUE,
rect = TRUE,
rect_border = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"),
rect_fill = TRUE)## Warning in if (color == "cluster") color <- "default": 条件が長さが 2 以上
## なので、最初の 1 つだけが使われます
# テーマ設定も可能
hc %>%
fviz_dend(k = 4,
cex = 0.5,
k_colors = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"),
color_labels_by_k = TRUE,
ggtheme = theme_gray()
)
樹形図(応用)
＜ポイント＞
・水平型の方が場合によって見やすいことがある
・循環型は系列の重複を避けるために利用される
・樹木型プロットは系列の重複を避けるために利用される# 水平型デンドログラム
hc %>%
fviz_dend(k = 4,
cex = 0.4,
horiz = TRUE,
k_colors = "jco",
rect = TRUE,
rect_border = "jco",
rect_fill = TRUE)
# 循環型デンドログラム
hc %>%
fviz_dend(cex = 0.5,
k = 4,
k_colors = "jco",
type = "circular")
# 樹木型プロット
# --- {igraph}が必要
hc %>%
fviz_dend(k = 4,
k_colors = "jco",
type = "phylogenic",
repel = TRUE)
樹形図の拡大
＜ポイント＞
・デンドログラムはxlimとylimを指定して拡大することができる
・ylimは距離を示し、xlimは系列ベクトルを示すhc %>%
fviz_dend(xlim = c(1, 20),
ylim = c(1, 8))
サブツリー
＜ポイント＞
・サブツリーを利用して大規模デンドログラムの視認性を改善することができる
＜ステップ＞
・1：デンドログラムを作成&格納
・2：デンドログラムをカット
・3：サブツリーの表示# ステップ1：デンドログラムを作成&格納
dend_plot <- hc %>% fviz_dend(k = 4, cex = 0.5, k_colors = "jco")
dend_plot
# ステップ2：デンドログラムをカット
dend_cuts <- dend_plot %>% attr("dendrogram") %>% cut(h = 10)
dend_cuts## $upper
## 'dendrogram' with 2 branches and 2 members total, at height 13.51624
##
## $lower
## $lower[[1]]
## 'dendrogram' with 2 branches and 19 members total, at height 6.461866
##
## $lower[[2]]
## 'dendrogram' with 2 branches and 31 members total, at height 7.188189# ステップ3：サブツリーの表示
dend_cuts$upper %>% fviz_dend()## Warning in min(-diff(our_dend_heights)): min の引数に有限な値がありません:
## Inf を返します
dend_cuts$lower[[1]] %>% fviz_dend()
dend_cuts$lower[[2]] %>% fviz_dend()
# 参考：プロット形式を変更することも可能
dend_cuts$lower[[2]] %>% fviz_dend(type = "circular")
樹形図の可視化２
準備
＜ポイント＞
・USArrestsデータセットを使用
・階層クラスタリングを実行(ユーグリッド距離 + ウォード法)# データロード
data("USArrests")
# データ確認
USArrests %>% as_tibble() %>% head()## # A tibble: 6 x 4
## Murder Assault UrbanPop Rape
## <dbl> <int> <int> <dbl>
## 1 13.2 236 58 21.2
## 2 10.0 263 48 44.5
## 3 8.10 294 80 31.0
## 4 8.80 190 50 19.5
## 5 9.00 276 91 40.6
## 6 7.90 204 78 38.7# クラスタリング
# --- 距離測定 ：ユーグリッド距離
# --- リンケージ：ウォード法
hc <-
USArrests %>%
scale() %>%
dist(method = "euclidean") %>%
hclust(method = "ward.D2")plot.hclust()
＜ポイント＞
・hclust()によって生成されたオブジェクトからプロットを作成
＜構文＞
plot(x,
labels = NULL,
hang = 0.1,
main = "Cluster dendrogram",
sub = NULL,
xlab = NULL,
ylab = "Height", ...)
＜引数＞
x ：hclust()によって生成されたオブジェクト
labels：ツリーの葉のラベルの文字ベクトル(デフォルト値は行名)
hang ：ラベルがプロットの残りの部分の下に掛かるプロットの高さの割合
負の値を指定すると、ラベルは0からハングダウン
main sub xlab ylab：タイトルの文字列# とりあえずプロット
hc %>% plot()
# ラベル揃え + 文字縮小
hc %>% plot(hang = -1, cex = 0.6)
plot.dendrogram()
＜ポイント＞
・hclust()の結果を直接でなく、as.dendrogram()を間にかませる
・plot.hclust()よりもプロットの表現が柔軟になる
＜構文＞
plot(x,
type = c("rectangle", "triangle"),
horiz = FALSE
)
＜構文＞
x ：クラス樹状図のオブジェクト
tyep ：プロットのタイプ
horiz ：樹状図が水平にするか# 樹形図を四辺形で表示
hc %>%
as.dendrogram() %>%
plot(type = "rectangle", ylab = "Height")
# 樹形図を三角形で表示
hc %>%
as.dendrogram() %>%
plot(type = "triangle", ylab = "Height")
# ズーム
hc %>%
as.dendrogram() %>%
plot(xlim = c(1, 20), ylim = c(1,8))
ヒートマップ
概要
＜ポイント＞
・ヒートマップは階層的クラスタリングを可視化するための別の方法
--- サンプルとフィーチャーのクラスターを同時に視覚化することができる
＜作成方法＞
・stats::heatmap()
・gplotsR::heatmap2()
・pheatmap::pheatmap()
・d3heatmap::d3heatmap()データ準備
＜ポイント＞
・mtcarsデータセットのデータをそのまま使用する
--- hclust()はなし# データ確認
mtcars %>% as_tibble()## # A tibble: 32 x 11
## mpg cyl disp hp drat wt qsec vs am gear carb
## * <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 21.0 6.00 160 110 3.90 2.62 16.5 0 1.00 4.00 4.00
## 2 21.0 6.00 160 110 3.90 2.88 17.0 0 1.00 4.00 4.00
## 3 22.8 4.00 108 93.0 3.85 2.32 18.6 1.00 1.00 4.00 1.00
## 4 21.4 6.00 258 110 3.08 3.22 19.4 1.00 0 3.00 1.00
## 5 18.7 8.00 360 175 3.15 3.44 17.0 0 0 3.00 2.00
## 6 18.1 6.00 225 105 2.76 3.46 20.2 1.00 0 3.00 1.00
## 7 14.3 8.00 360 245 3.21 3.57 15.8 0 0 3.00 4.00
## 8 24.4 4.00 147 62.0 3.69 3.19 20.0 1.00 0 4.00 2.00
## 9 22.8 4.00 141 95.0 3.92 3.15 22.9 1.00 0 4.00 2.00
## 10 19.2 6.00 168 123 3.92 3.44 18.3 1.00 0 4.00 4.00
## # ... with 22 more rows# 基準化
df <- mtcars %>% scale()
df %>% summary()## mpg cyl disp hp
## Min. :-1.6079 Min. :-1.225 Min. :-1.2879 Min. :-1.3810
## 1st Qu.:-0.7741 1st Qu.:-1.225 1st Qu.:-0.8867 1st Qu.:-0.7320
## Median :-0.1478 Median :-0.105 Median :-0.2777 Median :-0.3455
## Mean : 0.0000 Mean : 0.000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.4495 3rd Qu.: 1.015 3rd Qu.: 0.7688 3rd Qu.: 0.4859
## Max. : 2.2913 Max. : 1.015 Max. : 1.9468 Max. : 2.7466
## drat wt qsec vs
## Min. :-1.5646 Min. :-1.7418 Min. :-1.87401 Min. :-0.868
## 1st Qu.:-0.9661 1st Qu.:-0.6500 1st Qu.:-0.53513 1st Qu.:-0.868
## Median : 0.1841 Median : 0.1101 Median :-0.07765 Median :-0.868
## Mean : 0.0000 Mean : 0.0000 Mean : 0.00000 Mean : 0.000
## 3rd Qu.: 0.6049 3rd Qu.: 0.4014 3rd Qu.: 0.58830 3rd Qu.: 1.116
## Max. : 2.4939 Max. : 2.2553 Max. : 2.82675 Max. : 1.116
## am gear carb
## Min. :-0.8141 Min. :-0.9318 Min. :-1.1222
## 1st Qu.:-0.8141 1st Qu.:-0.9318 1st Qu.:-0.5030
## Median :-0.8141 Median : 0.4236 Median :-0.5030
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 1.1899 3rd Qu.: 0.4236 3rd Qu.: 0.7352
## Max. : 1.1899 Max. : 1.7789 Max. : 3.2117heatmap()
＜構文＞
heatmap(x, scale = "row")
＜引数＞
x ：数値行列
scale ：値が行方向または列方向のいずれかにセンタリング&スケーリングされるか# デフォルトプロット
df %>% heatmap(scale = "none")
# カラー取得
col <- colorRampPalette(c("red", "white", "blue"))(256)
col %>% glimpse()## chr [1:256] "#FF0000" "#FF0202" "#FF0404" "#FF0606" ...# カラーパレット
col <- colorRampPalette(brewer.pal(10, "RdYlBu"))(256)
col %>% glimpse()## chr [1:256] "#A50026" "#A60126" "#A80326" "#AA0526" ...# ヒートマップにカラーパレットを適用
df %>%
heatmap(scale = "none",
col = col,
RowSideColors = rep(c("blue", "pink"), each = 16),
ColSideColors = c(rep("purple", 5), rep("orange", 6)))
heatmap.2()
＜構文＞
heatmap.2(x, scale = "row")
＜引数＞
x ：数値行列
scale ：値が行方向または列方向のいずれかにセンタリング&スケーリングされるかdf %>%
gplots::heatmap.2(scale = "none",
col = bluered(100),
trace = "none",
density.info = "none")
クラスタリングの評価と検証
傾向の評価
準備
＜ポイント＞
・irisデータセットの数値部分を使用
・ランダムデータで同様のサイズのデータセットを作成# データ確認
iris %>% head(3)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa# Species列以外を使用
df <- iris[, -5] %>% scale()
# 使用データと同じランダム数値のデータセットを作成
random_df <-
df %>%
apply(2, function(x){runif(length(x), min(x), (max(x)))}) %>%
as_tibble() %>%
scale()目視検査
＜ポイント＞
・PCAでPC1とPC2を抽出して散布図に表示
--- クラスタがある場合は散布図に傾向が出る# irisデータセット
df %>%
prcomp() %>%
fviz_pca_ind(title = "PCA - Iris data",
habillage = iris$Species,
palette = "jco",
geom = "point",
ggtheme = theme_classic(),
legend = "bottom")
# ランダムデータ
random_df %>%
prcomp() %>%
fviz_pca_ind(title = "PCA - Random data",
geom = "point",
ggtheme = theme_classic())
傾向評価の理由
＜ポイント＞
・クラスタが適切に分かれているかに関らずクラスター分析は可能
--- 事前の傾向評価をしていないと不適切なあてはめを吸いる可能性がある
# 乱数設定
set.seed(123)
# k-meansを実行
km.res1 <- df %>% kmeans(3)
km.res2 <- random_df %>% kmeans(3)
# k-means結果：iris
# --- クラスタの重複が比較的少ない
list(data = df,
cluster = km.res1$cluster) %>%
fviz_cluster(list(data = df,
cluster = km.res1$cluster),
ellipse.type = "norm",
geom = "point",
stand = FALSE,
palette = "jco",
ggtheme = theme_classic())
# k-means結果：ランダム
# --- クラスタの重複が多い
list(data = random_df,
cluster = km.res2$cluster) %>%
fviz_cluster(ellipse.type = "norm",
geom = "point",
stand = FALSE,
palette = "jco",
ggtheme = theme_classic())
# 参考：ランダムの階層クラスタリング
random_df %>%
dist() %>%
hclust() %>%
fviz_dend(k = 3,
k_colors = "jco",
as.ggplot = TRUE,
show_labels = FALSE)
統計的評価
# ホプキンス統計量
# --- iris
res <- df %>% get_clust_tendency(n = nrow(df)-1, graph = FALSE)
res$hopkins_stat## [1] 0.1815219# ホプキンス統計量
# --- ランダム
res <- random_df %>% get_clust_tendency(n = nrow(random_df)-1, graph = FALSE)
res$hopkins_stat## [1] 0.4878529# 参考：get_clust_tendency()の出力
res %>% attributes()## $names
## [1] "hopkins_stat" "plot"視覚的評価
＜ポイント＞
・クラスター傾向(VAT)アプローチによって視覚的に評価する
＜手順＞
・ユークリッド距離でオブジェクト間の距離行列を計算
・類似オブジェクトが互いに近くなるように距離行列を並び替え
--- 順序付き非類似行列(ODM)# iris
df %>%
dist() %>%
fviz_dist(show_labels = FALSE) +
labs(title = "Iris data")
# ランダム
random_df %>%
dist() %>%
fviz_dist(show_labels = FALSE) +
labs(title = "Random data")
最適クラスタ数の決定
はじめに
＜ポイント＞
・k-meansのような分割型クラスタリングを行う上で最適クラスタ数の決定が必要
・最適クラスタ数は主観的なアプローチとなる
＜使用関数＞
・factoextra::fviz_nbclust()
--- 最適クラスター数を決定するためのプロットを作成
--- fviz_nbclust
・NbClust::NbClust
--- クラスター数の決定メソッドを同時に計算
データ準備
# データロード
data("USArrests")
# 基準化
df <- USArrests %>% scale()
df %>% summary()## Murder Assault UrbanPop Rape
## Min. :-1.6044 Min. :-1.5090 Min. :-2.31714 Min. :-1.4874
## 1st Qu.:-0.8525 1st Qu.:-0.7411 1st Qu.:-0.76271 1st Qu.:-0.6574
## Median :-0.1235 Median :-0.1411 Median : 0.03178 Median :-0.1209
## Mean : 0.0000 Mean : 0.0000 Mean : 0.00000 Mean : 0.0000
## 3rd Qu.: 0.7949 3rd Qu.: 0.9388 3rd Qu.: 0.84354 3rd Qu.: 0.5277
## Max. : 2.2069 Max. : 1.9948 Max. : 1.75892 Max. : 2.6444# データ確認
df %>% head()## Murder Assault UrbanPop Rape
## Alabama 1.24256408 0.7828393 -0.5209066 -0.003416473
## Alaska 0.50786248 1.1068225 -1.2117642 2.484202941
## Arizona 0.07163341 1.4788032 0.9989801 1.042878388
## Arkansas 0.23234938 0.2308680 -1.0735927 -0.184916602
## California 0.27826823 1.2628144 1.7589234 2.067820292
## Colorado 0.02571456 0.3988593 0.8608085 1.864967207fviz_nbclust関数
＜構文＞
fviz_nbclust(x,
FUNcluster,
method = c("silhouette", "wss", "gap_stat"))
＜引数＞
・x ：数値行列又はデータフレーム
・FUNCluster：パーティション関数(kmeans / pam / clara)
・method ：最適クラスタ数の決定メソッド(silhouette / wss / gap_stat)エルボー分析
＜ポイント＞
・クラスターは４が適切であることを示唆# エルボー分析
df %>%
fviz_nbclust(kmeans, method = "wss") +
geom_vline(xintercept = 4, linetype = 2) +
labs(subtitle = "Elbow method")
シルエット分析
＜ポイント＞
・クラスターは２が適切であることを示唆# シルエット分析の可視化
df %>%
fviz_nbclust(kmeans, method = "silhouette") +
geom_vline(xintercept = 4, linetype = 2) +
labs(subtitle = "Silhouette method")
ギャップ統計量
＜ポイント＞
・エルボー分析とシルエット法の欠点は、クラスタリング特性のみを測定すること
・ギャップ統計量は統計的視点での最適数を測定することができる
・クラスターは４が適切であることを示唆# 乱数設定
set.seed(123)
# ギャップ統計量の可視化
df %>%
fviz_nbclust(kmeans, nstart = 25, method = "gap_stat", nboot = 50) +
labs(subtitle = "Gap statistic method")
NbClust()
＜ポイント＞
・最適なクラスタ数を選択するための30のインデックスを算出
#
nb <-
df %>%
NbClust::NbClust(distance = "euclidean",
min.nc = 2,
max.nc = 10,
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:
## * 10 proposed 2 as the best number of clusters
## * 2 proposed 3 as the best number of clusters
## * 8 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
## * 2 proposed 10 as the best number of clusters
##
## ***** Conclusion *****
##
## * According to the majority rule, the best number of clusters is 2
##
##
## *******************************************************************# 出力データ
nb %>% attributes()## $names
## [1] "All.index" "All.CriticalValues" "Best.nc"
## [4] "Best.partition"# インデックス
nb$All.index %>% colnames()## [1] "KL" "CH" "Hartigan" "CCC" "Scott"
## [6] "Marriot" "TrCovW" "TraceW" "Friedman" "Rubin"
## [11] "Cindex" "DB" "Silhouette" "Duda" "Pseudot2"
## [16] "Beale" "Ratkowsky" "Ball" "Ptbiserial" "Frey"
## [21] "McClain" "Dunn" "Hubert" "SDindex" "Dindex"
## [26] "SDbw"# 可視化
nb %>% fviz_nbclust()## Among all indices:
## ===================
## * 2 proposed 0 as the best number of clusters
## * 10 proposed 2 as the best number of clusters
## * 2 proposed 3 as the best number of clusters
## * 8 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
## * 2 proposed 10 as the best number of clusters
##
## Conclusion
## =========================
## * According to the majority rule, the best number of clusters is 2 .
クラスタ検証統計
ベストクラスター化アルゴリズムの選択
階層クラスタのP値の計算
データ準備
＜ポイント＞
・67の腫瘍を含む73の肺組織の916遺伝子の遺伝子発現プロファイルのデータ
・列はサンプルであり、行は遺伝子
・分析では30サンプル(30列)のデータセットを使用# データロード
data("lung")
# データ概要
lung[, 1:4] %>% as_tibble()## # A tibble: 916 x 4
## fetal_lung `232-97_SCC` `232-97_node` `68-96_Adeno`
## * <dbl> <dbl> <dbl> <dbl>
## 1 - 0.400 4.28 3.68 - 1.35
## 2 - 2.22 5.21 4.75 - 0.910
## 3 - 1.35 -0.840 -2.88 3.35
## 4 0.680 0.560 -0.450 - 0.200
## 5 NA 4.14 3.58 - 0.400
## 6 - 3.23 -2.84 -2.72 - 0.830
## 7 0.490 0.230 0.620 - 1.15
## 8 0.300 2.86 1.63 0.160
## 9 0.500 2.74 2.02 0.290
## 10 - 1.25 0.240 1.06 NA
## # ... with 906 more rows# 30行のデータセットを作成
df <- sample(1:73, 30) %>% lung[, .]
df %>% dim()## [1] 916 30pvclust()
＜構文＞
pvclust(data,
method.hclust = "average",
method.dist = "correlation",
nboot = 1000)
＜引数＞
data ：数値データマトリックスまたはデータフレーム。
method.hclust ：階層的クラスタリングで使用される凝集メソッド
method.dist ：使用される距離測定値
nboot ：ブートストラップ数(デフォルトは1000)# クラスター分析とp値測定
res.pv <-
df %>%
pvclust(method.dist = "cor",
method.hclust = "average",
nboot = 10)## Bootstrap (r = 0.5)... Done.
## Bootstrap (r = 0.6)... Done.
## Bootstrap (r = 0.7)... Done.
## Bootstrap (r = 0.8)... Done.
## Bootstrap (r = 0.9)... Done.
## Bootstrap (r = 1.0)... Done.
## Bootstrap (r = 1.1)... Done.
## Bootstrap (r = 1.2)... Done.
## Bootstrap (r = 1.3)... Done.
## Bootstrap (r = 1.4)... Done.# データ属性
res.pv %>% attributes()## $names
## [1] "hclust" "edges" "count" "msfit" "nboot" "r" "store"
##
## $class
## [1] "pvclust"# データ抽出
clusters <- res.pv %>% pvpick()
clusters## $clusters
## $clusters[[1]]
## [1] "306-99_node" "306-99_Adeno"
##
## $clusters[[2]]
## [1] "222-97_normal" "219-97_normal"
##
## $clusters[[3]]
## [1] "184-96_Adeno" "184-96_node"
##
## $clusters[[4]]
## [1] "313-99MT_Adeno" "313-99PT_Adeno"
##
## $clusters[[5]]
## [1] "314-99_SCLC" "315-99_node"
##
## $clusters[[6]]
## [1] "59-96_SCC" "139-97_LCLC"
##
## $clusters[[7]]
## [1] "68-96_Adeno" "137-96_Adeno"
##
## $clusters[[8]]
## [1] "220-97_SCC" "219-97_SCC" "3-SCC" "232-97_node"
## [5] "239-97_SCC" "246-97_SCC_c" "245-97_node" "157-96_SCC"
##
##
## $edges
## [1] 1 2 3 4 6 8 10 17# プロット作成
res.pv %>% plot(hang = -1, cex = 0.5)
res.pv %>% pvrect()
ParPvclust()
＜ポイント＞
・並列計算バージョンの関数を使うとブートストラップにかかる時間が大幅に短縮できる
＜構文＞
parPvclust(cl = NULL,
data, method.hclust = "average",
method.dist = "correlation",
nboot = 1000,
iseed = NULL)
＜引数＞
cl ：並列クラスター
data ：数値データマトリックスまたはデータフレーム。
data, method.hclust ：階層的クラスタリングで使用される凝集メソッド
method.dist ：使用される距離測定値
nboot ：ブートストラップ数(デフォルトは1000)# 並列準備
cl <- makeCluster(2, type = "PSOCK")
# parallel version of pvclust
res.pv <-
df %>%
parPvclust(cl = cl,
data = .,
nboot = 1000)## Warning in parPvclust(cl = cl, data = ., nboot = 1000): "parPvclust" has been integrated into pvclust (with "parallel" option).
## It is available for back compatibility but will be unavailable in the future.## Multiscale bootstrap... Done.# 並列終了
cl %>% stopCluster()
# データ属性
res.pv %>% attributes()## $names
## [1] "hclust" "edges" "count" "msfit" "nboot" "r" "store"
##
## $class
## [1] "pvclust"