情報科学3
クラスタリング
矢内 勇生
2021-11-05（改訂：2021-11-09）
準備
必要なパッケージを読み込む。
pacman::p_load(tidyverse,
patchwork,
GGally,
clustrd)\(K\)平均法の仕組み
データの準備
df1 <- tibble(
x1 = c(1, 1, 2, 3, 6, 7, 7, 9, 3, 3),
x2 = c(4, 8, 7, 8, 4, 2, 3, 4, 2, 1)
)散布図を作る。
p <- ggplot(df1, aes(x = x1, y = x2)) +
geom_point()
plot(p)
\(K\)平均法
df1 に含まれる観測点を、\(k=3\)の\(K\)平均法で3つのクラスタに分けてみよう。
(1)
まず、ランダムに3点を選び、それらを別のクラスに分類する。
set.seed(2021-11-06)
df2 <- df1 %>%
mutate(cluster = sample(c(1:3, rep(NA, n() - 3)), size = n()),
cluster = factor(cluster))選んだ3点を表示する。
p1 <- p +
geom_point(data = filter(df2, !is.na(cluster)),
aes(x = x1, y = x2, color = cluster),
shape = 4, size = 4) +
theme(legend.position = "none")
plot(p1)
3点までの距離を測り、最も近いグループに分類する。今回は、ユークリッド距離を利用する。 平面上でユークリッド距離を測る関数を作る。
dist_euc <- function(x, y, a, b) {
# (x, y) と (a, b) のユークリッド距離
sqrt((x - a)^2 + (y - b)^2)
}クラスタ1〜3の座標を取り出す。
# クラスタ1
x1_c1 <- df2 %>% filter(cluster == 1) %>% pull(x1)
x2_c1 <- df2 %>% filter(cluster == 1) %>% pull(x2)
# クラスタ2
x1_c2 <- df2 %>% filter(cluster == 2) %>% pull(x1)
x2_c2 <- df2 %>% filter(cluster == 2) %>% pull(x2)
# クラスタ3
x1_c3 <- df2 %>% filter(cluster == 3) %>% pull(x1)
x2_c3 <- df2 %>% filter(cluster == 3) %>% pull(x2)各点とそれぞれのクラスタの「中心」との距離を測る。
df2 <- df2 %>%
mutate(dist_1 = dist_euc(x1, x2, x1_c1, x2_c1),
dist_2 = dist_euc(x1, x2, x1_c2, x2_c2),
dist_3 = dist_euc(x1, x2, x1_c3, x2_c3))各点を、最も距離が短いクラスに分類する。
df2 <- df2 %>%
mutate(cluster = case_when(
dist_1 <= dist_2 & dist_1 < dist_3 ~ 1,
dist_2 < dist_1 & dist_2 <= dist_3 ~ 2,
TRUE ~ 3))途中経過を可視化する。
p2 <- ggplot(df2, aes(x = x1, y = x2,
color = as.factor(cluster),
shape = as.factor(cluster))) +
geom_point() +
theme(legend.position = "none")
plot(p2)
(2)
ここまでの分類結果から、各クラスタの重心の座標を求める。
# クラスタ1
x1_c1 <- df2 %>%
filter(cluster == 1) %>%
pull(x1) %>%
mean()
x2_c1 <- df2 %>%
filter(cluster == 1) %>%
pull(x2) %>%
mean()
# クラスタ2
x1_c2 <- df2 %>%
filter(cluster == 2) %>%
pull(x1) %>%
mean()
x2_c2 <- df2 %>%
filter(cluster == 2) %>%
pull(x2) %>%
mean()
# クラスタ3
x1_c3 <- df2 %>%
filter(cluster == 3) %>%
pull(x1) %>%
mean()
x2_c3 <- df2 %>%
filter(cluster == 3) %>%
pull(x2) %>%
mean()
## tibble にまとめる
df2_center <- tibble(
x1 = c(x1_c1, x1_c2, x1_c3),
x2 = c(x2_c1, x2_c2, x2_c3),
cluster = factor(1:3))各クラスタの重心を図に加える。
p2_c <- p2 +
geom_point(data = df2_center,
aes(x = x1, y = x2,
color = cluster,
shape = cluster),
shape = 4, size = 4)
plot(p2_c)
各点とそれぞれのクラスタの重心との距離を測る。
df3 <- df2 %>%
mutate(dist_1 = dist_euc(x1, x2, x1_c1, x2_c1),
dist_2 = dist_euc(x1, x2, x1_c2, x2_c2),
dist_3 = dist_euc(x1, x2, x1_c3, x2_c3))各点を、最も距離が短いクラスに分類する。
df3 <- df3 %>%
mutate(cluster = case_when(
dist_1 <= dist_2 & dist_1 < dist_3 ~ 1,
dist_2 < dist_1 & dist_2 <= dist_3 ~ 2,
TRUE ~ 3))途中経過を可視化する。
p3 <- ggplot(df3, aes(x = x1, y = x2,
color = as.factor(cluster),
shape = as.factor(cluster))) +
geom_point() +
theme(legend.position = "none")
plot(p3)
(3)
ここまでの分類結果から、各クラスタの重心の座標を求める。
# クラスタ1
x1_c1 <- df3 %>%
filter(cluster == 1) %>%
pull(x1) %>%
mean()
x2_c1 <- df3 %>%
filter(cluster == 1) %>%
pull(x2) %>%
mean()
# クラスタ2
x1_c2 <- df3 %>%
filter(cluster == 2) %>%
pull(x1) %>%
mean()
x2_c2 <- df3 %>%
filter(cluster == 2) %>%
pull(x2) %>%
mean()
# クラスタ3
x1_c3 <- df3 %>%
filter(cluster == 3) %>%
pull(x1) %>%
mean()
x2_c3 <- df3 %>%
filter(cluster == 3) %>%
pull(x2) %>%
mean()
## tibble にまとめる
df3_center <- tibble(
x1 = c(x1_c1, x1_c2, x1_c3),
x2 = c(x2_c1, x2_c2, x2_c3),
cluster = factor(1:3))各クラスタの重心を図に加える。
p3_c <- p3 +
geom_point(data = df3_center,
aes(x = x1,
y = x2,
color = cluster,
shape = cluster),
shape = 4, size = 4)
plot(p3_c)
各点とそれぞれのクラスタの重心との距離を測る。
df4 <- df3 %>%
mutate(dist_1 = dist_euc(x1, x2, x1_c1, x2_c1),
dist_2 = dist_euc(x1, x2, x1_c2, x2_c2),
dist_3 = dist_euc(x1, x2, x1_c3, x2_c3))各点を、最も距離が短いクラスに分類する。
df4 <- df4 %>%
mutate(cluster = case_when(
dist_1 <= dist_2 & dist_1 < dist_3 ~ 1,
dist_2 < dist_1 & dist_2 <= dist_3 ~ 2,
TRUE ~ 3))途中経過を可視化する。
p4 <- ggplot(df4, aes(x = x1,
y = x2,
color = as.factor(cluster),
shape = as.factor(cluster))) +
geom_point() +
theme(legend.position = "none")
plot(p4)
(4)
ここまでの分類結果から、各クラスタの重心の座標を求める。
# クラスタ1
x1_c1 <- df4 %>%
filter(cluster == 1) %>%
pull(x1) %>%
mean()
x2_c1 <- df4 %>%
filter(cluster == 1) %>%
pull(x2) %>%
mean()
# クラスタ2
x1_c2 <- df4 %>%
filter(cluster == 2) %>%
pull(x1) %>%
mean()
x2_c2 <- df4 %>%
filter(cluster == 2) %>%
pull(x2) %>%
mean()
# クラスタ3
x1_c3 <- df4 %>%
filter(cluster == 3) %>%
pull(x1) %>%
mean()
x2_c3 <- df4 %>%
filter(cluster == 3) %>%
pull(x2) %>%
mean()
## tibble にまとめる
df4_center <- tibble(
x1 = c(x1_c1, x1_c2, x1_c3),
x2 = c(x2_c1, x2_c2, x2_c3),
cluster = factor(1:3))各クラスタの重心を図に加える。
p4_c <- p4 +
geom_point(data = df4_center,
aes(x = x1,
y = x2,
color = cluster,
shape = cluster),
shape = 4, size = 4)
plot(p4_c)
各点とそれぞれのクラスタの重心との距離を測る。
df5 <- df4 %>%
mutate(dist_1 = dist_euc(x1, x2, x1_c1, x2_c1),
dist_2 = dist_euc(x1, x2, x1_c2, x2_c2),
dist_3 = dist_euc(x1, x2, x1_c3, x2_c3))各点を、最も距離が短いクラスに分類する。
df5 <- df5 %>%
mutate(cluster = case_when(
dist_1 <= dist_2 & dist_1 < dist_3 ~ 1,
dist_2 < dist_1 & dist_2 <= dist_3 ~ 2,
TRUE ~ 3))途中経過を可視化する。
p5 <- ggplot(df5, aes(x = x1,
y = x2,
color = as.factor(cluster),
shape = as.factor(cluster))) +
geom_point() +
theme(legend.position = "none")
plot(p5)
p4 と p5 を比べてみよう。
plot(p4 | p5)
グルーピングに変化がないことがわかる。よって、ここでクラスタリングを終了する。
Rによるクラスタリングの実行手順
clustrdパッケージに含まれる macro データを使って実習を行う。
data(macro)
glimpse(macro)## Rows: 20
## Columns: 6
## $ GDP <dbl> 4.8, 3.2, 3.9, 2.3, 3.6, 4.1, 4.1, 2.9, 3.2, 2.3, 2.8, 1.1, 1.4, 1…
## $ LI <dbl> 8.4, 2.5, -1.0, 0.7, 2.5, 1.1, 1.4, 1.6, 0.6, 5.6, -7.5, 0.6, -0.1…
## $ UR <dbl> 8.1, 8.4, 11.8, 11.7, 19.0, 8.9, 4.5, 4.2, 10.3, 3.2, 4.9, 4.7, 9.…
## $ IR <dbl> 5.32, 5.02, 3.60, 3.69, 4.83, 4.20, 5.59, 3.69, 11.70, 20.99, 4.84…
## $ TB <dbl> 0.7, 1.6, 8.8, 3.9, 1.2, 7.0, -1.4, 7.0, -8.3, 0.0, -8.7, -0.6, 4.…
## $ NNS <dbl> 4.7, 5.2, 7.7, 7.3, 9.6, 4.0, 7.0, 15.8, 8.0, 12.7, 14.0, 9.4, 12.…このデータセットは、GDPをはじめとする6つの経済指標を20か国について記録している。データに含まれる国は、row.names() で確認できる。
row.names(macro)## [1] "Australia" "Canada" "Finland" "France" "Spain"
## [6] "Sweden" "USA" "Netherlands" "Greece" "Mexico"
## [11] "Portugal" "Austria" "Belgium" "Denmark" "Germany"
## [16] "Italy" "Japan" "Norway" "Switzerland" "UK"階層的クラスタリング
階層的クラスタリング (hierarchical clustering) はhclust() で実行できる。 method 引数で、グループ間の距離の測り方を指定する。選択肢は、
"single"：単リンク法（最短距離法）"complete"：完全リンク法（最長距離法）"average"：群平均法"centroid"：重心法"ward.D2"：ウォード法
である。
まず、単リンク法でクラスタリングを実行してみよう。hclust() の第1引数は、dist() の結果として出力される非類似度行列（dissimilarity matrix）なので、dist() の結果を保存してから実行しよう。
macro_dist <- dist(macro)
c_single <- hclust(macro_dist, method = "single")この結果を plot() すると、樹形図（dendrogram, デンドログラム）が表示される。
plot(c_single)
各ステップでのグループ間の距離は、樹形図の枝の高さで示されている。 この高さは、次のようにして確認できる。
c_single$height## [1] 2.158333 3.195638 4.375214 4.389032 4.426839 4.505830 5.050980
## [8] 5.138959 5.206726 5.444089 5.458470 5.490902 5.591386 5.934172
## [15] 6.516134 8.491148 11.001459 12.502400 15.850050この結果を元に、クラスタの数を\(k\)に指定するにはcutree()を使う。\(k=3\)にしてみよう。
cutree(c_single, k = 3)## Australia Canada Finland France Spain Sweden
## 1 1 1 1 1 1
## USA Netherlands Greece Mexico Portugal Austria
## 1 1 1 2 3 1
## Belgium Denmark Germany Italy Japan Norway
## 1 1 1 1 1 1
## Switzerland UK
## 1 1メキシコとポルトガルがそれぞれ単独でグループを作り、残りの18か国が1つのグループに分けられれた。
次に、完全リンク法でクラスタリングを行ってみよう。
c_complete <- hclust(macro_dist, method = "complete")樹形図を示す。
plot(c_complete)
\(k=3\)でクラスタに分ける。
cutree(c_complete, k = 3)## Australia Canada Finland France Spain Sweden
## 1 1 1 1 1 1
## USA Netherlands Greece Mexico Portugal Austria
## 1 1 2 3 2 1
## Belgium Denmark Germany Italy Japan Norway
## 1 1 1 1 1 1
## Switzerland UK
## 1 1メキシコは相変わらず単独でクラスタを構成しているが、ポルトガルはギリシャと同じクラスタになった。残りの17か国は同じクラスタに分類されている。
実習課題
- 距離の測り方を変えてクラスタリングを実行し、距離の測り方によって結果がどのように異なるか考察しなさい。
- \(k\)の値を4や5にすると、距離の測り方による違いはどう変化するだろうか。
非階層的クラスタリング：\(K\)平均法
\(K\)平均法は、kmeans() で実行できる。クラスタの数は centers 引数で指定する。 \(K=3\)で分類してみよう。
set.seed(2021-11-07)
c_k3 <- kmeans(macro, centers = 3)結果を確認しよう。
c_k3$cluster## Australia Canada Finland France Spain Sweden
## 2 2 2 2 2 2
## USA Netherlands Greece Mexico Portugal Austria
## 3 1 3 3 3 3
## Belgium Denmark Germany Italy Japan Norway
## 2 2 2 2 1 1
## Switzerland UK
## 1 3階層的クラスタリングとは異なる結果が得られている。
kmeans() 関数は、最初のステップにおけるグループの中心をランダムに選ぶ。 nstart 引数を1より大きい値に設定すると、複数の初期値をランダムに選び、最もよい結果のみを返す。 試しに、nstart = 50 にしてみよう。
c_k3b <- kmeans(macro, centers = 3, nstart = 50)結果を確認しよう。
c_k3b$cluster## Australia Canada Finland France Spain Sweden
## 2 2 2 2 2 2
## USA Netherlands Greece Mexico Portugal Austria
## 2 3 1 1 1 2
## Belgium Denmark Germany Italy Japan Norway
## 3 2 2 2 3 3
## Switzerland UK
## 3 2先ほどとは異なる結果が得られた。ただし、グループのラベルの付け方が異なる。（注意：日本やノルウェーなどを含むグループのラベルが1から3に変わっているが、ラベルの変更は結果の違いとは言えない。）
各グループのばらつき具合は、クラスタ内総平方和で測られる。先程の結果と比べてみよう。
rbind(c_k3$tot.withinss,
c_k3b$tot.withinss)## [,1]
## [1,] 895.9534
## [2,] 863.4262クラスタ内総平方和が小さくなっており、先程よりも各グループが似た者同士の集まりになっていると考えられる。
結果を可視化するために、元のデータにクラスタリングの結果を加えた tibble を作る。
df_k3b <- macro %>%
mutate(cluster = factor(c_k3b$cluster))GGallyパッケージの ggpairs() で散布図行列を作る。
ggpairs(df_k3b,
columns = 1:6,
aes(color = cluster),
legend = 1) +
theme(legend.position = "bottom")
、