Rで主成分分析する基本

Rで主成分分析を行う方法を自分用にまとめておきます。
理論は省略。

メイン関数stats::prcomp

  • statsは基本ロードされるので個別にライブラリを入れる必要は無し。
  • 第一引数に分析対象のデータフレーム。自動的に全変数を主成分分析するのでカテゴリ変数を含むとエラー
  • オプションscaleでスケール合わせ

prcomp(iris[-5], scale = TRUE)

注意

主成分分析は変数間のスケールが合っていなければ解に意味を見出しにくいので、普通はスケール合わせを行う必要がある。
オプションscale = TRUEでスケール合わせも同時に処理する。

結果取得

  • 関数prcompの結果prcompクラスのオブジェクトが得られる。
  • $rotation: loading vector(主成分のベクトル)が主成分ごとにまとまった行列(変数×主成分)
  • $x: スコア(観測値の主成分方向への射影値)の各主成分ごとにまとまった行列(観測×主成分)
  • summaryでそれぞれの主成分の重要度を取得することができる。
pr.out <- prcomp(iris[-5], scale = TRUE)

pr.out$rotation
#                    PC1         PC2        PC3        PC4
#Sepal.Length  0.5210659 -0.37741762  0.7195664  0.2612863
#Sepal.Width  -0.2693474 -0.92329566 -0.2443818 -0.1235096
#Petal.Length  0.5804131 -0.02449161 -0.1421264 -0.8014492
#Petal.Width   0.5648565 -0.06694199 -0.6342727  0.5235971

pr.out$x %>% head
#           PC1        PC2         PC3          PC4
#[1,] -2.257141 -0.4784238  0.12727962  0.024087508
#[2,] -2.074013  0.6718827  0.23382552  0.102662845
#[3,] -2.356335  0.3407664 -0.04405390  0.028282305
#[4,] -2.291707  0.5953999 -0.09098530 -0.065735340
#[5,] -2.381863 -0.6446757 -0.01568565 -0.035802870
#[6,] -2.068701 -1.4842053 -0.02687825  0.006586116


#scoreが理屈通り計算されているか
pr.out$x[1,1]
#       PC1
# -2.257141

iris1.scaled <- (iris[1,-5] - pr.out$center) / pr.out$scale
sum(iris1.scaled * pr.out$rotation[,1]) 
# [1] -2.257141

summary(pr.out)
#Importance of components:
#                          PC1    PC2     PC3     PC4
#Standard deviation     1.7084 0.9560 0.38309 0.14393
#Proportion of Variance 0.7296 0.2285 0.03669 0.00518
#Cumulative Proportion  0.7296 0.9581 0.99482 1.00000

注意

主成分分析の結果は符号の任意性があるため、実行するたびに逆転してもびっくりしないこと

可視化

  • 関数biplot一発でbiplotを描くことができる。
    • オプションscale=0もしておく。(scoreベクトルのスケール合わせ)
  • 汎用関数plotでscree plotを描けるけどあんまり美しくない・・・
biplot(pr.out, scale = 0)

#見栄えや解釈のため符号を反転させてもよい
pr.out$rotation <- pr.out$rotation * -1
pr.out$x <- pr.out$x * -1
biplot(pr.out, scale = 0)

#scree plot
plot(pr.out)

#もうちょっとまともなscree plot
library(tibble)
library(dplyr)
library(ggplot2)
summary(pr.out)$importance[1,] %>%
  enframe %>%
  ggplot(aes(name, value)) +
    geom_line(aes(group = 1)) +
    geom_point()

相対的な重要度も下記のようにちょっとコードを変更すれば見れる。

summary(pr.out)$importance[2,] %>%
  enframe %>%
  ggplot(aes(name, value)) +
    geom_line(aes(group = 1)) +
    geom_point()
相対重要度
summary(pr.out)$importance[3,] %>%
  enframe %>%
  ggplot(aes(name, value)) +
    geom_line(aes(group = 1)) +
    geom_point()
累積相対重要度

参考文献

コメントを残す