PCA, Principal Component Analysis（主成分分析）の解析例として、前回まで iris のデータセットを使ってきました [1][2][3]。今回もひきつづき iris のデータセットを使って主成分分析をしますが、クラスター分析を利用して解析する集団を定義してみます。
R スクリプト
まず、本記事で使用したスクリプトを以下に示しました。
library(ggplot2) library(randomForest) library(factoextra) library(plyr) library(reshape2) # ----------------------------------------------------------------------------- # MAIN # ----------------------------------------------------------------------------- # data preparation for PCA df <- as.data.frame(iris) row.names(df) <- paste(df$Species, row.names(df), sep="_") df$Species <- NULL # ----------------------------------------------------------------------------- # PCA df_pca <- prcomp(df) df_out <- as.data.frame(df_pca$x) df_out$Species <- sapply(strsplit(as.character(row.names(df_out)), "_"), "[[", 1) # scatter by Species p <- ggplot(df_out, aes(x = PC1, y = PC2, color = Species)) p <- p + geom_point() plot(p) # ----------------------------------------------------------------------------- # Hierarchical clustering d <- dist(df_out[, 1:4], method = "euclidean") res.hc <- hclust(d, method = "ward.D2" ) plot(res.hc, cex = 0.5, hang = -1) grp <- cutree(res.hc, k = 3) table(grp) fviz_cluster(list(data = df_out[, 1:4], cluster = grp)) # apply cluster to PCA field df_clust <- df_out[, 1:4] df_clust$Cluster <- grp[row.names(df_clust)] df_clust$Cluster <- as.factor(df_clust$Cluster) p <- ggplot(df_clust, aes(x = PC1, y = PC2, color = Cluster)) p <- p + scale_color_manual(values=c("#ff0000", "#00a000", "#0000ff")) p <- p + geom_point() plot(p) # apply cluster to raw data df2 <- df df2$Cluster <- grp[row.names(df_clust)] head(df2) # choose cluster 1 and 3 df2_1_3 <- df2[df2$Cluster != 2, ] df2_1_3$Cluster <- as.factor(df2_1_3$Cluster) # Random Forest set.seed(71) df_rf <- randomForest(Cluster ~ ., data = df2_1_3, importance = TRUE, proximity = TRUE) round(importance(df_rf), 2) # variable importance df_bar <- as.data.frame(importance(df_rf)) sortlist <- order(df_bar$MeanDecreaseAccuracy) df_bar <- df_bar[sortlist,] df_bar$parameter <- row.names(df_bar) p <- ggplot(data = df_bar, aes(x = parameter, y = MeanDecreaseAccuracy)) p <- p + geom_bar(stat = "identity", fill="steelblue") p <- p + scale_x_discrete(limits = row.names(df_bar)) p <- p + coord_flip() plot(p) # boxplot dfmelt <- melt(df2_1_3, measure.vars = 1:4) p <- ggplot(dfmelt, aes(Cluster, value, fill = variable)) p <- p + geom_jitter(aes(color = variable)) + geom_boxplot() p <- p + facet_wrap(~variable, ncol = 4) plot(p)
視覚的なクラス分けの問題点
前回は、単純に setosa と versicolor + virginica のグループに分けて、クラス A と B としました。はたしてこれは妥当だったのでしょうか。
「妥当か」というファジィな問いに対して、満足できる答えを出すことは難しいでしょう。
しかし PC1 から PC4 まである四次元の空間のうち、主要なばらつき軸である PC1 と PC2 の空間で明らかに異なる集団である二集団をみているのだから、そこそこ妥当であろうという人がいるかもしれません。あるいは「種」の情報を前提にして setosa と versicolor + virginica という二つの集団に分けることはあまりに意図的すぎるため、当たり前の結果しか生まないという意見が出るかもしれません。
なんにせよ、あまり人に依存せず、できることなら再現性があるクラス分けをしたいものです。
階層型クラスタリングの利用
データセット iris を主成分分析して得た主成分の空間（PCA 空間と呼ぶことにします）にあるデータに対して、hclust で階層型クラスター分析をおこなって、得られたデンドログラムから、クラスター数を決めます。そこから２つのデータ集団を定義して違いを調べることにします。
以下が hclust によるクラスタリングの結果です。この樹形図状のプロットをみると、大きくは赤の破線で示した三つの集団に分けられます。今回の解析例では大きな集団を見るだけで十分ですのでこの結果からクラスター数を 3 として解析を進めます。
3 つのクラスターをプロットしました。赤 (1) の集団は setosa、緑 (2) は versicolor と virginica が混ざった集団、最後の青 (3) は virginica だけの集団になっています。
Species の代わりに、クラスター分析で得た 3 つのクラスターでデータを色分けした、PC1 と PC2 の散布図のプロットを示しました。
ここではクラスター 1（赤）とクラスター 3（青）の違いを調べてみます。クラスターの情報を元のデータに対応させます。
> head(df2)
Sepal.Length Sepal.Width Petal.Length Petal.Width Cluster
setosa_1 5.1 3.5 1.4 0.2 1
setosa_2 4.9 3.0 1.4 0.2 1
setosa_3 4.7 3.2 1.3 0.2 1
setosa_4 4.6 3.1 1.5 0.2 1
setosa_5 5.0 3.6 1.4 0.2 1
setosa_6 5.4 3.9 1.7 0.4 1
Random Forest で特徴量を算出
クラスター ２ を除いたデータフレームで、クラスター 1 と 3 の違いに対する、パラメータ（「がく」と「花弁」の長さと幅）の重要度（寄与度）を Random Forest で計算すると、以下のようになります。
> round(importance(df_rf), 2)
1 3 MeanDecreaseAccuracy MeanDecreaseGini
Sepal.Length 15.84 16.01 16.12 14.70
Sepal.Width 0.00 0.00 0.00 0.00
Petal.Length 14.39 14.51 14.65 12.87
Petal.Width 15.17 15.39 15.51 13.93
以下に、各パラメータの MeanDecreaseAccuracy の値をパラメータの重要度とみなして棒グラフで表しました。主成分 PC1 におけるクラスター 1 と 3 の差は、sepal（がく）の幅を除いた他のパラメータ、すなわち、がくの長さ、花弁の長さと幅の違いによるものであると結論できるでしょう。
前回の結果と異なる
クラスター 1 は全ての setosa のデータに対応しているので、前回の解析で定義した Class A と同じ集団です。一方、クラスター 3 は、前回の Class B に対して、部分集合になります。そのため、特徴量については同じような傾向を持った結果が出るだろうと予想していたのですが、そうではありませんでした。
ボックスプロットにして、クラスター 1 と 3 における各パラメータの値（分布）を比較すると、Random Forest の導き出した結果が（ファジィですが）妥当であることを確認できます。パラメータの数が多い時には全部を確認できませんので、重要度の大きなパラメータに絞って確認することになります。
まとめ
大量のパラメータを持った大きなデータが、ある一定条件では同じ分布になる、というのが理想状態ですが、現実はそうではありません。差が生じる原因を探索するため、日々データと格闘しています。解析手順が同じであれば、同じデータセットに対し、誰がやっても同じ結果が出て欲しいと常々願っているのですが、確固たる解析手順や方法があるほどに枯れた分野ではないので、それが解析に主観を加えることにつながり、結果として人が違えば結論が異なるという事態になります。
今回は、少しでも解析者の主観を排除できそうな工夫を加えてみました。
本ブログは、自分の考え方を整理するため、簡単な例に置き換え書き留めておくためのメモのようなものです。あとになって自分が読み直すことを想定していますので、内容に誤りがあれば随時修正するようにしています。
参考サイト
- bitWalk's: R で主成分分析 (1) [2018-07-10]
- bitWalk's: R で主成分分析 (2) [2018-07-21]
- bitWalk's: R で主成分分析 (3) [2018-07-23]
にほんブログ村
