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()
