R システムを用いた主成分分析

このページでは，R システムを用いた主成分分析のための操作手順の一例を図解で説明するここでは，R の princomp() 関数を使う．

【補足】princomp() 関数では，特に指定しない限り，不偏分散共分散行列あるいは不偏相関係数行列が使われます（このことは R の princomp() 関数のマニュアルに書いてあります）．

あらかじめ決めておく事項

主成分分析を行う CSV ファイル名と属性名を決めておいてください．

このページでは，次のように書く．

散布図を作成する CSV ファイル名: C:\R\Book1.csv
属性名：seq, date, USD, EUR, AUD
属性名は英語になっていること．CSV ファイルの第一行目に書いていること（そうなっているものとして，手順を書く）．

前準備

R システムのインストール

【サイト内の関連ページ】

Windows での R のプログラム例 のインストール手順は，別ページ »で説明
なお，Windows での Microsoft R Open のインストール手順は，別ページ »で説明なお，Windows での RStudio のインストール手順は，別ページ »で説明
Ubuntu での R のプログラム例 のインストール手順は，別ページ »で説明
なお，Ubuntu での RStudio のインストール手順は，別ページ »で説明

princomp() 関数を使って主成分分析を行う手順

princomp() 関数では，

不偏相関係数行列から主成分を求めるか
これは，単位が違う変数が混ざっていて，かつ，各変数が含む誤差の分散が同じであるような場合に向く，とは言えます．
不偏分散共分散行列から主成分を求めるか，
これは，単位が違う変数が混ざっていない場合に向くとは言えます．

を選ぶことができる．どちらが「良い」とは簡単に言えません．このページでは，どちらか片方を選ぶことはせず，両方の手順を併記します．

CSVファイルを読み込み，テーブルに格納

(前準備) 使用する CSV ファイルの作成
* ここでは Book1.csv をダウンロードし，分かりやすいディレクトリに置く
（参考：「外国為替データ（時系列データ）の情報源の紹介」の Web ページ）
以下の説明では、
- Windows の場合：　データファイル名： C:\R\Book1.csv
- Linuxの場合：　データファイル名： /tmp/Book1.csv
として説明を続ける．

* 自前の CSV ファイルを使うときの注意： read.table() 関数を使うので，属性名は英語になっていること．属性名は，CSV ファイルの第一行目に書いていること．
使用する CSV ファイルの確認
属性名が CSV ファイルの1行目に書かれていることを確認する．
R の起動
次のコマンドを実行．
◆ 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 のオプション
- X <- ・・・変数 X に読み込むという意味
- C:/R/Book1.csv, "/tmp/Book1.csv" ・・・読み込む CSV ファイル名．Windows では区切りには「/」を使うことに注意．
- header="TRUE" または header="FALSE" ・・・列ラベルが設定されているか
- seq="," や seq="\" や seq=" " や・・・列を区切る記号（CSV ファイルのときは「seq=","」）
- na.string="NA" ・・・ Not a Number には "NA" を使うという意味
- dec="." 　・・・　ファイルで使われている小数点記号（既定値は，ピリオド）
- strip.white=TRUE ・・・個々のデータの先頭や末尾にある「空白文字」を取り除いて読み込む
- skip=＜行数＞・・・読み飛ばし行数
- nrow=＜行数＞・・・読み込み行数
- (その他のオプション) dec: ファイルで使われている小数点記号を指定できる
オブジェクト X の確認
次のコマンドを実行．
edit(X);
次のコマンドを実行．
str(X)

主成分分析（不偏相関係数行列から主成分を求める場合）

主成分分析（不偏相関係数行列から主成分を求める場合）
PC <- princomp(~USD+EUR+AUD, cor=TRUE, data=X) PC
（オプション）確認のため PC の構造の表示
固有ベクトルの表示
* "loadings" という名前が付いているので，「主成分負荷量」のことでは無いだろうかと勘違いしやすいので注意してください．princomp() 関数の場合には，固有ベクトルが求まります．
unclass(loadings(PC))
* 結果が3行3列の行列になっているのは間違いではない． princomp() 関数が扱うデータフレームは，列数を変数の数，行数を標本数とみなし，i 行 j 列目の要素は，第 j 変数の第 i 番目の標本という意味になる．その主成分分析の結果である主成分は，変数の数を行数及び列数とする行列をなす．
* 第１列（一番左の列）が第１主成分，第２列が第２主成分，第３列が第３主成分．主成分の数は元の変数の数と一致する．
固有値の表示
PC$sd^2
固有値の総和は３になる．
主成分負荷量
主成分負荷量は，構想係数とも呼ばれます．相関係数行列から主成分を求めた場合，「主成分負荷量＝（固有値）の平方根×固有ベクトル」が成り立ちます．次の式で求めることができる．
t( PC$sd * t( PC$loadings ) )[, drop = FALSE]
主成分分析（不偏分散共分散行列から主成分を求める場合）
princomp()関数では，「cor=FALSE」と指定する．これは，不偏分散共分散行列を使うという意味．
PC <- princomp(~USD+EUR+AUD, cor=FALSE, data=X) PC
- 固有ベクトルの表示
  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]

主成分分析の結果のグラフ表示

相関係数行列の固有ベクトルがなす行列から２列を選ぶ
２次元の散布図を書きたいので、そのために、主成分分析で得られた主成分の第１列と第２列を選ぶ方法を練習しておきます．
* [,c(1, 2)] と書いて，第１列と第２列を選ぶ．
PC <- princomp(~USD+EUR+AUD, cor=TRUE, data=X) unclass(loadings(PC))[,c(1, 2)]
主成分スコアを求める
「unclass(loadings(PC))[,c(1, 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)] )

biplot() 関数を用いた主成分分析結果のバイプロット
R の biplot() 関数を使うと，第一主成分と第二主成分から２次元の散布図を描くことが簡単にできる．
* 上の plot() 関数での結果と同様のグラフになっているので，確かに，biplot() 関数は，第一主成分と第二主成分から２次元の散布図を描いているということが確認できる．
biplot( PC )
screeplot() 関数
screeplot( PC )

（オプション）より基本的な関数を用いた主成分分析

主成分分析を行う手法としては，大きく分けて，相関係数行列による方法と，分散共分散行列による方法が知られています．

princomp() 関数を使う方法だと、中身がブラックボックスです．より理解を深めるために，固有値分解を行うeigen() 関数と，不偏相関係数行列を求める cor() 関数と，不偏分散共分散行列を求める cor() 関数を組み合わせて，主成分を求めてみます．つまり，以下の通りです．

R の cor() 関数・・・不偏相関係数行列
R の cov() 関数・・・不偏分散共分散行列

* R による分散共分散行列の Web ページで，不偏相関係数行列と不偏分散共分散行列を説明している．

データの読み込み

下記のデータを変数 X に読み込む

X <-read.table("C:/R/Book1.csv", header=TRUE, sep=",", na.strings="NA", strip.white=TRUE )
X

不偏相関係数行列の固有値と固有ベクトルを求める
固有値の総和は３ $vectorsの各列が主成分．つまり，１番左の列が第一主成分
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の各列が主成分．つまり，１番左の列が第一主成分
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]