トップページデータ処理主成分分析,次元削減R システムを用いた主成分分析

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

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

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

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

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

このページでは,次のように書く.

前準備

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

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

princomp() 関数では,

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

CSVファイルを読み込み,データテーブルに格納

  1. (前準備) 使用する CSV ファイルの作成

    ※ ここでは Book1.csv をダウンロードし,分かりやすいディレクトリに置く

    (参考: 「外国為替データ(時系列データ)の情報源の紹介」の Web ページ

    以下の説明では、

    として説明を続ける.

    自前の CSV ファイルを使うときの注意: read.table() 関数を使うので, 属性名は英語になっていること.属性名は,CSV ファイルの第一行目に書いていること.

  2. 使用する CSV ファイルの確認

    属性名が CSV ファイルの1行目に書かれていることを確認する.

  3. 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: ファイルで使われている小数点記号を指定できる
  4. オブジェクト X の確認

    端末で,次のコマンドを実行.

    edit(X);
    

    端末で,次のコマンドを実行.

    str(X)
    

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

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

    PC <- princomp(~USD+EUR+AUD, cor=TRUE, data=X)
    PC
    
  2. (オプション)確認のため PC の構造の表示
  3. 固有ベクトルの表示

    ※ "loadings" という名前が付いているので,「主成分負荷量」のことでは無いだろうかと勘違いしやすいので注意してください.princomp() 関数の場合には,固有ベクトルが求まります.

    unclass(loadings(PC))
    

    結果が3行3列の行列になっているのは間違いではない. princomp() 関数が扱うデータフレームは, 列数を変数の数,行数を標本数とみなし,i 行 j 列目の要素は,第 j 変数の第 i 番目の標本という意味になる.その主成分分析の結果である主成分は,変数の数を行数及び列数とする行列をなす.

    ※ 第1列(一番左の列)が第1主成分,第2列が第2主成分, 第3列が第3主成分.主成分の数は元の変数の数と一致する.

  4. 固有値の表示

    PC$sd^2
    

    固有値の総和は3になる.

  5. 主成分負荷量

    主成分負荷量は,構想係数とも呼ばれます. 相関係数行列から主成分を求めた場合,「主成分負荷量=(固有値)の平方根×固有ベクトル」が成り立ちます.次の式で求めることができる.

    t( PC$sd * t( PC$loadings ) )[, drop = FALSE]
    
  6. 主成分分析不偏分散共分散行列から主成分を求める場合)

    princomp()関数では,「cor=FALSE」と指定する.これは,不偏分散共分散行列を使うという意味.

    PC <- princomp(~USD+EUR+AUD, cor=FALSE, data=X)
    PC
    

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

  1. 相関係数行列の固有ベクトルがなす行列から2列を選ぶ

    2次元の散布図を書きたいので、そのために、主成分分析で得られた主成分の第1列と第2列を選ぶ方法を練習しておきます.

    ※ [,c(1, 2)] と書いて,第1列と第2列を選ぶ.

    PC <- princomp(~USD+EUR+AUD, cor=TRUE, data=X)
    unclass(loadings(PC))[,c(1, 2)]
    
  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); 
    
  3. 第一主成分と第二主成分から2次元の散布図を描く.

    plot( data.matrix( X[,c(3,5,11)] ) %*% unclass(loadings(PC))[,c(1, 2)] )
    
  4. biplot() 関数を用いた主成分分析結果のバイプロット

    R の biplot() 関数を使うと, 第一主成分と第二主成分から2次元の散布図を描くことが簡単にできる.

    ※ 上の plot() 関数での結果と同様のグラフになっているので,確かに,biplot() 関数は, 第一主成分と第二主成分から2次元の散布図を描いているということが確認できる.

    biplot( PC )
    
  5. screeplot() 関数

    screeplot( PC )
    

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

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

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

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

  1. データの読み込み

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

    X <-read.table("C:/R/Book1.csv", header=TRUE, sep=",", na.strings="NA", strip.white=TRUE )
    X
    
  2. 不偏相関係数行列の固有値と固有ベクトルを求める

    固有値の総和は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]
    
  3. 不偏分散共分散行列の固有値と固有ベクトルを求める

    $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]
    

princomp2() 関数を用いた主成分分析

http://aoki2.si.gunma-u.ac.jp/R/princomp2.html で公開されている princomp2() 関数が便利です.実行例を以下に示しておきます.