2018-07-10

R で主成分分析 (1)

PCA, Principal Component Analysis(主成分分析)は、線形独立であるか自明でない多数の変数を、ばらつきが大きい順から線形独立な主成分 (PC, Principal Component ) と呼ばれる変数に変換(合成)します。データの次元を削減するために用いられます。

次元が削減されてなにが嬉しいかというと、それは、多数の変数=多次元の情報を、少ない次元で端的に表現できるということです。主成分はデータのばらつき=変動の大きい順に PC1, PC2, ... と変換されています。そのため、主成分のいくつかだけでデータのふるまい(特徴)を概ね説明できてしまう場合があります。

最近、PCA を駆使して何百、何千という変数を持つデータの解析をする業務が多くなりました。社内とは言え、他人が開発をしたシステムを使っているので、いざとなれば自力である程度検証できるようにしておかなければなりません。

ということで、簡単なデータを使って R で主成分分析をしてみます。

iris データの解析

R に組み込まれているデータセット iris を使って主成分分析をしてみます。sepal は花びらの「がく」、petal は「花弁」のことで、それぞれ長さと幅を測定したデータが並んでいます。また species は「種」で、アヤメ属 (iris) の花 setosa, versicolor, virginica の三種類あります。

> head(iris)
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa
> unique(iris$Species)
[1] setosa     versicolor virginica 
Levels: setosa versicolor virginica
> 

そこで、この三種類の花がどのような特徴を持っているか調べてみることにします。「種」を除けば変数が4つしかないので「がく」と「花弁」の長さと幅をボックスプロットで比較する方法と、同じ方法で主成分分析を使うとどのようになるかを比較してみます。

R スクリプト

本記事で使用したスクリプトを以下に示しました。このスクリプトでは ggplot2 と qcc の二つのパッケージを利用しています。お使いの環境にインストールされていなければ、install.packages('ggplot2') のようにしてインストールしてください。

library(ggplot2)
library(qcc)
# GLOBAL
theme_default <- theme(panel.background = element_blank(),
                       panel.border = element_rect(fill = NA),
                       panel.grid.major = element_blank(),
                       panel.grid.minor = element_blank(),
                       strip.background = element_blank(),
                       axis.text.x = element_text(colour = "black"),
                       axis.text.y = element_text(colour = "black"),
                       axis.ticks = element_line(colour = "black"),
                       plot.margin = unit(c(1, 1, 1, 1), "line"))
# -----------------------------------------------------------------------------
#  myboxplot
# -----------------------------------------------------------------------------
myboxplot <- function(pp, ymin, ymax) {
  pp <- pp + geom_boxplot()
  pp <- pp + scale_y_continuous(limits = c(ymin, ymax))
  pp <- pp + theme_default
  return(pp)
}
# -----------------------------------------------------------------------------
#  BOXPLOT for RAW DATA
# -----------------------------------------------------------------------------
df <- as.data.frame(iris)
#
y.min <- min(df[, 1:4])
y.max <- max(df[, 1:4])
# -----------------------------------------------------------------------------
# Sepal.Length
p <- ggplot(df, aes(x = Species, y = Sepal.Length))
p <- myboxplot(p, y.min, y.max)
plot(p)
# Sepal.Width
p <- ggplot(df, aes(x = Species, y = Sepal.Width))
p <- myboxplot(p, y.min, y.max)
plot(p)
# Petal.Length
p <- ggplot(df, aes(x = Species, y = Petal.Length))
p <- myboxplot(p, y.min, y.max)
plot(p)
# Petal.Width
p <- ggplot(df, aes(x = Species, y = Petal.Width))
p <- myboxplot(p, y.min, y.max)
plot(p)
# -----------------------------------------------------------------------------
#  BOXPLOT for PCA DATA
# -----------------------------------------------------------------------------
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)), "_"), "[[", 1)
#
y.min <- min(df_out[, 1:4])
y.max <- max(df_out[, 1:4])
# -----------------------------------------------------------------------------
# PC1
p <- ggplot(df_out, aes(x = Species, y = PC1))
p <- myboxplot(p, y.min, y.max)
plot(p)
# PC2
p <- ggplot(df_out, aes(x = Species, y = PC2))
p <- myboxplot(p, y.min, y.max)
plot(p)
# PC3
p <- ggplot(df_out, aes(x = Species, y = PC3))
p <- myboxplot(p, y.min, y.max)
plot(p)
# PC4
p <- ggplot(df_out, aes(x = Species, y = PC4))
p <- myboxplot(p, y.min, y.max)
plot(p)
# Pareto Chart
PCA <- df_pca$sdev
names(PCA) <- c("PC1", "PC2", "PC3", "PC4")
pareto <- pareto.chart(PCA, ylab = "Variation")
plot(pareto)

生データのボックスプロット

以下に、「がく」と「花弁」の長さと幅について、「種」別に描画したボックスプロットを示しました。議論が分かれるところですが、ここでは全体を比較するため縦軸のスケールを揃えてあります。

これらのボックスプロットをじ〜っと見較べて判る差異を以下にまとめました。

  • 「がく」と「花弁」の大きさは、概ね setosa, versicolor, virsinica の順に大きい。
  • 「花弁 (petal)」の長さの差異が、他の数値に比べて際立っている。
  • 「がく (sepal)」の幅の差異が、他の差異に比べて小さい。

主成分分析後のボックスプロット

「がく」と「花弁」の長さと幅について、prcomp で主成分分析をして、主成分を「種」別に描画したボックスプロットを示しました。生データのボックスプロットの時と同様に、全体を比較するため縦軸のスケールを揃えてあります。

主成分のばらつき(変動)の大きさをパレート図にしたものを下記に示しました。

パレート図の内訳をみると、主成分 PC1 と PC2 で全体の変動の 80% 以上を説明できている事がわかります。

> pareto
     
Pareto chart analysis for PCA
        Frequency   Cum.Freq.  Percentage Cum.Percent.
  PC1   2.0562689   2.0562689  68.9345126   68.9345126
  PC2   0.4926162   2.5488851  16.5145035   85.4490161
  PC3   0.2796596   2.8285447   9.3753300   94.8243460
  PC4   0.1543862   2.9829309   5.1756540  100.0000000
> 

参考サイト

  1. irisの正体 (R Advent Calendar 2012 6日目) - どんな鳥も
  2. R plot PCA using ggplot2 - Boqiang Hu

 

ブログランキング・にほんブログ村へ
にほんブログ村