PCA, Principal Component Analysis（主成分分析）は、線形独立であるか自明でない多数の変数を、ばらつきが大きい順から線形独立な主成分 (PC, Principal Component ) と呼ばれる変数に変換（合成）します。これは本記事のシリーズの最初で説明したことです [1]。
主成分で表す空間は直交空間であるため、主成分は互いに独立しています。また変動（ばらつき）が大きい順に PC1, PC2, PC3 のように主成分の番号が割り当てられています。データ集団の差を調べる際に、PC1 から順に調べていくことは、変動の大きな、すなわち差が大きな要因から調べていくことになるので合理的です。
ある主成分におけるデータ集団の差を調べるということは、結局のところ元データでどのパラメータ（変数）がその差に影響を及ぼしているかを調べることになります。データセット iris ではパラメータが「がく」と「花弁」の長さと幅、計４種類しかありません。このような小規模な解析では大げさですが、主成分のプロットから二種類の差異がある集団を選び、その差が生じる内訳を、主成分分析する前のパラメータ（この例では「がく」と「花弁」の長さと幅）の重要度（寄与度）で表してみます（特徴選択）。
R スクリプト
まず、本記事で使用したスクリプトを以下に示しました。今回はパラメータの重要度を計算するために Random Forest のアルゴリズムを使用しています。
library(ggplot2) library(randomForest) # classification function classify <- function(x) { if (x[1] == "setosa") "A" else "B" } # ----------------------------------------------------------------------------- # 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) df_out$Class <- sapply(df_out$Species, function(x) classify(x)) # scatter by Species p <- ggplot(df_out, aes(x = PC1, y = PC2, color = Species)) p <- p + geom_point() plot(p) # scatter by Class p <- ggplot(df_out, aes(x = PC1, y = PC2, color = Class)) p <- p + scale_color_manual(values=c("#ff0000", "#0000ff")) p <- p + geom_point() plot(p) # data preparation for classification with Random Forest method df_class <- iris df_class$Class <- sapply(df_class$Species, function(x) classify(x)) df_class$Class <- as.factor(df_class$Class) df_class$Species <- NULL head(df_class) # Random Forest set.seed(71) df_rf <- randomForest(Class ~ ., data = df_class, 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)
集団のクラス分け
主成分 PC1 と PC2 の散布図（「種 (Species)」別に色分け）
以下は主成分分析した結果から主成分 PC1 と PC2 で相関図を描き、「種 (Species)」別に色分けしたものです。PC1 の軸からデータの集団を眺めると、大雑把に、赤の setosa の集団と、緑の versicolor、青の virginica を合わせた集団の、２つのグループを形成していると見て取れます。
主成分 PC1 と PC2 の散布図（「クラス (Class)」別に色分け）
そこで、いったん、「種」別の色分けを外し、上記二つの集団をあらためて、赤 (A) と、青 (B) で分類し直します。データセット iris は 150 行ほどの小さなデータ量ですので、以下のように全てのデータを使ってクラス分けをしました。
- 赤を全ての setosa
- 青を versicolor と virginica を合わせた全て
> head(df_class) Sepal.Length Sepal.Width Petal.Length Petal.Width Class 1 5.1 3.5 1.4 0.2 A 2 4.9 3.0 1.4 0.2 A 3 4.7 3.2 1.3 0.2 A 4 4.6 3.1 1.5 0.2 A 5 5.0 3.6 1.4 0.2 A 6 5.4 3.9 1.7 0.4 A
あらためて、主成分 PC1 と PC2 で相関図を描き、赤(A)と青(B)で色付けし直しました。
Random Forest で特徴量を算出
A と B の違いに対する、パラメータ（「がく」と「花弁」の長さと幅）の重要度（寄与度）を Random Forest で計算すると、以下のようになります。
> round(importance(df_rf), 2)
A B MeanDecreaseAccuracy MeanDecreaseGini
Sepal.Length 2.96 2.63 3.95 7.12
Sepal.Width 3.33 2.49 3.62 0.39
Petal.Length 21.23 20.82 21.58 28.85
Petal.Width 22.22 21.78 22.72 29.83
以下に、各パラメータの MeanDecreaseAccuracy の値をパラメータの重要度とみなして棒グラフで表しました。主成分 PC1 におけるクラス A と B の差は、主に petal（花弁）の長さと幅の違いによるものであると結論できるでしょう。
まとめ
データセット iris の解析では、主成分分析や機械学習のアルゴリズムを用いずとも、丁寧にパラメータ間を調べるプロットを作成して注意深く眺めれば、「種」の差とパラメータ（「がく」と「花弁」の長さと幅）の関係を見つけることができるでしょう。
しかし、数百、数千、あるいはもっと多くのパラメータを持ったデータを解析する場合はどうでしょう。たくさんのプロットを描いて目視でカギとなる特徴を探し出し、その関係を理解するのは不可能に近い作業になります。データ集団のバラツキが大きい部分から軸を選び直し、主成分という名前のパラメータへ変換して主成分の直交空間を作り、その主成分間の関係を順番に解析して、差異をもたらす事象に関与する生パラメータを炙り出していく作業は、回り道をしているようですが計算の大部分をコンピュータに任せられるため、大量のデータを解析するのに向いています。
参考サイト
- bitWalk's: R で主成分分析 (1) [2018-07-10]
- bitWalk's: R で主成分分析 (2) [2018-07-21]
にほんブログ村
0 件のコメント:
コメントを投稿