R言語で主成分分析を実行する関数やその実行例を解説していきます。
多変量解析において各変数の寄与率を調べるためや多変量データの要約に用いられる主成分分析をR言語で実行する流れを詳しくみていきます。
実行例では、データセットを用いた主成分の算出だけでなく、バイプロットなどの主成分分析に関するグラフを作成する手順についてもみていきます。
この記事で扱うRスクリプトは以下よりダウンロードできます。
主成分分析
主成分分析の概要とR言語の関数をまとめました。
この記事で扱う統計解析の概要と関数を紹介します。
主成分分析の概要
主成分分析
観測値\(p\)次元母集団からの\(n\)個の観測値が得られたとき、標本平均が\(\boldsymbol{0}\)になるように座標をずらしたデータを\(n\times p\)行列\(\boldsymbol{X} =(\boldsymbol{x}_1, \ldots, \boldsymbol{x}_{n})^T\)で表す。\(\boldsymbol{X}\)の標本共分散行列を\(\boldsymbol{S} =n^{-1} \boldsymbol{X}^T\boldsymbol{X}\)とするとき、線形結合\(\boldsymbol{x}_i^T\boldsymbol{w}_i,\ i = 1, 2, \ldots, n\)の標本共分散\(\boldsymbol{x}_i^T\boldsymbol{w}_i\)が最大かつ\(i\neq j\)に対し\(\boldsymbol{w}_i\)と\(\boldsymbol{w}_j\)が直交し、ノルムが\(\|\boldsymbol{w}_i\|^2 = 1\)ようなベクトル\(\boldsymbol{w}_i\)を考える。\(\boldsymbol{w}_i\)は次のようにして求めることができる。
ここに\(\hat{\boldsymbol{X}}\)は
直交行列\(\boldsymbol{W} = (\boldsymbol{w}_1 , \ldots, \boldsymbol{w}_p)\)を定義する。\(X\)を\(w_i,\ i = 1, 2, \ldots, p\)で変換したデータ
は主成分得点と呼ばれ、行列\(\boldsymbol{W}\)は回転行列と呼ばれる（標本共分散行列の固有ベクトル）。さらに変換後の標本共分散行列は次で表される。
\(l_1 \geq l_2 \geq \cdots \geq l_p \geq 0\)は標本共分散行列の固有ベクトル\(\boldsymbol{w}_1, \ldots, \boldsymbol{w}_p\)に対応する固有値である。
関数prcomp
次の関数prcompを用います。
prcomp(x, ...)
## S3 method for class 'formula'
prcomp(formula, data = NULL, subset, na.action, ...)
## Default S3 method:
prcomp(x, retx = TRUE, center = TRUE, scale. = FALSE,
tol = NULL, rank. = NULL, ...)
## S3 method for class 'prcomp'
predict(object, newdata, ...)
prcompで算出した固有値は次の特異値分解に基づくものとなっています。
また、関数prcompの引数は以下の通りです。
| formula | 目的変数の表現をもたないformula型のオブジェクト。numeric型の変数から構成される。 |
| data | numeric型のデータフレーム。formulaで用いた変数が含まれるデータフレーム。 |
| subset | dataのうちどの行を用いるか表すベクトル。 |
| na.action | dataに欠測値NAがあるときにどうするか。 |
| ... | 他のメソッドに適用する引数。formulaを用いる際にscaleやtolを指定するときに用いる。 |
| x | numeric型の行列。主成分分析に用いるデータ。 |
| retx | logical型。固有ベクトルを算出するかどうか。 |
| center | logical型。各変数を0が中心になるようにシフトするかどうか。 |
| scale | logical型。主成分分析を実行する前に、各変数を標本分散が1になるように縮尺を変えるかどうか。 |
| tol | numeric型。tolで指定した値以下の成分は除去される。 |
関数princomp
prcompの他にも次の関数princompでも主成分分析を実行することができます。
princomp(x, …)
# S3 method for formula
princomp(formula, data = NULL, subset, na.action, …)
# S3 method for default
princomp(x, cor = FALSE, scores = TRUE, covmat = NULL, subset = rep_len(TRUE, nrow(as.matrix(x))), fix_sign = TRUE, …)
# S3 method for princomp
predict(object, newdata, …)
princompで計算される固有値は次のスペクトル分解に基づいています。
以下、関数princompの引数となります。
| formula | 目的変数の表現をもたないformula型のオブジェクト。numeric型の変数から構成される。 |
| data | numeric型のデータフレーム。formulaで用いた変数が含まれるデータフレーム。 |
| subset | dataのうちどの行を用いるか表すベクトル。 |
| na.action | dataに欠測値NAがあるときにどうするか。 |
| x | numeric型の行列。主成分分析に用いるデータ。 |
| cor | logical型。固有値の計算の際に相関行列を用いるかどうか。 |
| scores | logical型。各主成分得点を計算するかどうか。 |
| covmat | cov.wtが返す共分散行列または共分散のリスト。 |
| fix_sign | logical型。第一主成分の負荷量とスコアが負とならないように選ぶかどうか。 |
| ... | 他のメソッドに適用する引数。formulaを用いる際にscaleやtolを指定するときに用いる。 |
| object | princompから受け継ぐクラスのオブジェクト。 |
| newdata | 予測したいデータが含まれるデータフレームまたは行列。元のオブジェクトがデータフレームまたは行列の列名を含むformulaを用いている場合、同じ名前の列名をもつデータを含まなければならない。 |
実行例
主成分分析の手順をみていきます。
先ほど紹介した関数prcompやprincompを用いた主成分分析の実行方法や特異値分解との関係について解説します。
グラフの作成のためにパッケージggplot2。あらかじめにインストールしておいてください。
1 2 | install.packages("ggplot2") library(ggplot2) |
主成分分析による変数の寄与率の計算
まず、主成分分析を用いた多変量データの変数の寄与率の算出の流れを紹介します。
データセットに次のirisを用います。
関数headを実行すれば分かるように、irisは4つのアヤメの長さに関するデータと種類（Species）から構成されています。
1 2 3 4 5 6 7 8 9 | > dataset <- iris > 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 |
このirisのデータのSepal.Length、Sepal.Width、Petal.Length、Petal.Widthの4つの変数の寄与率を主成分分析で計算していきます。
上で紹介した関数prcompまたはprincompで主成分分析を実行することができます。
次に示すように関数の引数にデータセットを渡すことで簡単に主成分分析が行えます。prcomでは引数scale、princompでは引数corをTRUEにしていますが、これは標本分散が1になるようにデータのスケールを変更するためです。
また、princompでfix_sign=FALSEとしているのは、計算される第一主成分の負荷量の符号をprcompで計算されるものと同じにするためです。
1つのデータセットに対しデータセットの正規化なしで実行しても問題ないですが、2つ以上の群からのデータセットに対し主成分分析を実行し、変数の寄与率を解釈する際には注意が必要です。
1 2 | PCAResult_SVD <- prcomp(subset(dataset, select = -Species), scale = TRUE) #特異値分解 PCAResult_SD <- princomp(subset(dataset, select = -Species), cor = TRUE, fix_sign = FALSE) #スペクトル分解 |
PCAResult_SVDを参照すると次のように主成分分岐の結果（固有値や固有ベクトル）がコンソールに出力されます。
1 2 3 4 5 6 7 8 9 10 | > PCAResult_SVD Standard deviations (1, .., p=4): [1] 1.7083611 0.9560494 0.3830886 0.1439265 Rotation (n x k) = (4 x 4): 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 |
これらの結果を取得するには、次のようにオブジェクトPCAResult_SVDの後ろに各変数の名前を付けます。xで主成分得点、sdevで固有値の平方根、centerで各変数の標本平均、rotationで固有ベクトルを取得することができます。寄与率を取得するには5行目のようにsdevから計算してあげます。固有値の比から寄与率を計算することができます。
1 2 3 4 5 6 7 8 | scores <- PCAResult_SVD$x #主成分得点 sdevs <- PCAResult_SVD$sdev #標準偏差（固有値の平方根） centers <- PCAResult_SVD$center #標本平均 rotations <- PCAResult_SVD$rotation #回転行列（固有ベクトル） contribRatio <- sdevs^2 / sum(sdevs^2) #寄与率 PCAResult_SD$sdev PCAResult_SD$center |
これらの結果をcsvファイルに保存したいときは次のようにデータフレームに保存しwrite.csvでcsvファイルに出力します。下の画像のcsvファイルが作業ディレクトリに作成されます。
高次元のデータを要約したいときはxの主成分得点を用いますが、今用いているデータセットは4次元のirisなので、各変数の寄与率に関する結果のみ出力しました。
1 2 | resultTable <- data.frame(principal.sdev = sdevs, rotation = rotations, contrib = contribRatio) write.csv(resultTable, "主成分分析.csv") |
主成分分析結果
また、関数screeplotで各変数の固有値をプロットすることができます。各変数の寄与率をグラフで表現することが可能です。
prcompのオブジェクトをそのままscreeplotに渡すと、固有値の平方根であるsdevに名前がないためx軸の目盛りが表示されません。そのため1行目と2行目でsdevに名前を与えています。（princompの場合最初から名前が付いているため、名前を付ける必要はありません。）
1 2 3 | names(sdevs) <- paste0("Comp.", seq_len(length(sdevs ))) PCAResult_SVD$sdev <- sdevs screeplot(PCAResult_SVD) |
screeplot
パッケージggplot2でスクリープロットを描きたいときは次のように関数geom_colを用います。
1 2 3 | ggplot(data.frame(variance = PCAResult_SVD$sdev), aes(x = as.factor(names(PCAResult_SVD$sdev)), y = variance)) + geom_col() + labs(x = "") |
screeplot (ggplot2)
主成分分析による次元の削減
続いて主成分分析を用いた高次元データの要約について見ていきます。
主成分分析により寄与率の大きい主成分で次元の削減を行っていきます。
データセットに次のパッケージMASSのBostonを用います。Bostonは次のように住宅価格に関する変数を持つデータとなります。
1 2 3 4 5 6 7 8 9 10 11 | > #次元の削減 > library(MASS) > dataset <- Boston > head(dataset) crim zn indus chas nox rm age dis rad tax ptratio black lstat medv 1 0.00632 18 2.31 0 0.538 6.575 65.2 4.0900 1 296 15.3 396.90 4.98 24.0 2 0.02731 0 7.07 0 0.469 6.421 78.9 4.9671 2 242 17.8 396.90 9.14 21.6 3 0.02729 0 7.07 0 0.469 7.185 61.1 4.9671 2 242 17.8 392.83 4.03 34.7 4 0.03237 0 2.18 0 0.458 6.998 45.8 6.0622 3 222 18.7 394.63 2.94 33.4 5 0.06905 0 2.18 0 0.458 7.147 54.2 6.0622 3 222 18.7 396.90 5.33 36.2 6 0.02985 0 2.18 0 0.458 6.430 58.7 6.0622 3 222 18.7 394.12 5.21 28.7 |
早速、関数prcompを用いてBostonの数値データについて主成分分析を行います。
1 2 3 4 5 6 7 | vars <- c("crim", "zn", "indus", "nox", "rm", "age", "dis", "tax", "ptratio", "black", "lstat", "medv") PCAResult <- prcomp(dataset[, vars], scale = TRUE) scores <- PCAResult$x contribRatio <- PCAResult$sdev^2 / sum(PCAResult$sdev^2) loadings <- as.data.frame(PCAResult$rotation) colnames(loadings) <- paste0("PC", seq(ncol(loadings))) |
次に示すように関数biplotを実行することで、第一主成分と第二主成分の散布図と主成分の負荷量をプロットすることができます。
赤色の矢印で示される負荷量のベクトルと主成分軸から12個の変数がどのように2変数（第一主成分と第二主成分）に集約されているのかを視覚的にとらえることができます。
1 | biplot(PCAResult, cex = c(0.5, 1)) |
biPlot
ggplot2でbiplotをプロットするには次のように実行します。パッケージggbiplotというbiplotのためのパッケージがありますがRのバージョンによってインストールできないことがあるので、通常のggplot2を用いた方法となります。
1 2 3 4 5 6 7 8 9 10 11 12 | xRatio <- (max(PCAResult$x[, "PC1"]) - min(PCAResult$x[, "PC1"])) / (max(loadings[, "PC1"]) - min(loadings[, "PC1"])) yRatio <- (max(PCAResult$x[, "PC2"]) - min(PCAResult$x[, "PC2"])) / (max(loadings[, "PC2"]) - min(loadings[, "PC2"])) ggplot(data.frame(PCAResult$x, rad = as.factor(dataset$rad))) + geom_point(aes(x = PC1, y = PC2, color = rad)) + stat_ellipse(geom = "polygon", aes(x = PC1, y = PC2, fill = rad), alpha = 0.25, level = 0.95) + labs(x = paste0("PC1 (", round(100 * contribRatio[1], 2), "%)"), y = paste0("PC2 (", round(100 * contribRatio[2], 2), "%)")) + geom_segment(data = data.frame(PC1 = xRatio * loadings$PC1, PC2 = yRatio * loadings$PC2), aes(x = 0, y = 0, xend = PC1, yend = PC2), arrow = arrow(length = unit(0.05, "npc")), lwd = 1) + geom_text(data = data.frame(PC1 = xRatio * loadings$PC1, PC2 = yRatio * loadings$PC2), aes(x = PC1, y = PC2 + 0.5, label = vars)) |
biplot (ggplot)
バイプロットを描くことで高次元データを2次元のグラフに落とし込むことができ、データの解釈が容易になります。
グラフの矢印は各12個の変数の固有ベクトルと平行なベクトルであり、第一主成分PC1が大きいデータはage（年齢）、nox（一酸化窒素濃度）、indus（小売業以外の商業施設の割合）、tax（10000ドルあたりの不動産税率の総計）、crim（人工一人当たりの犯罪数）、lstat（低所得者の割合）、ptratio（児童と教師の割合）が大きい傾向にある。
反対に第二主成分PC1が小さいデータはrm（平均部屋数）、medv（住宅価格）、black（黒色人種の割合）、zn（25,000平方フィートを超える面積の住宅地の割合）、dis（雇用センターまでの距離）が小さい傾向にある。第二主成分についても同様のことがいえる。
第一主成分について12個のベクトルが左右に分かれているため、データを第一主成分が正負について2つのグループに分けると、第一主成分が負であるデータは比較的良い地域のものであり、正であるデータは悪い地域であることが確認できる。
またrad（高速道路へのアクセスの指数）別でプロットすることでどのような地域のデータであるのか把握しやすくなっている。radが小さいほど第一主成分は小さくなりdis（雇用センター）への距離が小さく、特にradが8だとmedv（住宅価格）が大きくなる傾向があるの分かる。
特異値分解との関係
最後に特異値分解との関係について簡単にみていきます。
主成分分析の概要で述べたように、主成分分析の結果の値は特異値分解の各行列に対応しています。
次のようにprcompやprincompを使わなくても、関数svdで主成分分析を行うこともできます。
1 2 3 4 5 6 7 8 | #特異値分解との関係 dataset <- iris stdDataset <- scale(subset(dataset, select = -Species)) SVD <- svd(cor(stdDataset)) PCA_loadings <- SVD$v PCA_rotation <- SVD$u PCA_scores <- stdDataset %*% PCA_rotation |
上で計算された値は、 固有ベクトルの符号を除いて一致しています。Sepal.Lengthに対応する固有ベクトルの符号が逆転しているため主成分得点の符号も逆転していますが、それ以外はすべて一致しています。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 | > PCA_rotation [,1] [,2] [,3] [,4] [1,] -0.5210659 -0.37741762 0.7195664 0.2612863 [2,] 0.2693474 -0.92329566 -0.2443818 -0.1235096 [3,] -0.5804131 -0.02449161 -0.1421264 -0.8014492 [4,] -0.5648565 -0.06694199 -0.6342727 0.5235971 > > PCAResult_SVD$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 > > PCA_loadings [,1] [,2] [,3] [,4] [1,] -0.5210659 -0.37741762 0.7195664 0.2612863 [2,] 0.2693474 -0.92329566 -0.2443818 -0.1235096 [3,] -0.5804131 -0.02449161 -0.1421264 -0.8014492 [4,] -0.5648565 -0.06694199 -0.6342727 0.5235971 > > PCAResult_SD$loadings Loadings: Comp.1 Comp.2 Comp.3 Comp.4 Sepal.Length 0.521 -0.377 0.720 0.261 Sepal.Width -0.269 -0.923 -0.244 -0.124 Petal.Length 0.580 -0.142 -0.801 Petal.Width 0.565 -0.634 0.524 Comp.1 Comp.2 Comp.3 Comp.4 SS loadings 1.00 1.00 1.00 1.00 Proportion Var 0.25 0.25 0.25 0.25 Cumulative Var 0.25 0.50 0.75 1.00 > > head(PCA_scores) [,1] [,2] [,3] [,4] 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 > > head(PCAResult_SVD$x) 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 |
まとめ
R言語で主成分析を実行する関数とその実行例について見てきました。
関数prcompやprincompで主成分分析を実行することができます。prcompは特異値分解、princompはスペクトル分解を基に計算しているため、結果に微妙な違いが出ることがあります。
また、主成分分析のグラフであるバイプロットは関数biplot、スクリープロットは関数screeplotでプロットすることができます。
