第11講 クラスタ分析

非階層的方法と分析の評価

Published

December 5, 2025

準備

以下で利用する共通パッケージを読み込む．

library(conflicted) 
conflicts_prefer(   
    dplyr::filter(),
    dplyr::select(),
    dplyr::lag(),
    )
library(tidyverse)  
library(GGally)
library(broom)      
library(ggfortify)
library(ggrepel)
library(cluster)
library(ggdendro)
#' macOSのための日本語表示の設定
if(Sys.info()["sysname"] == "Darwin") { # macOSか調べる
    jp_font <- "HiraMaruProN-W4"
    theme_update(text = element_text(family = jp_font))
    update_geom_defaults("text", list(family = theme_get()$text$family))
    update_geom_defaults("label", list(family = theme_get()$text$family))
    update_geom_defaults("text_repel", list(family = theme_get()$text$family))
    update_geom_defaults("label_repel", list(family = theme_get()$text$family))
    par(family = jp_font)
} else {jp_font <- NULL}

データセット

以下では2つのデータセットを使用する．

japan_social2.csv の概要

総務省統計局より取得した都道府県別の社会生活統計指標の一部

japan_social.csv を日本語化したもの

  • 都道府県名 :
  • 地方区分 :
  • 森林面積割合: 森林面積割合 (%) 2014年
  • 農業産出額 : 就業者１人当たり農業産出額(販売農家）(万円) 2014年
  • 人口割合 : 全国総人口に占める人口割合 (%) 2015年
  • 土地生産性 : 土地生産性（耕地面積１ヘクタール当たり）(万円) 2014年
  • 商品販売額 : 商業年間商品販売額［卸売業＋小売業］（事業所当たり）(百万円) 2013年
  • 参考 : https://www.e-stat.go.jp/SG1/estat/List.do?bid=000001083999&cycode=0

データの読み込み方の例．

js_data <-
    read_csv("data/japan_social2.csv",
             col_types = list(地方区分 = col_factor())) # 地方区分を因子化
omusubi.csv の概要

「ごはんを食べよう国民運動推進協議会」(平成30年解散) による 好きなおむすびの具に関するwebアンケートを都道府県別にまとめたもの - (閉鎖) http://www.gohan.gr.jp/result/09/anketo09.html

アンケートの概要は以下のとおりで，ファイルは質問2の結果を整理したもの

  • 【応募期間】 2009年1月4日～2009年2月28日
  • 【応募方法】 インターネット、携帯ウェブ
  • 【質問内容】
    1. おむすびを最近１週間に、何個食べましたか？ そのうち市販のおむすびは何個でしたか？
    2. おむすびの具では何が一番好きですか？
      • A.梅 B.鮭 C.昆布 D.かつお E.明太子 F.たらこ Ｇ.ツナ H.その他
    3. おむすびのことをあなたはなんと呼んでいますか？
      • A.おにぎり B.おむすび C.その他
    4. おむすびといえば、どういう形ですか？
      • A.三角形 B.丸形 C.俵形 D.その他
  • 【回答者数】
    • 男性 9,702人 32.0%
    • 女性 20,616人 68.0%
    • 総数 30,318人 100.0%

データの読み込み方の例．都道府県名や地方区分を日本語にする場合．

om_data <- 
    bind_cols( 
        read_csv(file = "data/omusubi.csv"),       # データの読み込み
        read_csv(file = "data/prefecture.csv")) |> # 日本語表記・地方の情報を追加
    select(jp,area_jp,ume:etc) |>
    set_names(c("都道府県名","地方区分",
                "梅","鮭","昆布","鰹","明太子","鱈子","ツナ","その他"))

k-means による非階層的クラスタリング

k-meansに関連した関数

関数 stats::kmeans() 
kmeans(x, centers, iter.max = 10, nstart = 1,
       algorithm = c("Hartigan-Wong", "Lloyd", "Forgy",
                     "MacQueen"), trace = FALSE)
#' x: データフレーム
#' centers: クラスタ数
#' iter.max: 最大繰り返し数
#' nstart: 初期値の候補数
#' algorithm: 最適化法の指定．既定値は 'Hartigan-Wong'
#' 詳細は '?stats::kmeans' を参照
関数 ggfortify::autoplot() 
autoplot(object, data = NULL, colour = "cluster", ...)
#' object: stats::kmeans(), cluster::pam() などの返値
#' data: クラスタリングに用いたデータ (kmeansの場合に必要)
#' 詳細は '?ggfortify::autoplot.kmeans' を参照

問題

  • データセット japan_social2.csv を用いて以下を確認しなさい．
    • 関数 stats::kmeans() を用いて各変数平均0，分散1に標準化 (関数 base::scale() を利用)したデータを7クラスタに分割しなさい．
    • 各クラスタ内の県名を表示しなさい．
    • 2次元散布図に各クラスタを表示しなさい．
  • データセット omusubi.csv でも確認しなさい．

解答欄

解答例

データの読み込みを行う． データ点のラベルを行名から取得して保存する仕様の base R の関数のために data.frame クラスも用意しておく．

js_data <-                                            # tibble クラス
    read_csv("data/japan_social2.csv",
             col_types = list(地方区分 = col_factor())) # 地方区分を因子化
js_df <-                                              # data.frame クラス
    js_data |>
    column_to_rownames(var = "都道府県名") |>           # 都道府県名を行名に変換
    select(-地方区分)                                   # 距離計算に不要な地方区分は除く

関数 stats::kmeans() を用いて k-平均法を実行する． 結果は確率的なので，必要に応じて乱数のシード値を指定する

set.seed(1234)          # 乱数のシード値を指定
k <- 7                  # 分割数を指定
js_kmeans <-
    js_df |>
    scale() |>          # 標準化(平均0 分散1)
    kmeans(centers = k, # クラスタ数
           nstart = 20) # 初期値を20回変更して試す

関数 stats::kmeans() の返値 cluster の情報を利用して．各クラスター内の県名を表示する． kmeans.object$cluster または kmeans.object[["cluster"]] を用いればよい．

for(i in 1:k) {
    cat("=== cluster",i,"===\n")
    which(js_kmeans[["cluster"]] == i) |> names() |> print()
}
=== cluster 1 ===
 [1] "岩手県" "宮城県" "秋田県" "山形県" "福島県" "新潟県" "富山県" "石川県"
 [9] "福井県" "三重県" "滋賀県" "兵庫県" "島根県" "岡山県" "広島県" "山口県"
=== cluster 2 ===
[1] "東京都"
=== cluster 3 ===
[1] "宮崎県"   "鹿児島県"
=== cluster 4 ===
[1] "北海道"
=== cluster 5 ===
 [1] "青森県" "茨城県" "栃木県" "群馬県" "静岡県" "香川県" "佐賀県" "長崎県"
 [9] "熊本県" "沖縄県"
=== cluster 6 ===
 [1] "山梨県"   "長野県"   "岐阜県"   "京都府"   "奈良県"   "和歌山県"
 [7] "鳥取県"   "徳島県"   "愛媛県"   "高知県"   "大分県"  
=== cluster 7 ===
[1] "埼玉県"   "千葉県"   "神奈川県" "愛知県"   "大阪府"   "福岡県"

関数 ggfortify::autoplot() を用いて2次元でのクラスタ表示を行う． 前回の主成分分析を利用したクラスタの表示とほぼ同様な記述で行うことができる．

js_kmeans |>
    autoplot(data = scale(js_df),       # PCAを計算するデータの指定が必須
             frame = TRUE,              # クラスタ毎に枠を付ける
             frame.type = "convex",     # 凸包 "convex"・楕円 "norm,t" が指定できる
             label = TRUE,              # ラベルを付加
             label.repel = TRUE,        # 重なりを回避(ラベルが消える場合もあるので注意)
             label.size = 3,            # ラベルの大きさ
             label.show.legend = FALSE) # 凡例の中のアルファベットを除く

関数 cluster::clusplot() を用いる場合は以下のようにすればよい．

clusplot(x = js_df, 
         clus = js_kmeans$cluster,              # クラスタ番号
         stand = TRUE,                          # データの標準化を行う
         lines = 0, labels = 3,                 # 表示の指定
         main = NULL, sub = NULL, cex = 0.8,    # タイトルなどの調整
         col.p = rainbow(k)[js_kmeans$cluster], # 虹色で色付け
         col.clus = "orange", shade = FALSE)    # クラスタ囲みの指定

データセット omusubi.csv でも試してみる． ここでは tibble クラスのみでの実装の例を示す．

om_data <-
    bind_cols( # 日本語表記・地方の情報を追加
        read_csv(file = "data/omusubi.csv"),
        read_csv(file = "data/prefecture.csv")) |>
    select(jp,area_jp,ume:etc) |>
    set_names(c("都道府県名","地方区分",
                "梅","鮭","昆布","鰹","明太子","鱈子","ツナ","その他"))

k-平均法を実行する．

k <- 6 # 6分割で分析
om_kmeans <-
    om_data |>
    select(where(is.double)) |> # 計算対象の列を抽出
    sqrt() |>                   # Hellinger距離
    kmeans(centers = k,         # クラスタ数
           nstart = 20)         # 初期値を20回変更して試す

各クラスター内の県名を表示する．

for(i in 1:k) {
    cat("=== cluster",i,"===\n")
    om_data |>
        mutate(cluster = om_kmeans[["cluster"]]) |>
        filter(cluster == i) |>
        pull(都道府県名) |>
        print() 
}
=== cluster 1 ===
[1] "北海道" "青森県" "秋田県"
=== cluster 2 ===
 [1] "石川県"   "福井県"   "京都府"   "大阪府"   "兵庫県"   "奈良県"  
 [7] "和歌山県" "島根県"   "広島県"   "徳島県"   "愛媛県"  
=== cluster 3 ===
[1] "福島県" "群馬県" "山梨県" "長野県" "鳥取県" "佐賀県"
=== cluster 4 ===
[1] "岩手県" "山形県" "新潟県" "沖縄県"
=== cluster 5 ===
 [1] "宮城県"   "茨城県"   "栃木県"   "埼玉県"   "千葉県"   "東京都"  
 [7] "神奈川県" "静岡県"   "愛知県"   "三重県"   "高知県"   "宮崎県"  
=== cluster 6 ===
 [1] "富山県"   "岐阜県"   "滋賀県"   "岡山県"   "山口県"   "香川県"  
 [7] "福岡県"   "長崎県"   "熊本県"   "大分県"   "鹿児島県"

2次元でのクラスタ表示を行う．

om_kmeans |>
    autoplot(data = om_data |>
                 mutate(across(where(is.double),\(x)c(sqrt(x)))), 
             frame = TRUE, 
             frame.type = "norm", # 楕円で囲む
             label = TRUE, 
             label.label = "都道府県名", 
             label.repel = TRUE, 
             label.size = 3, 
             label.show.legend = FALSE)

関数 cluster::clusplot() を用いて表示する．

clusplot(x = om_data |>
             column_to_rownames(var = "都道府県名") |> 
             select(where(is.double)) |>
             sqrt(),
         clus = om_kmeans$cluster, # クラスタ番号
         stand = FALSE, # データは標準化しない
         lines = 0, labels = 3, # 表示の指定
         main = NULL, sub = NULL, cex = 0.8, # タイトルなどの調整
         col.p = rainbow(k)[om_kmeans$cluster], # 虹色で色付け
         col.clus = "orange", shade = FALSE)     # クラスタ囲みの指定

k-medoids による非階層的クラスタリング

k-medoids に関連した関数

関数 cluster::pam() 
pam(x, k, diss = inherits(x, "dist"),
    metric = c("euclidean", "manhattan"), 
    medoids = if(is.numeric(nstart)) "random",
    nstart = if(variant == "faster") 1 else NA,
    stand = FALSE, cluster.only = FALSE,
    do.swap = TRUE,
    keep.diss = !diss && !cluster.only && n < 100,
    keep.data = !diss && !cluster.only,
    variant = c("original", "o_1", "o_2", "f_3", "f_4", "f_5", "faster"),
    pamonce = FALSE, trace.lev = 0)
#' x: データフレーム，または距離行列 
#' k: クラスタの数
#' metric: 距離の指定(xがデータフレームの場合)
#' stand: 標準化(平均0，絶対偏差1)
#' 詳細は '?cluster::pam' を参照

問題

  • データセット japan_social2.csv を用いて以下を確認しなさい．

    • 関数 cluster::pam() を用いて変数平均0，絶対偏差1に標準化したデータを7クラスタに分割しなさい．
    • 各クラスタ内の県名を表示しなさい．
    • 2次元散布図に各クラスタを表示しなさい．

  • データセット omusubi.csv でも確認しなさい．

解答欄

解答例

関数 cluster::pam() を用いて k-medoids を実行する．

k <- 7
js_pam <-
    js_df |>  
    pam(stand = TRUE, # 正規化(平均0 絶対偏差1)
        k = k)        # クラスタ数の指定

各クラスター内の県名を表示する．

for(i in 1:k){
    cat("=== cluster",i,"===\n")
    which(js_pam[["clustering"]] == i) |> names() |> print()
}
=== cluster 1 ===
[1] "北海道"
=== cluster 2 ===
 [1] "青森県"   "茨城県"   "栃木県"   "群馬県"   "静岡県"   "佐賀県"  
 [7] "長崎県"   "熊本県"   "宮崎県"   "鹿児島県" "沖縄県"  
=== cluster 3 ===
 [1] "岩手県" "宮城県" "山形県" "新潟県" "石川県" "長野県" "岐阜県" "三重県"
 [9] "京都府" "兵庫県" "岡山県" "広島県" "香川県" "愛媛県"
=== cluster 4 ===
[1] "秋田県" "福島県" "富山県" "福井県" "滋賀県" "奈良県" "鳥取県" "島根県"
[9] "山口県"
=== cluster 5 ===
[1] "埼玉県"   "千葉県"   "神奈川県" "愛知県"   "大阪府"   "福岡県"  
=== cluster 6 ===
[1] "東京都"
=== cluster 7 ===
[1] "山梨県"   "和歌山県" "徳島県"   "高知県"   "大分県"

2次元でのクラスタ表示を行う

js_pam |> 
    autoplot(frame = TRUE, 
             frame.type = "convex", 
             label = TRUE, 
             label.repel = TRUE, 
             label.size = 3, 
             label.show.legend = FALSE)

Manhattan距離だとどのようになるか試してみる． いくつかのクラスタでメンバが変わっていることが確認できる．

js_df |>  
    pam(stand = TRUE, 
        metric = "manhattan",
        k = k) |>
    autoplot(frame = TRUE, 
             frame.type = "convex", 
             label = TRUE, 
             label.repel = TRUE, 
             label.size = 3, 
             label.show.legend = FALSE)

関数 cluster::clusplot() を用いる場合も kmeans とほぼ同様に行うことができる．

clusplot(x = js_df,
         clus = js_pam$clustering,
         stand = TRUE,
         lines = 0, labels = 3, 
         main = NULL, sub = NULL, cex = 0.8,
         col.p = rainbow(k)[js_pam$clustering],
         col.clus = "orange", shade = FALSE)

データセット omusubi.csv でも試してみる． k-medoids を実行する．

k <- 6
om_pam <- 
    om_data |>
    column_to_rownames(var = "都道府県名") |> 
    select(where(is.double)) |>
    sqrt() |>
    pam(k = k)

各クラスター内の県名を表示する．

for(i in 1:k){
    cat("=== cluster",i,"===\n")
    om_data |>
        mutate(cluster = om_pam[["clustering"]]) |>
        filter(cluster == i) |>
        pull(都道府県名) |>
        print() 
}
=== cluster 1 ===
[1] "北海道" "青森県" "秋田県"
=== cluster 2 ===
[1] "岩手県" "山形県" "新潟県" "沖縄県"
=== cluster 3 ===
 [1] "宮城県"   "茨城県"   "栃木県"   "群馬県"   "埼玉県"   "千葉県"  
 [7] "東京都"   "神奈川県" "山梨県"   "長野県"   "宮崎県"  
=== cluster 4 ===
[1] "福島県" "岐阜県" "愛知県" "鳥取県" "岡山県" "香川県" "佐賀県"
=== cluster 5 ===
 [1] "富山県"   "静岡県"   "三重県"   "滋賀県"   "京都府"   "大阪府"  
 [7] "兵庫県"   "奈良県"   "和歌山県" "島根県"   "愛媛県"   "高知県"  
[13] "福岡県"   "長崎県"   "熊本県"   "大分県"   "鹿児島県"
=== cluster 6 ===
[1] "石川県" "福井県" "広島県" "山口県" "徳島県"

2次元でのクラスタ表示を行う．

om_pam |> 
    autoplot(frame = TRUE, 
             frame.type = "convex", 
             label = TRUE, 
             label.repel = TRUE, 
             label.size = 3, 
             label.show.legend = FALSE)

関数 cluster::clusplot() で表示する場合は， クラスタリングの結果が om_pam$clustering に保管されていることに注意する．

clusplot(x = om_data |>
             column_to_rownames(var = "都道府県名") |> 
             select(where(is.double)) |>
             sqrt(),
         clus = om_pam$clustering,
         stand = TRUE,
         lines = 0, labels = 3, 
         main = NULL, sub = NULL, cex = 0.8,
         col.p = rainbow(k)[om_pam$clustering],
         col.clus = "orange", shade = FALSE)

凝集係数による距離の検討

凝集係数の計算

関数 cluster::agnes() の返値の情報

以下を利用して banner plot を描くことができる．

object[["ac"]] # 凝集係数の取得 (object$ac でも良い)
object[["height"]] # デンドログラムの高さ(結合時のクラスタ距離)
object[["order.lab"]] # デンドログラムの並び順のラベル
#' object: cluster::agnes() の返値

関数 summary(object) はこれらの情報をまとめて表示する．

base R での視覚化では，以下の関数を用いることができる．

#' cluster::plot.agnes() 系統樹および凝集係数の表示
plot(x, ask = FALSE, which.plots = NULL, main = NULL,
     sub = paste("Agglomerative Coefficient = ",round(x$ac, digits = 2)),
     adj = 0, nmax.lab = 35, max.strlen = 5, xax.pretty = TRUE, ...)
#' x: cluster::agnes() の返値
#' which.plots: 1 (banner plot), 2 (dendrogram)
#' banner plot の表示には以下の cluster::bannerplot() が呼出される
bannerplot(x, w = rev(x$height), fromLeft = TRUE,
           main=NULL, sub=NULL, xlab = "Height",  adj = 0,
           col = c(2, 0), border = 0, axes = TRUE, frame.plot = axes,
           rev.xax = !fromLeft, xax.pretty = TRUE,
           labels = NULL, nmax.lab = 35, max.strlen = 5,
           yax.do = axes && length(x$order) <= nmax.lab,
           yaxRight = fromLeft, y.mar = 2.4 + max.strlen/2.5, ...)

問題

  • データセット japan_social.csv を用いて以下を検討しなさい．

    • 関数 agnes() を用いて階層的クラスタリングを行いなさい． - 標準化: 行う - データ距離: ユークリッド距離，およびマンハッタン距離 - クラスタ距離: 群平均法
    • 凝集係数を用いて2つのデータ距離の評価を行いなさい．
    • 凝集係数が低いいくつかのデータを削除して評価しなさい．

解答欄

解答例

データ間をユークリッド距離で測り階層的クラスタリングを行う．

js_agnes_euc <- agnes(js_df,
                      metric="euclidean", # データ距離
                      stand=TRUE,       # 標準化
                      method="average")   # クラスタ距離
js_agnes_euc |> as.dendrogram() |>  
    ggdendrogram(rotate = FALSE, 
                 theme_dendro = FALSE) + 
    labs(title = "Euclidean distance",
         x = "Prefecture", y = "Height") +
    theme(axis.text.y = element_text(size = 9))

マンハッタン距離で階層的クラスタリングを行う．

js_agnes_man <- agnes(js_df,
                      metric="manhattan",
                      stand=TRUE,
                      method="average")
js_agnes_man |> as.dendrogram() |>  
    ggdendrogram(rotate = FALSE, 
                 theme_dendro = FALSE) + 
    labs(title = "Manhattan distance",
         x = "Prefecture", y = "Height") +
    theme(axis.text.y = element_text(size = 9))

それぞれの凝集係数を確認する．

js_agnes_euc[["ac"]]
[1] 0.8817664
js_agnes_man[["ac"]]
[1] 0.8750056

ユークリッド距離の方がわずかに良いことがわかる．

データ毎の凝集係数を表示する．1

1 一部のデータの距離が大きいと凝集係数は大きくなりがちであるが，その理由を考えてみよう．

tibble(x = js_agnes_euc[["height"]],
       y = length(js_agnes_euc[["height"]]):1) |>
    ggplot() +
    #' 各枝のheightと最大値で矩形を描く
    geom_rect(aes(xmin = x, xmax = max(x),
                  ymin = y, ymax = y+1), 
              fill = "orange", alpha = 0.6) + # 塗り潰し色と透明度を指定
    #' y軸にラベルを表示する
    scale_y_continuous(breaks = length(js_agnes_euc[["order.lab"]]):1,
                       expand = expansion(add = 0.5),
                       labels = js_agnes_euc[["order.lab"]]) +
    labs(title = "Euclidean distance") +
    theme(axis.text.y = element_text(size = 9))

tibble(x = js_agnes_man[["height"]],
       y = length(js_agnes_man[["height"]]):1) |>
    ggplot() +
    geom_rect(aes(xmin = x, xmax = max(x),
                  ymin = y, ymax = y+1), 
              fill = "orange", alpha = 0.6) + 
    scale_y_continuous(breaks = length(js_agnes_man[["order.lab"]]):1,
                       expand = expansion(add = 0.5),
                       labels = js_agnes_man[["order.lab"]]) +
    labs(title = "Manhattan distance") +
    theme(axis.text.y = element_text(size = 9))

bannerplot の情報を整理するための関数を定義してもよい．

my_bannerplot <- function(object, # agnesの返値
                          fill = "orange", alpha = 0.6, ...) {
    p <- tibble(x = object[["height"]],
                y = length(object[["height"]]):1) |>
        ggplot() +
        geom_rect(aes(xmin = x, xmax = max(x),
                      ymin = y, ymax = y+1), 
                  fill = fill, alpha = alpha) + 
        scale_y_continuous(breaks = length(object[["order.lab"]]):1,
                           expand = expansion(add = 0.5),
                           labels = object[["order.lab"]])
    p
}
my_bannerplot(js_agnes_euc,
              fill = "red", alpha = 0.6)

2 返値の要素を抽出する場合にはパイプ(|>)が使えない場合があるので注意する．

北海道，東京，宮崎，鹿児島を除いて再計算する．2

agnes(js_df |> slice(-c(1,13,45,46)),
      metric = "euclidean",
      stand = TRUE,
      method = "average")[["ac"]]
[1] 0.8069585
agnes(js_df |> slice(-c(1,13,45,46)),
      metric = "manhattan",
      stand = TRUE,
      method = "average")[["ac"]]
[1] 0.7815715

いずれにせよユークリッド距離の方が凝集係数は大きいことがわかる．

個別の係数を確認してみる．

js_df |> slice(-c(1,13,45,46)) |>
    agnes(metric = "euclidean",
          stand = TRUE,
          method = "average") |>
    my_bannerplot() + # 注で定義した関数を利用
    labs(title = "Euclidean distance") +
    theme(axis.text.y = element_text(size = 9))

  js_df |> slice(-c(1,13,45,46)) |>
    agnes(metric = "manhattan",
          stand = TRUE,
          method = "average") |>
    my_bannerplot() +
    labs(title = "Manhattan distance") +
    theme(axis.text.y = element_text(size = 9))

base R の関数を利用する場合は以下のようにすればよい．

デンドログラムを表示する．

plot(js_agnes_euc, which.plots=2, 
     main="euclidean") 

plot(js_agnes_man, which.plots=2, 
     main="manhattan")

データ毎の凝集係数を表示する．

plot(js_agnes_euc, which.plots=1, # banner plotの表示
     nmax.lab=50,   # 表示するラベルの上限 (標準は40)
     max.strlen=5,  # 表示するラベルの文字数の上限
     main="euclidean")

plot(js_agnes_man, which.plots=1,
     nmax.lab=50,  
     max.strlen=5,
     main="manhattan")

シルエット係数によるクラスタ数の検討

シルエット係数の計算

関数 cluster::pam() の返値の情報 
object[["silinfo"]] # シルエット係数に関する様々な情報
object[["silinfo"]][["widths"]] # 各データのシルエット係数
object[["silinfo"]][["clus.avg.widths"]] # 各クラスタのシルエット係数の平均
object[["silinfo"]][["avg.width"]] # シルエット係数の平均
#' object: cluster::agnes() の返値

関数 summary(object) はこれらの情報をまとめて表示する．

補助的な関数 
#' シルエット係数の取得
silhouette(x, ...)
#' x: cluster::pam() などの返値
silhouette(x, dist, dmatrix, ...)
#' x: stats::cutree() などの返値
#' dist/dmatrix: 距離行列または解離度を表す行列など
視覚化のための関数 (tidyverse) 
#' ggfortify::autoplot.silhouette() シルエット係数の表示
autoplot(object,
         colour = "red", linetype = "dashed", size = 0.5, bar.width = 1, ...)
#' object: cluster::silhouette() の返値
#' colour/linetype/size: reference line の設定

base R での視覚化では，以下の関数を用いることができる．

#' cluster::plot.partition() 2次元クラスタおよびシルエット係数の表示
plot(x, ask = FALSE, which.plots = NULL,
     nmax.lab = 40, max.strlen = 5, data = x$data, dist = NULL,
     stand = FALSE, lines = 2,
     shade = FALSE, color = FALSE, labels = 0, plotchar = TRUE,
     span = TRUE, xlim = NULL, ylim = NULL, main = NULL, ...)
#' x: cluster::pam() などの返値
#' which.plots: 1 (cluster plot), 2 (silhouette plot)
#' silhouette plot の表示には以下の cluster::plot.silhouette() が呼出される
plot(x, nmax.lab = 40, max.strlen = 5,
     main = NULL, sub = NULL, xlab = expression("Silhouette width "* s[i]),
     col = "gray",  do.col.sort = length(col) > 1, border = 0,
     cex.names = par("cex.axis"), do.n.k = TRUE, do.clus.stat = TRUE, ...)
#' x: cluster::silhouette() の返値

問題

データセット omusubi.csv を用いて以下を検討しなさい

  • Hellinger距離を用いて距離行列を計算しなさい
Note

\boldsymbol{p},\boldsymbol{q} を確率ベクトルとして定義される確率分布の距離 \begin{equation} d_{\mathrm{hel}}(\boldsymbol{p},\boldsymbol{q}) = \frac{1}{\sqrt{2}}d_{\mathrm{euc}}(\sqrt{\boldsymbol{p}},\sqrt{\boldsymbol{q}}) \end{equation}

  • クラスタ数4-10のシルエット係数を比較しなさい
  • 適当と思われるクラスタ数による分析を行いなさい
  • Euclid距離を用いて同様な分析を行いなさい

解答欄

解答例

クラスタ数 4-10 で平均シルエット係数を確認する．

om_df_hel <- # Hellinger距離を計算しやすくデータフレームを用意
    om_data |>
    column_to_rownames(var = "都道府県名") |> 
    select(where(is.double)) |>
    sqrt() 
for(k in 4:10){
    cat(pam(om_df_hel, k = k)$silinfo$avg.width,
        " (k = ", k, ")\n", sep="")
}
0.1723566 (k = 4)
0.1516094 (k = 5)
0.175647 (k = 6)
0.1996024 (k = 7)
0.2037189 (k = 8)
0.2064451 (k = 9)
0.1948242 (k = 10)

pam(om_df_hel, k = k)[["silinfo"]][["avg.width"]] と書いても良い．

k=7,8,9 (上位3つ) のシルエット係数を視覚化する．

for(k in 7:9) {
    gg <- pam(om_df_hel, k = k) |>
        silhouette() |>
        autoplot() +
        labs(title = paste("k =", k))
    print(gg) # ggplotはfor文内では明示的にprintする必要がある
}

ラベルも表示したい場合は繁雑だが以下のように指定することができる．

foo <-
    om_df_hel |>
    pam(k = 8) |>
    silhouette()
foo |>
    autoplot() +
    scale_x_continuous(
        breaks = 1:47,
        labels = foo |>
            as.data.frame() |>
            rownames_to_column() |>
            as_tibble() |>
            group_by(cluster) |>
            arrange(sil_width, .by_group = TRUE) |>
            pull(rowname)) # シルエット係数の並び順の都道府県名

または並び順は異なるが必要な情報をデータフレームにして， 自身で描画することもできる．

foo <- 
    om_df_hel |>
    pam(k = 8)
foo |>
    silhouette() |>
    as.data.frame() |>
    rownames_to_column() |>
    as_tibble() |>
    mutate(rowname = fct_rev(as_factor(rowname)),
           cluster = as_factor(cluster)) |>
    ggplot(aes(x = rowname, y = sil_width, fill = cluster)) +
    geom_bar(stat = "identity") +
    geom_hline(yintercept = foo[["silinfo"]][["avg.width"]],
               colour = "orange", linetype = "dashed") +
    labs(x = "都道府県名", y = "シルエット係数") +
    coord_flip()

悪いシルエット係数が少ないという意味で k=8 が良さそう．

k <- 8
om_df_hel |>
    pam(k = k) |>
    autoplot(frame = TRUE, 
             frame.type = "convex", 
             label = TRUE, 
             label.repel = TRUE, 
             label.size = 3, 
             label.show.legend = FALSE)

同様な描画は以下でも可能．

om_pam <-
    om_df_hel |>
    pam(k = k)
plot(om_pam,
     which.plot = 1, # cluster::clusplot のオプションを参考
     stand = TRUE,
     lines = 0, labels = 3, 
     main = "", sub = NULL, cex = 0.8, # タイトルと文字の大きさの調整
     col.p = rainbow(k)[om_pam$clustering], # クラスタ番号ごとに色付け
     col.clus = "orange", shade = FALSE) # クラスタを楕円で表示

Euclid距離による分析を行う． k = 5-10 で検証する．

for(k in 5:10) {
    foo <-
        om_data |>
        column_to_rownames(var = "都道府県名") |> 
        select(where(is.double)) |>
        pam(k = k)
    p <- foo |>
        silhouette() |>
        autoplot() +
        labs(title = paste("k =", k, "(Silhouette coef. =",
                           round(foo$silinfo$avg.width, digits = 3), ")"))
    print(p) 
}

悪いシルエット係数が少ないという意味で k=7 が良さそう．

k <- 7
om_data |>
    column_to_rownames(var = "都道府県名") |> 
    select(where(is.double)) |>
    pam(k = k) |>
    autoplot(frame = TRUE, 
             frame.type = "norm", 
             label = TRUE, 
             label.repel = TRUE, 
             label.size = 3, 
             label.show.legend = FALSE)

階層的クラスタリングでもシルエット係数を計算することができる．

om_agnes <- agnes(om_df_hel)
om_agnes |> as.dendrogram() |>
    ggdendrogram()

silhouette(cutree(om_agnes, k = k), 
           daisy(om_df_hel)) |> # 距離行列が必要
    autoplot()

シルエット係数のグラフはクラスタ毎に降順(大きいものが上)に並べ替えられている． グラフに合わせた要素名を取り出すには例えば以下のようにすれば良い．

silhouette(cutree(om_agnes, k = k), daisy(om_df_hel)) |>
    as_tibble() |>
    mutate(prefecture=rownames(om_df_hel)) |>
    arrange(desc(cluster), desc(sil_width))
# A tibble: 47 × 4
   cluster neighbor sil_width prefecture
     <dbl>    <dbl>     <dbl> <chr>     
 1       7        5     0     大分県    
 2       6        5     0     高知県    
 3       5        3     0.440 島根県    
 4       5        6     0.399 徳島県    
 5       5        3     0.392 兵庫県    
 6       5        6     0.391 福井県    
 7       5        7     0.356 広島県    
 8       5        6     0.340 大阪府    
 9       5        3     0.306 愛媛県    
10       5        3     0.306 富山県    
# ℹ 37 more rows