統計データ解析

第10講 サンプルコード

Published

November 26, 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("都道府県名","地方区分",
                "梅","鮭","昆布","鰹","明太子","鱈子","ツナ","その他"))

階層的クラスタ分析に用いる主な関数

クラスタ分析に関しては base R 系に強力なパッケージが用意されている．

  • stats : base R の基本的な統計に関するパッケージ
    • 関数 dist(), hclust(), kmeans() など
  • cluster : Kaufman and Rousseeuw (1990) にもとづくパッケージ
    • 関数 daisy(), agnes(), pam() など

一方 tidyverse ではこれらを利用して視覚化を行うための関数が提供されている．

  • ggfortify : ggplotによる描画を補助するパッケージ
    • 関数 autoplot() など
  • ggdendro : ggplotによるデンドログラム描画のパッケージ

クラスタ分析においては，データ間の距離の計算が基本となる．

関数 stats::dist() 
dist(x, method = "euclidean", diag = FALSE, upper = FALSE, p = 2)
#' x: データフレーム
#' method: 距離 (標準はユークリッド距離，他は"manhattan","minkowski"など)
#' diag: 対角成分を持たせるか 
#' upper: 上三角成分を持たせるか (標準は下三角成分のみ)
#' 返値は dist クラス
#' 詳細は '?stats::dist' を参照
関数 cluster::daisy() 
daisy(x, metric = c("euclidean", "manhattan", "gower"),
      stand = FALSE, type = list(), weights = rep.int(1, p),
      warnBin = warnType, warnAsym = warnType, warnConst = warnType,
      warnType = TRUE)
#' x: データフレーム
#' metric: 距離 (標準はユークリッド距離，他は"manhattan"など)
#' stand: 正規化(平均と絶対偏差の平均による)の有無
#' 返値は dissimilarity クラス
#' 詳細は '?cluster::daisy' を参照

階層的クラスタ分析を行うには関数 stats::hclust() または関数 cluster::agnes() を用いる．

階層的クラスタ分析 
hclust(d, method = "complete", members = NULL)
#' d: 距離行列
#' method: 分析法 (標準は最長距離法，他は"single","average"など)
#' 詳細は '?stats::hclust' を参照
agnes(x, diss = inherits(x, "dist"), metric = "euclidean",
      stand = FALSE, method = "average", par.method,
      keep.diss = n < 100, keep.data = !diss, trace.lev = 0)
#' x: データフレーム，または距離行列
#' metric: 距離 (標準はユークリッド距離，他は 'manhattan' など)
#' stand: 正規化(平均と絶対偏差の平均による)の有無
#' method: 分析法 (標準は群平均法，他は 'single', 'complete' など)
#' 詳細は '?cluster::agnes' を参照

分析や視覚化のためにいくつかの補助的な関数が用意されている．

補助的な関数

デンドログラムに基づいてクラスタの分割を行うことができる．

cutree(tree, k = NULL, h = NULL)
#' tree: stats::hclust() の返値
#' k: クラスタの数を指定して分割
#' h: クラスタの高さを指定して分割
#' 詳細は '?stats::cutree' を参照

デンドログラムの視覚化を行うには以下を用いる．

ggdendrogram(data,
             segments = TRUE, labels = TRUE, leaf_labels = TRUE,
             rotate = FALSE, theme_dendro = TRUE, ...)
#' data: stats::hclust(), stats::dendrogram() などの返値
#' 詳細は '?ggdend::ggdendrogram' を参照

base R のグラフィクスによるデンドログラムの表示は以下の関数を用いる．

plot(x, labels = NULL, hang = 0.1, check = TRUE,
     axes = TRUE, frame.plot = FALSE, ann = TRUE,
     main = "Cluster Dendrogram",
     sub = NULL, xlab = NULL, ylab = "Height", ...)
#' x: stats::hclust() の返値
#' 詳細は `stats::hclust` を参照

関数 cutree() とほぼ同様な仕様で， クラスタの分割表示を行うこともできる．

rect.hclust(tree, k = NULL, which = NULL, x = NULL, h = NULL,
            border = 2, cluster = NULL)
#' tree: stats::hclust() の返値
#' 詳細は 'stats::rect.hclust' を参照

2次元の散布図の作成には主成分分析などを利用し， クラスタの分割には関数 cutree() を利用すればよい．

2次元の視覚化 
autoplot(object, data = NULL, frame = FALSE,
         scale = 1, x = 1, y = 2,
         variance_percentage = TRUE, ...)
#' object: stats::prcomp() などの返値
#' data: 描画に必要な追加データ
#' frame: クラスタごと(colourなどで指定)に凸包または楕円を描画
#' 詳細は '?ggfortify::autoplot.pca_common()/ggbiplot()'を参照

base R のグラフィクスによる視覚化には以下を用いればよい．

clusplot(x, clus, diss = FALSE, stand = FALSE,
         lines = 2, shade = FALSE, color = FALSE,
         labels= 0, plotchar = TRUE,
         col.p = "dark green", col.txt = col.p, col.clus = 5，...)
#' x: データフレーム
#' clus: クラスタ分割
#' stand: 正規化の有無
#' lines: クラスタ間の繋がりの表示 (0:無，1:外，2:中心)
#' shade: 網掛けの有無
#' labels: ラベルの表示 (0:無，2:データとクラスタ, 3:データ, 4:クラスタ, など)
#' col.p/txt/clue: データ点・文字・クラスタの色指定
#' 詳細は '?cluster::clusplot.default' を参照
#' 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

距離の計算

問題

都道府県別の社会生活統計指標 (japan_social2.csv)を用いて以下を確認しなさい．

  • 正規化せずにユークリッド距離とマンハッタン距離の計算を行いなさい．
  • 正規化して上記と同様の計算を行いなさい．
  • 関数 cluster::daisy() による正規化を用いて，関東の都県同士の距離を表示しなさい．
  • 大阪と四国の間の距離を表示しなさい．
  • ユークリッド距離とマンハッタン距離の散布図を描き比較しなさい．

解答欄

解答例

データを読み込む．

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

以下の処理では js_data のまま扱うこともできるが， 距離・類似度を計算する関数は js_df のような形式 (関数 base::data.frame() で作成する data.frame クラス) を想定しているため用意しておく．

ユークリッド距離とマンハッタン距離を計算する．

js_dist_euc <- dist(js_df, method = "euclidean")
js_dist_man <- dist(js_df, method = "manhattan")
js_daisy_euc <- daisy(js_df, metric = "euclidean")
js_daisy_man <- daisy(js_df, metric = "manhattan")

両者が同じことを確認するために表示する． そのまま表示すると距離行列全体の下三角部分が表示されてしまうので， 行列に変換して一部(5行5列分)を表示する．

as.matrix(js_dist_euc)[1:5,1:5]
         北海道   青森県    岩手県   宮城県    秋田県
北海道   0.0000 718.5670 824.99141 855.2632 891.44207
青森県 718.5670   0.0000 115.20838 241.1784 198.89514
岩手県 824.9914 115.2084   0.00000 192.9282  90.65793
宮城県 855.2632 241.1784 192.92825   0.0000 217.03608
秋田県 891.4421 198.8951  90.65793 217.0361   0.00000
as.matrix(js_daisy_euc)[1:5,1:5]
         北海道   青森県    岩手県   宮城県    秋田県
北海道   0.0000 718.5670 824.99141 855.2632 891.44207
青森県 718.5670   0.0000 115.20838 241.1784 198.89514
岩手県 824.9914 115.2084   0.00000 192.9282  90.65793
宮城県 855.2632 241.1784 192.92825   0.0000 217.03608
秋田県 891.4421 198.8951  90.65793 217.0361   0.00000

正規化したユークリッド距離とマンハッタン距離を計算する． 関数 stats::dist() 単体では正規化する機能はないので，関数 base::scale() を併用する． 関数 cluster::daisy() には正規化のオプションがあるので，これを利用する．

js_dist_euc <- dist(scale(js_df), method = "euclidean")
js_dist_man <- dist(scale(js_df), method = "manhattan")
js_daisy_euc <- daisy(js_df, metric = "euclidean", stand = TRUE)
js_daisy_man <- daisy(js_df, metric = "manhattan", stand = TRUE)

関数 base::scale() と cluster::daisy() では正規化の方法が異なることに注意する．

as.matrix(js_dist_man)[10:15,10:15]
            群馬県    埼玉県    千葉県    東京都 神奈川県    新潟県
群馬県    0.000000  6.255665  4.392306 14.089741 6.478754  3.816135
埼玉県    6.255665  0.000000  2.916545 10.169610 2.761346  5.422816
千葉県    4.392306  2.916545  0.000000 11.392256 3.983992  7.228734
東京都   14.089741 10.169610 11.392256  0.000000 7.942358 15.016108
神奈川県  6.478754  2.761346  3.983992  7.942358 0.000000  7.240533
新潟県    3.816135  5.422816  7.228734 15.016108 7.240533  0.000000
as.matrix(js_daisy_man)[10:15,10:15]
            群馬県    埼玉県    千葉県   東京都  神奈川県    新潟県
群馬県    0.000000  8.461062  5.947642 22.32430  8.860566  5.406536
埼玉県    8.461062  0.000000  4.154908 17.01594  3.561477  7.208855
千葉県    5.947642  4.154908  0.000000 19.04507  5.590610  9.617467
東京都   22.324301 17.015938 19.045071  0.00000 14.103536 23.473777
神奈川県  8.860566  3.561477  5.590610 14.10354  0.000000  9.617582
新潟県    5.406536  7.208855  9.617467 23.47378  9.617582  0.000000

以下では関数 cluster::daisy() による正規化を用いる．

関東の都県同士の距離を表示する．

glimpse(js_daisy_euc)        # 距離行列のもつ情報を見る
 'dissimilarity' num [1:1081] 6.75 7.61 7.72 8.19 7.21 ...
 - attr(*, "Labels")= chr [1:47] "北海道" "青森県" "岩手県" "宮城県" ...
 - attr(*, "Size")= int 47
 - attr(*, "Metric")= chr "euclidean"
attr(js_daisy_euc, "Labels") # 県名を確認
 [1] "北海道"   "青森県"   "岩手県"   "宮城県"   "秋田県"   "山形県"  
 [7] "福島県"   "茨城県"   "栃木県"   "群馬県"   "埼玉県"   "千葉県"  
[13] "東京都"   "神奈川県" "新潟県"   "富山県"   "石川県"   "福井県"  
[19] "山梨県"   "長野県"   "岐阜県"   "静岡県"   "愛知県"   "三重県"  
[25] "滋賀県"   "京都府"   "大阪府"   "兵庫県"   "奈良県"   "和歌山県"
[31] "鳥取県"   "島根県"   "岡山県"   "広島県"   "山口県"   "徳島県"  
[37] "香川県"   "愛媛県"   "高知県"   "福岡県"   "佐賀県"   "長崎県"  
[43] "熊本県"   "大分県"   "宮崎県"   "鹿児島県" "沖縄県"  
#' attributes(js_daisy_euc)$Labels でも良い
#' もとのデータから地方区分列が '関東' となる県を抽出してもよい
#' which(js_data[["地方区分"]]=="関東")
as.matrix(js_daisy_euc)[8:14, 8:14]
            茨城県    栃木県    群馬県    埼玉県    千葉県   東京都  神奈川県
茨城県    0.000000  2.090955  3.038905  2.672169  2.143663 13.02122  4.091957
栃木県    2.090955  0.000000  2.333057  3.412450  3.634136 13.36895  4.777569
群馬県    3.038905  2.333057  0.000000  4.305088  3.568353 12.65315  4.708828
埼玉県    2.672169  3.412450  4.305088  0.000000  2.438165 11.62034  2.294246
千葉県    2.143663  3.634136  3.568353  2.438165  0.000000 12.22220  2.862094
東京都   13.021215 13.368948 12.653149 11.620338 12.222198  0.00000 11.176496
神奈川県  4.091957  4.777569  4.708828  2.294246  2.862094 11.17650  0.000000
as.matrix(js_daisy_man)[8:14, 8:14]
            茨城県    栃木県    群馬県    埼玉県    千葉県   東京都  神奈川県
茨城県    0.000000  3.632295  5.448720  4.223827  3.742307 21.18333  7.728868
栃木県    3.632295  0.000000  4.493538  6.322881  7.374602 22.86824  8.764702
群馬県    5.448720  4.493538  0.000000  8.461062  5.947642 22.32430  8.860566
埼玉県    4.223827  6.322881  8.461062  0.000000  4.154908 17.01594  3.561477
千葉県    3.742307  7.374602  5.947642  4.154908  0.000000 19.04507  5.590610
東京都   21.183328 22.868238 22.324301 17.015938 19.045071  0.00000 14.103536
神奈川県  7.728868  8.764702  8.860566  3.561477  5.590610 14.10354  0.000000

dist/dissimilarity オブジェクトは距離以外の様々な属性 (attributes) を持つ．

  • glimpse(obj) : オブジェクトの構造(structure)を見る(関数strでも良い)
  • attributes(obj) : 属性を表示(変更)する
  • attr(obj,属性名) : 特定の属性を表示(変更)する

大阪と四国の間の距離を表示する．1

1 1行または1列の部分行列を切り出す場合には，ベクトルとして扱われるので注意が必要．

as.matrix(js_daisy_euc)[27, 36:39, drop = FALSE] # 行列として表示
         徳島県   香川県   愛媛県   高知県
大阪府 7.114222 5.186444 6.256155 7.557471
as.matrix(js_daisy_man)["大阪府", # 1行なので標準ではベクトルとして扱われる
                        c("徳島県","香川県","愛媛県","高知県")]
   徳島県    香川県    愛媛県    高知県 
13.948439  8.653867 11.480234 15.265914

ユークリッド距離とマンハッタン距離の散布図を描く．

gg <-
    tibble(
        ユークリッド距離 = as.vector(js_daisy_euc),
        マンハッタン距離 = as.vector(js_daisy_man)) |>
    ggplot(aes(x = ユークリッド距離, y = マンハッタン距離)) +
    geom_point(colour = "blue")
gg + geom_abline(slope = 1, intercept = 0, colour = "red") # 基準線を付けて表示

gg + xlim(0,3) + ylim(0,3) # 原点近傍のみ表示

2点間の距離の大小関係を考えると，ユークリッド距離とマンハッタン距離ではいくつか順序が入れ替わっていることが確認できる．

階層的クラスタリング

問題

都道府県別の社会生活統計指標を用いて以下の分析を行いなさい．

  • 平均0，分散1に正規化したデータのユークリッド距離を用いて，群平均法による階層的クラスタリングを行いなさい
  • クラスタ数を5つとして分割を行いなさい

解答欄

解答例

正規化してユークリッド距離を測る

js_dist <-
    js_data |>
    column_to_rownames(var = "都道府県名") |> # 都道府県名を行名に変換
    select(-地方区分) |>                      # 距離計算に不要な地方区分は除く
    scale() |>                               # 標準化
    dist()

クラスタリングを実行する．

js_hclust <-
    js_dist |>
    hclust(method = "average") # 群平均法

デンドログラムを描画する．

js_hclust |> 
    as.dendrogram() |>  
    ggdendrogram(rotate = FALSE, # 横向きにしない
                 theme_dendro = FALSE) + # 'TRUE'とすれば無地
    labs(title = "ユークリッド距離 + 群平均法",
         x = "都道府県名", y = "Height") +
    theme(axis.text.y = element_text(size = 9))

少数のクラスタに分割してみる．

k <- 5 # 分割数を指定
js_clust <- cutree(js_hclust, k = k) # デンドログラムを分割
js_pref <- js_data |> pull(都道府県名) # 県名の取得
for(i in 1:k){
    cat("=== cluster",i,"===\n")
    print(js_pref[js_clust==i])
}
=== cluster 1 ===
[1] "北海道"
=== cluster 2 ===
 [1] "青森県"   "岩手県"   "宮城県"   "秋田県"   "山形県"   "福島県"  
 [7] "茨城県"   "栃木県"   "群馬県"   "新潟県"   "富山県"   "石川県"  
[13] "福井県"   "山梨県"   "長野県"   "岐阜県"   "静岡県"   "三重県"  
[19] "滋賀県"   "京都府"   "兵庫県"   "奈良県"   "和歌山県" "鳥取県"  
[25] "島根県"   "岡山県"   "広島県"   "山口県"   "徳島県"   "香川県"  
[31] "愛媛県"   "高知県"   "佐賀県"   "長崎県"   "熊本県"   "大分県"  
[37] "沖縄県"  
=== cluster 3 ===
[1] "埼玉県"   "千葉県"   "神奈川県" "愛知県"   "大阪府"   "福岡県"  
=== cluster 4 ===
[1] "東京都"
=== cluster 5 ===
[1] "宮崎県"   "鹿児島県"

関数 stats::hclust() では base R の関数 plot() を利用して下記のような方法で簡便にクラスタの分割を表示できる． ggplot で同様なグラフを描く方法は後述する．

plot(js_hclust,
     hang = -1, # ラベルを揃えて表示
     cex = 0.8, # 文字のサイズを調整
     sub = "", xlab = "", main = "")
rect.hclust(js_hclust, k = k, border = "orange")

主成分分析を併用してクラスタの表示を行う．

js_data |>
    select(where(is.double)) |>         # 数値列のみ対象
    prcomp() |>                         # 主成分分析
    autoplot(data = js_data |>          # 元データ                      
                 mutate(cluster = factor(js_clust)), # +クラスタ番号
             colour = "cluster",        # クラスタごとに色分け
             frame = TRUE,              # クラスタ毎に枠を付ける
             frame.type = "convex",     # 凸包 "convex"・楕円 "norm,t" が指定できる
             label = TRUE,              # ラベルを付加
             label.label = "都道府県名",  # 都道府県名をラベル
             label.repel = TRUE,        # 重なりを回避(ラベルが消える場合もあるので注意)
             label.size = 3,            # ラベルの大きさ
             label.show.legend = FALSE) # 凡例の中のアルファベットを除く

最大クラスタを取り出して，再評価する．

table(js_clust)                     # 最大を確認
js_clust
 1  2  3  4  5 
 1 37  6  1  2 
max_c <- which.max(table(js_clust)) # 最大クラスタの番号を取り出す
js_data_sub <- 
    js_data |>
    slice(which(js_clust == max_c))
js_hclust_sub <-
    js_data_sub |>
    column_to_rownames(var = "都道府県名") |> 
    select(-地方区分) |>                     
    scale() |>                             
    dist() |>
    hclust(method = "average") 
js_hclust_sub |>
    as.dendrogram() |>  
    ggdendrogram(rotate = FALSE, theme_dendro = FALSE) +
    labs(title = "ユークリッド距離 + 群平均法",
         x = "都道府県名", y = "Height") +
    theme(axis.text.y = element_text(size = 9))

更に分割を進めてみる．

js_clust_sub <- cutree(js_hclust_sub, k = 4) # 最大クラスタをさらに分割
js_data_sub |>
    select(where(is.double)) |>
    prcomp() |> 
    autoplot(data = js_data_sub |>
                 mutate(cluster = factor(js_clust_sub)),
             colour = "cluster",
             frame = TRUE, 
             frame.type = "convex", 
             label = TRUE,
             label.label = "都道府県名",
             label.repel = TRUE,
             label.size = 3,
             label.show.legend = FALSE)

階層的クラスタリング

問題

データ omusubi.csv を用いて以下の分析を行いなさい．

  • Hellinger距離を用いて距離行列を作成しなさい．

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

  • 群平均法による階層的クラスタリングを行いなさい

  • クラスタ数を定めて2次元でのクラスタ表示を作成しなさい

解答欄

解答例

データ omusubi.csv を読み込んで，整理する．

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("都道府県名","地方区分",
                "梅","鮭","昆布","鰹","明太子","鱈子","ツナ","その他"))

データの散布図を描き，確認する．

om_data |>
    ggpairs(columns = 3:10,
            upper = list(continuous = "cor"),
            diag = list(continuous = "barDiag"),
            lower = list(continuous = wrap("points", size = 0.5),
                         mapping = aes(colour = 地方区分)))

都道府県別の人気比率を描画する．

om_data |>
    select(-地方区分) |>
    pivot_longer(-都道府県名) |>
    mutate(都道府県名 = fct_rev(as_factor(都道府県名)), # 表示のため順番を反転する
           name = as_factor(name)) |> # ggplot(aes(x = 県名, y = value)) +
    ggplot(aes(y = 都道府県名, x = value)) +
    geom_bar(aes(fill = name),
             stat = "identity",
             position = position_stack(reverse=TRUE)) +
    labs(title = "おむすびの具 都道府県別人気アンケート (2009)",
         x = "人気比率", fill = "具材") +
    theme(axis.text.y = element_text(size = 9))

Hellinger 距離による階層的クラスタリングを行う．

om_agnes <-
    om_data |>
    column_to_rownames(var = "都道府県名") |> 
    select(where(is.double)) |>
    sqrt() |> 
    agnes()

1/2乗してEuclid距離を計算すればHellinger距離に比例した量が得られる. 定数倍まで厳密に計算するのであれば 1/sqrt(2)*daisy(sqrt(om_df/100)) とすればよい．

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

om_agnes |>
    as.dendrogram() |>  
    ggdendrogram(rotate = TRUE, theme_dendro = FALSE) +
    labs(title = "おむすびの具人気アンケート", x = "都道府県名", y = "Height") +
    theme(axis.text.y = element_text(size = 9))

クラスタ数7としてクラスタを表示する．

k <- 7
om_clust <- cutree(om_agnes, k = k)  # クラスタを作成
om_data |>
    select(where(is.double)) |>
    prcomp() |>
    autoplot(data = om_data |>
                 mutate(cluster = factor(om_clust)),
             colour = "cluster",
             label = TRUE,
             label.label = "都道府県名",
             label.repel = TRUE,
             label.size = 3,
             label.show.legend = FALSE)

クラスタを明示するための補助線を加える．

om_data |>
    select(where(is.double)) |>
    prcomp() |>
    autoplot(data = om_data |>
                 mutate(cluster = factor(om_clust)),
             colour = "cluster",
             frame = TRUE,
             frame.type = "norm", # 各クラスタを正規分布で近似して楕円を描く
             label = TRUE,
             label.label = "都道府県名",
             label.repel = TRUE,
             label.size = 3,
             label.show.legend = FALSE)

若干繁雑ではあるが， 下記のような図を描くこともできる．

#' クラスタごとの矩形の座標を構成するための自作関数
rect_agnes <- function(x, k) {
    rect <- left_join(label(dendro_data(x, type="rectangle")),
                      tibble(label = rownames(x$data),
                             cluster = cutree(x, k = k))) |>
        group_by(cluster) |>
        summarize(xmin = min(x)-0.3, xmax = max(x)+0.3) |>
        mutate(ymax = mean(sort(x$height, decreasing = TRUE)[k-0:1]))
}
om_agnes |>
    as.dendrogram() |>
    ggdendrogram(rotate = TRUE, theme_dendro = FALSE) +
    geom_rect(data = rect_agnes(om_agnes, k = k),
              aes(xmin = xmin, xmax = xmax, ymin = 0, ymax = ymax,
                  fill = as_factor(cluster)),
              alpha = 0.3, show.legend = FALSE) +
    labs(title = "クラスタの分割",
         x = "都道府県名", y = "Height") +
    theme(axis.text.x = element_text(size = 9))

関数化する前の愚直な実装は以下のようになる．

om_dendr <- dendro_data(om_agnes, type="rectangle") # ggplot用に変換
om_rect <- left_join(label(om_dendr),
                     tibble(label = rownames(om_agnes$data),
                            cluster = om_clust)) |>
    group_by(cluster) |>
    summarize(xmin = min(x)-0.3, xmax = max(x)+0.3)
om_ymax <- mean(sort(om_agnes$height, decreasing = TRUE)[k-0:1])
om_agnes |>
    as.dendrogram() |>
    ggdendrogram(rotate = TRUE, theme_dendro = FALSE) +
    labs(title = "クラスタの分割",
         x = "都道府県名", y = "Height") +
    theme(axis.text.x = element_text(size = 9)) +
    geom_rect(data = om_rect,
              aes(xmin = xmin, xmax = xmax, ymin = 0, ymax = om_ymax,
                  fill = as_factor(cluster)),
              alpha = 0.3, show.legend = FALSE)

package::cluster には base R 系の描画関数が含まれている． 例えば上記と同様なグラフは以下のようにして書ける．

if(Sys.info()["sysname"] == "Darwin") par(family = "HiraMaruProN-W4") # macOS用
plot(om_agnes, which.plot = 2, cex = 0.8, # 関数plot.agnes()を使う場合
     main = "Dendrogram of Omusubi Data")

plot(as.dendrogram(om_agnes), # 関数plot.dendrogram()を使う場合
     main = "Dendrogram of Omusubi Data")
rect.hclust(om_agnes, k = k, border = (1:k)+1)

clusplot(x = om_data |> 
             column_to_rownames(var = "都道府県名") |> 
             select(where(is.double)),
         clus = cutree(om_agnes, k = k),
         labels = 2,
         col.p = "green", col.txt = "blue", col.clus = "orange", cex = 0.8,
         main = "Cluster of Omusubi Data")