このページでは,R システムを用いた主成分分析のための操作手順の一例を図解で説明する ここでは,R の princomp() 関数を使う.
【補足】princomp() 関数では,特に指定しない限り, 不偏分散共分散行列 あるいは 不偏相関係数行列 が使われます (このことは R の princomp() 関数のマニュアルに書いてあります).
主成分分析を行う CSV ファイル名と属性名を決めておいてください.
このページでは,次のように書く.
属性名は英語になっていること.CSV ファイルの第一行目に書いていること(そうなっているものとして,手順を書く).
princomp() 関数では,
これは,単位が違う変数が混ざっていて,かつ,各変数が含む誤差の分散が同じであるような場合に向く,とは言えます.
これは,単位が違う変数が混ざっていない場合に向くとは言えます.
を選ぶことができる.どちらが「良い」とは簡単に言えません.このページでは,どちらか片方を選ぶことはせず, 両方の手順を併記します.
※ ここでは Book1.csv をダウンロードし,分かりやすいディレクトリに置く
(参考: 「外国為替データ(時系列データ)の情報源の紹介」の Web ページ)
以下の説明では、
として説明を続ける.
※ 自前の CSV ファイルを使うときの注意: read.table() 関数を使うので, 属性名は英語になっていること.属性名は,CSV ファイルの第一行目に書いていること.
属性名が CSV ファイルの1行目に書かれていることを確認する.
次のコマンドを実行.
◆ Windows での動作手順例
X <- read.table("C:/R/Book1.csv", header=TRUE, sep=",", na.strings="NA", dec=".", strip.white=TRUE);
X <- read.table("/tmp/Book1.csv", header=TRUE, sep=",", na.strings="NA", dec=".", strip.white=TRUE);
R の read.table のオプション
次のコマンドを実行.
edit(X);
次のコマンドを実行.
str(X)
PC <- princomp(~USD+EUR+AUD, cor=TRUE, data=X) PC
※ "loadings" という名前が付いているので,「主成分負荷量」のことでは無いだろうかと勘違いしやすいので注意してください.princomp() 関数の場合には,固有ベクトルが求まります.
unclass(loadings(PC))
※ 結果が3行3列の行列になっているのは間違いではない. princomp() 関数が扱うデータフレームは, 列数を変数の数,行数を標本数とみなし,i 行 j 列目の要素は,第 j 変数の第 i 番目の標本という意味になる.その主成分分析の結果である主成分は,変数の数を行数及び列数とする行列をなす.
※ 第1列(一番左の列)が第1主成分,第2列が第2主成分, 第3列が第3主成分.主成分の数は元の変数の数と一致する.
PC$sd^2
固有値の総和は3になる.
主成分負荷量は,構想係数とも呼ばれます. 相関係数行列から主成分を求めた場合,「主成分負荷量=(固有値)の平方根×固有ベクトル」が成り立ちます.次の式で求めることができる.
t( PC$sd * t( PC$loadings ) )[, drop = FALSE]
princomp()関数では,「cor=FALSE」と指定する.これは,不偏分散共分散行列を使うという意味.
PC <- princomp(~USD+EUR+AUD, cor=FALSE, data=X) PC
unclass(loadings(PC))
※ 第1列(一番左の列)が第1主成分,第2列が第2主成分, 第3列が第3主成分.主成分の数は元の変数の数と一致する.
PC$sd^2
固有値の総和は3にはならない 不偏分散共分散行列の対角成分の和になる.
分散共分散行列から主成分を求めた場合,「主成分負荷量=(固有値)の平方根×固有ベクトル/(変数の分散)の平方根」が成り立ちます.次の式で求めることができる.
t( PC$sd / sqrt( c( sd( X[,3] ), sd( X[,4] ), sd( X[,5] ) ) ) * t( PC$loadings ) )[, drop = FALSE]
2次元の散布図を書きたいので、そのために、主成分分析で得られた主成分の第1列と第2列を選ぶ方法を練習しておきます.
※ [,c(1, 2)] と書いて,第1列と第2列を選ぶ.
PC <- princomp(~USD+EUR+AUD, cor=TRUE, data=X) unclass(loadings(PC))[,c(1, 2)]
「unclass(loadings(PC))[,c(1, 2)]」で2列を選んでいるので,スコアは2個になる.
※ 「data.matrix( norm[,3:5] ) %*% eigen( cor( norm[,3:5] ) )$vectors[,c(1,2)]」と書いても同じ結果になる.
D <- data.matrix( X[,c(3,5,11)] ) %*% unclass(loadings(PC))[,c(1, 2)] edit(D);
plot( data.matrix( X[,c(3,5,11)] ) %*% unclass(loadings(PC))[,c(1, 2)] )
R の biplot() 関数を使うと, 第一主成分と第二主成分から2次元の散布図を描くことが簡単にできる.
※ 上の plot() 関数での結果と同様のグラフになっているので,確かに,biplot() 関数は, 第一主成分と第二主成分から2次元の散布図を描いているということが確認できる.
biplot( PC )
screeplot( PC )
主成分分析を行う手法としては,大きく分けて,相関係数行列による方法と,分散共分散行列による方法が知られています.
princomp() 関数を使う方法だと、中身がブラックボックスです. より理解を深めるために,固有値分解を行うeigen() 関数と, 不偏相関係数行列を求める cor() 関数と, 不偏分散共分散行列を求める cor() 関数を組み合わせて,主成分を求めてみます. つまり,以下の通りです.
※ R による分散共分散行列 の Web ページで,不偏相関係数行列と 不偏分散共分散行列を説明している.
下記のデータを変数 X に読み込む
X <-read.table("C:/R/Book1.csv", header=TRUE, sep=",", na.strings="NA", strip.white=TRUE ) X
固有値の総和は3 $vectorsの各列が主成分.つまり,1番左の列が第一主成分
eigen( cor( X[,c(3,5,11)] ) )
※ princomp() 関数と比べてみましょう.
同じ結果が得られていることが分かります.
「unclass(loadings(PC))」のところは固有ベクトル,
「PC$sd^2」のところは固有値,
「t( sqrt ( PC$sd ...」のところは主成分負荷量です.
PC <- princomp(~USD+EUR+AUD, cor=TRUE, data=X) unclass(loadings(PC)) PC$sd^2 t( PC$sd * t( PC$loadings ) )[, drop = FALSE]
$vectorsの各列が主成分.つまり,1番左の列が第一主成分
eigen( cov( X[,c(3,5,11)] ) )
※ princomp() 関数と比べてみましょう.
固有ベクトルは同じ値になっていますが,固有値は値が少し違うということが見て取れます.
「unclass(loadings(PC))」のところは固有ベクトル,
「PC$sd^2」のところは固有値,
「t( sqrt ( PC$sd ...」のところは主成分負荷量です.
PC <- princomp(~USD+EUR+AUD, cor=FALSE, data=X) unclass(loadings(PC)) PC$sd^2 t( PC$sd / sqrt( c( sd( X[,3] ), sd( X[,4] ), sd( X[,5] ) ) ) * t( PC$loadings ) )[, drop = FALSE]
本サイトは金子邦彦研究室のWebページである.
資料等の公開では,原則,「クリエイティブコモンズ BY NC SA」として公開するようにしている. PDFファイル,パワーポイントファイルなどには, 「クリエイティブコモンズ BY NC SA」を明記するとともに,ロゴを記載するようにしている(作業が間に合っていない分もあるのでご容赦ください).
公開している資料をご利用になる場合の,再配布の条件,剽窃の防止などについて,別ページ »で説明 再配布や資料改変の際には,そのページをご確認ください.
サイトマップは,サイトマップのページをご覧下さい. 本サイト内の検索は,サイト内検索のページをご利用下さい.
問い合わせ先: 金子邦彦(かねこ くにひこ)