Rによるデータサイエンス（08/21） - 「大人の教養・知識・気付き」を伸ばすブログ

　Rについて

Rによるデータサイエンス(第2版):データ解析の基礎から最新手法まで

作者:金明哲
森北出版

Amazon

をベースに学んでいく。
　今回はクラスター分析（PP.106-125）を扱う。

前回
9.　クラスター分析
参考文献
補足　スペック情報
次回

前回

power-of-awareness.com

9.　クラスター分析

　データ構造の側面から似ている個体(変数)を同じグループに仕分けて分類を行う方法のうち、データのパターンが似ている個体を同じグループにまとめる分析手法がクラスター分析である。どの個体がどのグループに属するかに関する事前情報を用いずにグルーピングする教師無し学習の1つである。
　クラスタリングは類似性という基準に基づいてデータをグループ化する手法であり、個体間の関係を理解するのに有用な手法である。

9.1　階層的クラスター分析

　個体間（とクラスター間）の類似度（または距離）に基づいて最も似ている個体から順に集めてクラスター(集団)を作成する方法を階層的クラスター分析という。階層的クラスター分析においてクラスターがつくられる様子は樹形図（デンドログラム）と呼ぶ。

	(1)	データから距離(または類似度)を計算する。
	(2)	クラスター分析の手法を選択する。
	(3)	選択された方法のコーフェン(cophenetic)行列を求める。
	(4)	コーフェン行列に基づいて樹形図を作成する。
	(5)	結果を検討する。

　階層的クラスター分析は、データから距離行列を求め、距離行列からcophenetic行列を求める。そしてcophenetic行列に基づいて樹形図を描く。

データ行列	$\Rightarrow$	距離行列	$\Rightarrow$	cophenetic行列	$\Rightarrow$	結果
$\begin{aligned}\begin{bmatrix}x_{11}&x_{12}&\cdots&x_{1n}\\x_{21}&x_{22}&\cdots&x_{2n}\\\vdots&\vdots&\ddots&\vdots\\x_{m1}&x_{m2}&\cdots&x_{mn}\end{bmatrix}\end{aligned}$	$\Rightarrow$	$\begin{aligned}\begin{bmatrix}0&&&\\d_{21}&0&&\\\vdots&\vdots&\ddots&\\d_{m1}&d_{m2}&\cdots&0\end{bmatrix}\end{aligned}$	$\Rightarrow$	$\begin{aligned}\begin{bmatrix}0&&&\\c_{21}&0&&\\\vdots&\vdots&\ddots&\\c_{m1}&c_{m2}&\cdots&0\end{bmatrix}\end{aligned}$	$\Rightarrow$	樹形図を描く

　距離行列からcophenetic行列を生成する方法はいくつかある。第1段階はすべて同じで、最も距離が近い2つの個体間の距離をcophenetic距離とする。この後新たにcophenetic距離を求める方法はクラスター分析の方法によって異なる。

9.1.1　個体間の類似度

　クラスター分析は個体間の類似度を「距離」を用いて測る*1。

9.1.2　クラスター間距離

　データ集合間における任意のデータ間の距離を求めた後、最も距離の小さい個体同士で1つのクラスターを形成する。するとクラスターと各個体との距離およびクラスター間の距離を定義する必要がある。
　いま観測された $n$ 個体の $p$ 次元データを $S=\{\boldsymbol{x}_1,\cdots,\boldsymbol{x}_n\}$ とし、この $S$ を分類対象とする。データ集合に含まれる任意の個体 $i,j\in S,i\neq j$ に対して

$\begin{aligned} d_{ij}=d(\boldsymbol{x}_i,\boldsymbol{x}_j),i\neq j,i,j\in\{1,2,\cdots,n\} \end{aligned}$

とする。この個体間距離に基づいて形成された2つのクラスターを $C_1,C_2$ とする。このとき

(1)	最短距離法	2つのクラスターに含まれる任意の個体間の距離のうち最も小さいものをそのクラスター間距離とする。すなわち $\begin{aligned}d(C_1,C_2)=\displaystyle{\min_{i,j}\left\{d(\boldsymbol{x}_i^{(1)},\boldsymbol{x}_j^{(2)})\right\}}\end{aligned}$ とする。
(2)	最長距離法	2つのクラスターに含まれる任意の個体間の距離のうち最も大きいものをそのクラスター間距離とする。すなわち $\begin{aligned}d(C_1,C_2)=\displaystyle{\max_{i,j}\left\{d(\boldsymbol{x}_i^{(1)},\boldsymbol{x}_j^{(2)})\right\}}\end{aligned}$ とする。
(3)	群平均法	2つのクラスターに含まれるすべての個体間距離の平均値をクラスター間の距離とする。すなわち $\begin{aligned}d(C_1,C_2)&=\displaystyle{\frac{1}{n_1 n_2}\sum_{i=1}^{n_1}\sum_{i=1}^{n_2}d(\boldsymbol{x}_i^{(1)},\boldsymbol{x}_j^{(2)})},\\\boldsymbol{x}_i^{(1)}&\in C_1,\boldsymbol{x}_j^{(2)}\in C_2\end{aligned}$ とする。
(4)	重心法	クラスターの標本平均ベクトル間の距離をクラスター間の距離とする。すなわち $\begin{aligned}d(C_1,C_2)&=d(\bar{\boldsymbol{x}}_1,\bar{\boldsymbol{x}}_2),\\\bar{\boldsymbol{x}}_k&=\displaystyle{\frac{1}{n_k}\sum_{i=1}^{n_k}\boldsymbol{x}_i^{(k)}},k=1,2\end{aligned}$ である。
(5)	メディアン法	重心法と同様に重心を決定する。ただし合成したクラスターの重心を中点とする。

といった方法がある。最短距離法や最長距離法、群平均法では任意の距離を利用できる一方で、重心法やメディアン法では(平方)Euclid距離に限られる。

9.1.3　Ward法

　クラスタリングにおいてデータの散らばりを定量的に測る指標に基づいて逐次クラスターを形成する分類法で(平方)Euclid距離に基づく方法をWard法という。

	(1)	各クラスター内のデータが1個*2の段階において散らばりの尺度 $W=0$ とする。
	(2)	散らばりの尺度 $W$ が最も小さくなるような2つのクラスターを融合させて、クラスターの数を1つ減らす。
	(3)	クラスターの数が1つならば終了させる。そうでなければ(2)に戻る。

　1次元データ集合 $S=\{x_1,\cdots,x_n\}$ をWard法で融合していくことを考える場合、偏差平方和

$\begin{aligned} \displaystyle{\sum_{i=1}^{n}(x_i\bar{x})^2},\bar{x}=\displaystyle{\frac{1}{n}\sum_{i=1}^{n}x_i} \end{aligned}$

を用いる。いま $S$ が $m$ 個のクラスター

$\begin{aligned} C_1=\{x_{1}^{(1)},\cdots,x_{n_1}^{(1)}\},\cdots,C_m=\{x_{1}^{(m)},\cdots,x_{n_m}^{(m)}\} \end{aligned}$

に分割されたとする。このとき各クラスター内の偏差平方和を、クラスター $j$ に所属する個体数を $n_j$ とすれば、

$\begin{aligned} W_j=\displaystyle{\sum_{i=1}^{n_j}(x_i^{(j)}-\bar{x}_j)^2},\bar{x}_j=\displaystyle{\frac{1}{n_j}\sum_{i=1}^{n_j}x_{i}^{(j)}} \end{aligned}$

で計算し、クラスター内全偏差平方和を

$\begin{aligned} W=\displaystyle{\sum_{k=1}^{m}W_k} \end{aligned}$

で定義する。ウォード法では2つのクラスター

$\begin{aligned} C_1=\{x_1^{(1)},\cdots,x_{n_1}^{(1)}\},C_2=\{x_1^{(2)},\cdots,x_{n_1}^{(2)}\} \end{aligned}$

があったときにこれらを融合した $C_{12}=C_1\cup C_2$ の偏差平方和を $W_{12}$ として

$\begin{aligned} d(C_1,C_2)=W_{12}-W_1-W_2 \end{aligned}$

で融合の水準を定義する。この融合の水準が最も小さいものを融合させる。
　分類対象が $n$ 個の $p$ 次元データ集合 $\{\boldsymbol{x}_1,\cdots,\boldsymbol{x}_n\}$ からなるときは、基本的に偏差平方和行列

$\begin{aligned} S&=(s_{jk})=\displaystyle{\sum_{i=1}^{n}(\boldsymbol{x}_i-\bar{\boldsymbol{x}})(\boldsymbol{x}_i-\bar{\boldsymbol{x}})^{\prime}},\\ \bar{\boldsymbol{x}}&=\displaystyle{\frac{1}{n}\sum_{i=1}^{n}\boldsymbol{x}_i} \end{aligned}$

を用いて定義する。このとき

$\begin{aligned} \mathrm{tr}(S)=\displaystyle{\sum_{i=1}^{n}s_{ii}},\left|S\right| \end{aligned}$

のいずれかで測ることができる。前者は偏差平方和の総和であり、後者は一般化分散と呼ばれる。
　いま $S$ が $m$ 個のクラスター

$\begin{aligned} C_1=\{\boldsymbol{x}_{1}^{(1)},\cdots,\boldsymbol{x}_{n_1}^{(1)}\},\cdots,C_m=\{\boldsymbol{x}_{1}^{(m)},\cdots,\boldsymbol{x}_{n_m}^{(m)}\} \end{aligned}$

に分割されたとする。このとき各クラスター内の偏差平方和行列を $S_j$ とすれば、

$\begin{aligned} W_j=\begin{cases}\mathrm{tr}(S_j),\\ \left|S_j\right| \end{cases} \end{aligned}$

で散らばりを計算し、クラスター内全体の散らばりを

$\begin{aligned} W=\begin{cases} \displaystyle{\sum_{j=1}^{m}\mathrm{tr}(S_j)},\\ \displaystyle{\sum_{j=1}^{m}\left|S_j\right|} \end{cases} \end{aligned}$

で定義する。各クラスターが1つのデータから成る初期状態では $W=0$ とする。
　2つのクラスター

$\begin{aligned} C_1=\{\boldsymbol{x}_1^{(1)},\cdots,\boldsymbol{x}_{n_1}^{(1)}\},C_2=\{\boldsymbol{x}_1^{(2)},\cdots,\boldsymbol{x}_{n_1}^{(2)}\} \end{aligned}$

があったときにこれらを融合した $C_{12}=C_1\cup C_2$ の偏差平方和を $W_{12}$ として

$\begin{aligned} d(C_1,C_2)=W_{12}-W_1-W_2=\begin{cases} \mathrm{tr}(S_{12})-\mathrm{tr}(S_{1})-\mathrm{tr}(S_{2}),\\ \left|S_{12}\right|-\left|S_1\right|-\left|S_2\right| \end{cases} \end{aligned}$

で融合の水準を定義する。この融合の水準が最も小さいものを融合させる。

9.2　非階層的クラスター分析

　階層的クラスター分析は個体数が多いと計算量が膨大になるため、大量のデータの解析には向いていない。大規模データには非階層的クラスター分析が用いられる。代表的なのは $k$ - $\mathrm{means}$ 法である。

	(1)	$k$ 個のクラスターに対して中心の初期値を適当に与える。
	(2)	すべてのデータについて $k$ 個のクラスター中心との距離を求め、各データを最も近いクラスターに分類する。
	(3)	形成されたクラスターの中心を求め直す。
	(4)	(3)において計算し直した中心がその前の値から変わらないのであれば終了、そうでなければ(2)に戻る。

9.3　モデルに基づいたクラスター分析

　混合分布によるクラスター分析、潜在クラスター分析と呼ばれる。観測データが異なる確率分布による混合分布に従うと仮定し、個体が属するクラスのラベルを隠れた変数として推定する。理論上は確率分布の指定が無いものの、現状では楕円分布（正規分布をより一般化したもの）を用いる。

9.3.1　混合正規モデルを用いたクラスター分析

　混合比率 $0\leq\pi_j\leq1,\displaystyle{\sum_{j=1}^{g}\pi_j}=1$ を用いて

$\begin{aligned} f(x|\boldsymbol{\theta})=\displaystyle{\sum_{j=1}^{g}\pi_j \frac{1}{\sqrt{2\pi\sigma_j^2}}\exp\left\{-\frac{(x-\mu_j)^2}{2\sigma_j^2}\right\}} \end{aligned}$

という確率モデルを考えたとき、これを混合正規モデルという。より一般に複数の確率分布の線形結合で表される確率分布を混合分布モデルという。
　 $\mathrm{Bayes}$ 統計学の枠組みに基づくと、混合分布モデルに基づきデータの分類ができる。第 $j$ 番目の群 $G_j$ を特徴づける混合分布 $f_j(\boldsymbol{x}|\boldsymbol{\theta}_j)$ とする。これはデータ $\boldsymbol{x}$ が第 $j$ 群から生成されたとする場合の条件付き尤度と言える。また対応する混合比率 $\pi_j$ は、その群を選択する事前確率と解釈できる。そこで $\mathrm{Bayes}$ の定理から、データ $\boldsymbol{x}$ が群 $G_j$ から生成されたという事後確率は

$\begin{aligned} P\{j|\boldsymbol{x}\}=\displaystyle{\frac{\pi_j f_j(\boldsymbol{x}|\boldsymbol{\theta}_j)}{\displaystyle{\sum_{k=1}^{g}\pi_k f_k(\boldsymbol{x}|\boldsymbol{\theta}_k)}}}=\displaystyle{\frac{\pi_j f_j(\boldsymbol{x}|\boldsymbol{\theta}_j)}{f(\boldsymbol{x}|\boldsymbol{\theta})}} \end{aligned}$

で与えられる。そこで、混合分布モデルが推定できれば、 $g$ 個の事後確率の中でその値が最も大きい群にデータ $\boldsymbol{x}$ に属するとする。このように分類の対象とするデータでモデルを推定し、事後確率最大化により各データを分類すればよい。

9.3.2　EMアルゴリズムによるモデル推定

　混合分布モデルの母数推定には $\mathrm{EM}$ アルゴリズムをはじめとする数値的な最適化法を用いる。
　混合比率、平均ベクトルおよび分散共分散行列 $\{(\pi_j,\boldsymbol{\mu}_j,\Sigma_j);j=1,2,\cdots,g\}$ に対して初期値を設定する。また事前に収束判定条件として $\varepsilon\gt0$ を定めておく。

$t=1$ とし、初期値 $\left(\pi_j^{(1)},\boldsymbol{\mu}_j^{(1)},\Sigma_j^{(1)}\right)$ を設定する。

$i$ 番目のデータ $\boldsymbol{x}_i$ が群 $G_j$ から生成されたという事後確率
$\begin{aligned}p^{(t+1)}(j|\boldsymbol{x}_i)=\displaystyle{\frac{\pi_j^{(t)}f_j(\boldsymbol{x}_i|\boldsymbol{\theta}_j^{(t)})}{\displaystyle{\sum_{k=1}^{g}\pi_k^{(t)}f_k(\boldsymbol{x}_i|\boldsymbol{\theta}_k^{(t)})}}}\end{aligned}$
を計算する。

$t=t+1$ とし、それに対応する混合比率、平均ベクトルおよび分散共分散行列を
$\begin{aligned}\pi_{j}^{(t+1)}&=\displaystyle{\frac{1}{n}\sum_{i=1}^{n}p^{(t+1)}(j|\boldsymbol{x}_i)},\\\boldsymbol{\mu}_{j}^{(t+1)}&=\displaystyle{\frac{1}{n\pi_j^{(t+1)}}\sum_{i=1}^{n}p^{(t+1)}(j|\boldsymbol{x}_i)\boldsymbol{x}_i},\\\Sigma_{j}^{(t+1)}&=\displaystyle{\frac{1}{n\pi_j^{(t+1)}}\sum_{i=1}^{n}p^{(t+1)}(j|\boldsymbol{x}_i)\left(\boldsymbol{x}_i-\boldsymbol{\mu}_j^{(t+1)}\right)\left(\boldsymbol{x}_i-\boldsymbol{\mu}_j^{(t+1)}\right)^{\prime}}\end{aligned}$
で更新する。

観測データに基づく尤度関数について
$\begin{aligned}\left|\displaystyle{\sum_{i=1}^{n}}\log\left\{\sum_{j=1}^{g}\pi_{j}^{(t+1)}f_j\left(\boldsymbol{x}_i\left|\right.\boldsymbol{\mu}_j^{(t+1)},\Sigma_j^{(t+1)}\right)\right\}-\displaystyle{\sum_{i=1}^{n}}\log\left\{\sum_{j=1}^{g}\pi_{j}^{(t)}f_j\left(\boldsymbol{x}_i\left|\right.\boldsymbol{\mu}_j^{(t)},\Sigma_j^{(t)}\right)\right\}\right|\end{aligned}$
が $\varepsilon$ 以下ならばステップを終了する。そうでなければ2.に戻る。

　ここまで、混合分布モデルの成分数 $g$ を固定した上でモデルを推定してきた。この $g$ がクラスター数に相当する。その決定は、 $g=2,3,\cdots$ として赤池情報量基準 $\mathrm{AIC}$ やベイズ情報量基準 $\mathrm{BIC}$ によって評価・選択する方法がある。

9.4　Rでのケーススタディ

######################
### クラスター分析 ### 
######################
library(ggplot2)
library(MASS)

# データの準備
df_scores <- data.frame("arithmetics" = c(89,57,80,40,78,55,90),
                        "science" = c(80,70,90,60,85,65,85),
                        "japanese" = c(67,80,35,50,45,80,88),
                        "english" = c(46,85,40,45,55,75,92),
                        "socialstudy" = c(50,90,50,55,60,85,95))
row.names(df_scores) <- c("田中","佐藤","鈴木","本田","川端","吉野","斎藤")

### 階層的クラスター分析
df_scores_dist <- dist(df_scores, method = "euclid")
ls_clust <- hclust(d = df_scores_dist, method = "ward.D")

summary(ls_clust)

#クラスター形成の履歴: (個体数-1)行2列
# 行番号：クラスター形成の順番
ls_clust$merge

# 樹形図の枝の長さ
ls_clust$height

# デンドログラムを図示
plot(ls_clust)
plot(ls_clust,hang = -1)
rect.hclust(ls_clust, k =2, border = c(2,4))

cutree(tree = ls_clust,k = 2) # 所属するクラスター

# ヒートマップでの各クラスの特徴分析
heatmap(as.matrix(df_scores))

# cophenetic行列の確認
cophenetic(ls_clust)

# 妥当性確認：距離行列とcophenetic行列の相関係数を見る
cor(df_scores_dist,cophenetic(ls_clust))

### 非階層的クラスター分析
# k-means法
km_scores <- kmeans(df_scores,2)
summary(km_scores)

km_scores$cluster # 分類されたクラスター


### モデルに基づいたクラスター分析
library("mclust")

data(iris)

# まずは指定するクラスター数を推計する
mclust_BIC_iris <- mclustBIC(iris[,-5])
plot(mclust_BIC_iris) # mclustではBICが「大きい」方を良いものとする=通常と逆

colnames(mclust_BIC_iris)[which(mclust_BIC_iris[3,]==max(mclust_BIC_iris[3,],na.rm = F))] # NoC＝3で一番BICが大きいモデル

# ただしhcにVEVが無いため、2番目を用いる
iris_hc <- hc(modelName = colnames(mclust_BIC_iris)[which(mclust_BIC_iris[3,]==max(mclust_BIC_iris[3,-12],na.rm = F))],
              data = iris[,-5])
iris_hcl <- hclass(iris_hc,3)

table(iris[,5],iris_hcl)

clPairs(data = iris[,-5],cl = iris_hcl)

混合正規分布( $\mathrm{VVV},k=3$ )に基づいた $\mathrm{iris}$ のクラスタリング結果の対散布図

参考文献

小西貞則(2010)「多変量解析入門　線形から非線形へ」（岩波書店）*3

補足　スペック情報

エディション	Windows 10 Home
バージョン	20H2
プロセッサ	Intel(R) Core(TM) i5-1035G4 CPU @ 1.10GHz 1.50 GHz
実装 RAM	8.00 GB
システムの種類	64 ビットオペレーティングシステム、x64 ベースプロセッサ
R　バージョン	3.6.3 (2020-02-29)
RStudio　バージョン	1.2.5033

次回

power-of-awareness.com

*1:距離については前回議論したのでここでは触れない。power-of-awareness.com

*2:一度も融合を行なっていない状態

*3:多変量解析入門 - 岩波書店参照。

前回

9. クラスター分析

9.1 階層的クラスター分析

9.1.1 個体間の類似度

9.1.2 クラスター間距離

9.1.3 Ward法

9.2 非階層的クラスター分析

9.3 モデルに基づいたクラスター分析

9.3.1 混合正規モデルを用いたクラスター分析

9.3.2 EMアルゴリズムによるモデル推定

9.4 Rでのケーススタディ