library(conflicted) # 関数名の衝突を警告
conflicts_prefer( # 優先的に使う関数を指定
dplyr::filter(),
dplyr::select(),
dplyr::lag(),
)
library(tidyverse)
library(GGally)
library(broom) # 解析結果を tibble 形式に集約
library(ggfortify)
library(cluster)
library(ggdendro)
#' macOSのための日本語表示の設定
if(Sys.info()["sysname"] == "Darwin") { # macOSか調べる
theme_update(text = element_text(family = "HiraMaruProN-W4"))
label_family <- "HiraMaruProN-W4"
} else {label_family <- NULL}統計データ解析
第10講 サンプルコード
準備
以下で利用する共通パッケージを読み込む．
距離の計算
問題
都道府県別の社会生活統計指標を用いて以下を確認しなさい．
#' データの読み込み方の例
js_df <- read_csv("data/japan_social.csv") |>
column_to_rownames(var = "Pref") |> # 'Pref'を行名に変換
select(-Area) # 地方名は除く- 正規化せずにユークリッド距離とマンハッタン距離の計算を行いなさい．
- 正規化して上記と同様の計算を行いなさい．
daisyによる正規化を用いて，関東の都県同士の距離を表示しなさい．- 大阪と四国の間の距離を表示しなさい．
- ユークリッド距離とマンハッタン距離の散布図を描き比較しなさい．
解答例
データを読み込む．
js_data <- read_csv("data/japan_social.csv") |>
mutate(Area = as_factor(Area))
js_df <- js_data |>
column_to_rownames(var = "Pref") |> # 'Pref'を行名に変換
select(-Area) # 距離計算に不要な地方名は除く以下の処理では js_data のまま扱うこともできるが，距離・類似度を計算する関数は js_df のような形式(base::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] Hokkaido Aomori Iwate Miyagi Akita
Hokkaido 0.0000 718.5670 824.99141 855.2632 891.44207
Aomori 718.5670 0.0000 115.20838 241.1784 198.89514
Iwate 824.9914 115.2084 0.00000 192.9282 90.65793
Miyagi 855.2632 241.1784 192.92825 0.0000 217.03608
Akita 891.4421 198.8951 90.65793 217.0361 0.00000as.matrix(js_daisy_euc)[1:5,1:5] Hokkaido Aomori Iwate Miyagi Akita
Hokkaido 0.0000 718.5670 824.99141 855.2632 891.44207
Aomori 718.5670 0.0000 115.20838 241.1784 198.89514
Iwate 824.9914 115.2084 0.00000 192.9282 90.65793
Miyagi 855.2632 241.1784 192.92825 0.0000 217.03608
Akita 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] Gumma Saitama Chiba Tokyo Kanagawa Niigata
Gumma 0.000000 6.255665 4.392306 14.089741 6.478754 3.816135
Saitama 6.255665 0.000000 2.916545 10.169610 2.761346 5.422816
Chiba 4.392306 2.916545 0.000000 11.392256 3.983992 7.228734
Tokyo 14.089741 10.169610 11.392256 0.000000 7.942358 15.016108
Kanagawa 6.478754 2.761346 3.983992 7.942358 0.000000 7.240533
Niigata 3.816135 5.422816 7.228734 15.016108 7.240533 0.000000as.matrix(js_daisy_man)[10:15,10:15] Gumma Saitama Chiba Tokyo Kanagawa Niigata
Gumma 0.000000 8.461062 5.947642 22.32430 8.860566 5.406536
Saitama 8.461062 0.000000 4.154908 17.01594 3.561477 7.208855
Chiba 5.947642 4.154908 0.000000 19.04507 5.590610 9.617467
Tokyo 22.324301 17.015938 19.045071 0.00000 14.103536 23.473777
Kanagawa 8.860566 3.561477 5.590610 14.10354 0.000000 9.617582
Niigata 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] "Hokkaido" "Aomori" "Iwate" "Miyagi" ...
- attr(*, "Size")= int 47
- attr(*, "Metric")= chr "euclidean"attr(js_daisy_euc, "Labels") # 県名を確認 [1] "Hokkaido" "Aomori" "Iwate" "Miyagi" "Akita" "Yamagata"
[7] "Fukushima" "Ibaraki" "Tochigi" "Gumma" "Saitama" "Chiba"
[13] "Tokyo" "Kanagawa" "Niigata" "Toyama" "Ishikawa" "Fukui"
[19] "Yamanashi" "Nagano" "Gifu" "Shizuoka" "Aichi" "Mie"
[25] "Shiga" "Kyoto" "Osaka" "Hyogo" "Nara" "Wakayama"
[31] "Tottori" "Shimane" "Okayama" "Hiroshima" "Yamaguchi" "Tokushima"
[37] "Kagawa" "Ehime" "Kochi" "Fukuoka" "Saga" "Nagasaki"
[43] "Kumamoto" "Oita" "Miyazaki" "Kagoshima" "Okinawa" #' attributes(js_daisy_euc)$Labels でも良い
#' もとのデータから Area 列が "Kanto" となる県を抽出してもよい
#' which(js_data[["Area"]]=="Kanto")
as.matrix(js_daisy_euc)[8:14, 8:14] Ibaraki Tochigi Gumma Saitama Chiba Tokyo Kanagawa
Ibaraki 0.000000 2.090955 3.038905 2.672169 2.143663 13.02122 4.091957
Tochigi 2.090955 0.000000 2.333057 3.412450 3.634136 13.36895 4.777569
Gumma 3.038905 2.333057 0.000000 4.305088 3.568353 12.65315 4.708828
Saitama 2.672169 3.412450 4.305088 0.000000 2.438165 11.62034 2.294246
Chiba 2.143663 3.634136 3.568353 2.438165 0.000000 12.22220 2.862094
Tokyo 13.021215 13.368948 12.653149 11.620338 12.222198 0.00000 11.176496
Kanagawa 4.091957 4.777569 4.708828 2.294246 2.862094 11.17650 0.000000as.matrix(js_daisy_man)[8:14, 8:14] Ibaraki Tochigi Gumma Saitama Chiba Tokyo Kanagawa
Ibaraki 0.000000 3.632295 5.448720 4.223827 3.742307 21.18333 7.728868
Tochigi 3.632295 0.000000 4.493538 6.322881 7.374602 22.86824 8.764702
Gumma 5.448720 4.493538 0.000000 8.461062 5.947642 22.32430 8.860566
Saitama 4.223827 6.322881 8.461062 0.000000 4.154908 17.01594 3.561477
Chiba 3.742307 7.374602 5.947642 4.154908 0.000000 19.04507 5.590610
Tokyo 21.183328 22.868238 22.324301 17.015938 19.045071 0.00000 14.103536
Kanagawa 7.728868 8.764702 8.860566 3.561477 5.590610 14.10354 0.000000dist/dissimilarity オブジェクトは距離以外の様々な属性 (attributes) を持つ．
glimpse(obj): オブジェクトの構造(structure)を見る(関数strでも良い)attributes(obj): 属性を表示(変更)するattr(obj,属性名): 特定の属性を表示(変更)する
大阪と四国の間の距離を表示する．1
1 1行または1列の部分行列を切り出す場合には，ベクトルとして扱われるので注意が必要．
as.matrix(js_daisy_euc)[27, 36:39, drop = FALSE] # 行列として表示 Tokushima Kagawa Ehime Kochi
Osaka 7.114222 5.186444 6.256155 7.557471as.matrix(js_daisy_man)["Osaka", # 1行なので標準ではベクトルとして扱われる
c("Tokushima","Kagawa","Ehime","Kochi")]Tokushima Kagawa Ehime Kochi
13.948439 8.653867 11.480234 15.265914 ユークリッド距離とマンハッタン距離の散布図を描く．
p <- tibble( # 列名に特殊な文字(空白など)を含む場合は `` で囲む
`Euclid dist.` = as.vector(js_daisy_euc),
`Manhattan dist.` = as.vector(js_daisy_man)) |>
ggplot(aes(x = `Euclid dist.`, y = `Manhattan dist.`)) +
geom_point(colour = "blue")
p + geom_abline(slope = 1, intercept = 0, colour = "red") # 基準線を付けて表示
p + xlim(0,3) + ylim(0,3) # 原点近傍のみ表示
2点間の距離の大小関係を考えると，ユークリッド距離とマンハッタン距離ではいくつか順序が入れ替わっていることが確認できる．
階層的クラスタリング
問題
都道府県別の社会生活統計指標を用いて以下の分析を行いなさい．
- 平均0，分散1に正規化したデータのユークリッド距離を用いて，群平均法による階層的クラスタリングを行いなさい
- クラスタ数を5つとして分割を行いなさい
解答例
正規化してユークリッド距離を測る
js_dist <- js_df |> scale() |> dist() # 'dist(scale(js_df))' でも良いクラスタリングを実行する．
js_hclust <- js_dist |>
hclust(method = "average") # 群平均法デンドログラムを描画する．
js_hclust |> as.dendrogram() |>
ggdendrogram(rotate = FALSE, # 横向きにしない
theme_dendro = FALSE) + # 'TRUE'とすれば無地
labs(title = "Euclidean + Average",
x = "Prefecture", y = "Height") +
theme(axis.text.y = element_text(size = 9))
少数のクラスタに分割してみる．
k <- 5 # 分割数を指定
js_clust <- cutree(js_hclust, k = k) # デンドログラムを分割
js_pref <- rownames(js_df) # 県名の取得
for(i in 1:k){
cat("=== cluster",i,"===\n")
print(js_pref[js_clust==i])
}=== cluster 1 ===
[1] "Hokkaido"
=== cluster 2 ===
[1] "Aomori" "Iwate" "Miyagi" "Akita" "Yamagata" "Fukushima"
[7] "Ibaraki" "Tochigi" "Gumma" "Niigata" "Toyama" "Ishikawa"
[13] "Fukui" "Yamanashi" "Nagano" "Gifu" "Shizuoka" "Mie"
[19] "Shiga" "Kyoto" "Hyogo" "Nara" "Wakayama" "Tottori"
[25] "Shimane" "Okayama" "Hiroshima" "Yamaguchi" "Tokushima" "Kagawa"
[31] "Ehime" "Kochi" "Saga" "Nagasaki" "Kumamoto" "Oita"
[37] "Okinawa"
=== cluster 3 ===
[1] "Saitama" "Chiba" "Kanagawa" "Aichi" "Osaka" "Fukuoka"
=== cluster 4 ===
[1] "Tokyo"
=== cluster 5 ===
[1] "Miyazaki" "Kagoshima"関数 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_df |>
prcomp() |> # 主成分分析
autoplot(data = tibble(cluster = factor(js_clust)),
colour = "cluster",
## 主成分分析の図をクラスタ毎に色分けするためのデータを追加
frame = TRUE, # クラスタ毎に枠を付ける
frame.type = "convex", # 凸包 "convex"・楕円 "norm,t" が指定できる
label = TRUE, # ラベルを付加
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 m <- which.max(table(js_clust)) # 最大クラスタの番号を取り出す
js_df_sub <- js_df[js_clust==m,]
js_hclust_sub <- js_df_sub |>
scale() |> dist() |> hclust(method = "average")
js_hclust_sub |> as.dendrogram() |>
ggdendrogram(rotate = FALSE, theme_dendro = FALSE) +
labs(title = "euclidean + average",
x = "prefecture", y = "distance") +
theme(axis.text.y = element_text(size = 9))
更に分割を進めてみる．
js_clust_sub <- cutree(js_hclust_sub, k = 4) # 最大クラスタをさらに分割
js_df_sub |>
prcomp() |> # 主成分分析
autoplot(data = tibble(cluster = factor(js_clust_sub)),
colour = 'cluster',
frame = TRUE,
frame.type = "convex",
label = TRUE,
label.repel = TRUE,
label.size = 3,
label.show.legend = FALSE)
階層的クラスタリング
問題
データ omusubi.csv を用いて以下の分析を行いなさい．
#' データの読み込み
om_data <- read_csv(file = "data/omusubi.csv")
om_df <- om_data |> column_to_rownames(var = "prefecture")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"))距離計算用のデータフレームを作成する．
om_df <- om_data |>
select(ume:etc,jp) |>
set_names(c("梅","鮭","昆布","鰹","明太子","鱈子","ツナ","その他","県名")) |>
column_to_rownames(var = "県名")データの散布図を描き，確認する．
om_data |>
select(ume:etc,jp,area_jp) |>
set_names(c("梅","鮭","昆布","鰹","明太子","鱈子","ツナ","その他","県名","地方")) |>
ggpairs(columns = 1:8,
upper = list(continuous = "cor"),
diag = list(continuous = "barDiag"),
lower = list(continuous = wrap("points", size = 0.5),
mapping = aes(colour = 地方)))
県別の人気比率を描画する．
om_data |>
select(ume:etc,jp) |>
set_names(c("梅","鮭","昆布","鰹","明太子","鱈子","ツナ","その他","県名")) |>
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_df |> 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 = "距離") +
theme(axis.text.y = element_text(size = 9))
クラスタ数7としてクラスタを表示する．
k <- 7
om_clust <- cutree(om_agnes, k = k) # クラスタを作成
om_df |>
prcomp() |>
autoplot(data = tibble(cluster = factor(om_clust)),
colour = 'cluster',
label = TRUE,
label.repel = TRUE,
label.size = 3,
label.family = label_family,
label.show.legend = FALSE)
クラスタを明示するための補助線を加える．
om_df |>
prcomp() |>
autoplot(data = tibble(cluster = factor(om_clust)),
colour = 'cluster',
frame = TRUE,
frame.type = "norm", # 各クラスタを正規分布で近似して楕円を描く
label = TRUE,
label.repel = TRUE,
label.size = 3,
label.family = label_family,
label.show.legend = FALSE)
ggplotでは繁雑となるが下記のような図を描くこともできる．
om_dendr <- dendro_data(om_agnes, type="rectangle") # ggplot用に変換
om_rect <- left_join(label(om_dendr),
tibble(label = rownames(om_df),
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 = "距離") +
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)
繁雑な処理を関数化するには例えば以下のようにすればよい．
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 = "距離") +
theme(axis.text.x = element_text(size = 9))
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_df,
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")