生物統計学

主成分分析に関するRの関数について

dist(変数)

データは数値行列。次のprincompで使う距離行列を計算する。デフォルトはユークリッド距離。

princomp(変数, 因子, 関数)

データにはあらかじめ計算した距離行列を与える必要がある。また、変数の数がケースの数より多い時はエラーになる。一般的に、次のprcompを使う方がいい。

prcomp(formula,data,scale.=TRUE)

formulaは従属変数を書かずに、分析に使う数値型の独立変数のみを書く。（例　~X1+X2+X3+X4） dataはデータフレーム。分析に使う列が連続している場合はdata[,3:10]などとするほうがformulaを使うよりも簡単。計測値の単位が異なるときなど、数値の桁や変動の範囲が変数ごとに大きく異なる場合はscale.=Tで標準化する。scaleのあとのピリオドに注意。デフォルトでは標準化しない。結果はデータフレームとして返され、$xに主成分得点が入っている。

plot(x,y,pch=番号)

２次元の散布図を描く。pchはプロットの文字を指定するか、記号を数字で示したもの。代表的な記号は0: □、1: ○、2: △、15: ■、16: ●、17: ▲など。詳しくはhelp(pch)を参照。pchにケースに応じた文字や数値のベクトルを与えれば、ケースごとに記号を変えることができる。

library(rgl)

library関数はRに組み込みでないパッケージを読み込んで使えるようにする。パッケージはあらかじめダウンロードしておく必要がある。この場合、3Dグラフィックス用のパッケージrglを読み込んでいる。

rgl.spheres(x,y,z,radius=0.1,col=色)

3Dの散布図を描く。x,y,zは各ケースの座標を表す数値ベクトル。rglパッケージの3Dグラフィックスはマウスを使って回転、拡大、縮小ができる。主成分分析で３つの主成分得点を同時に見ることができ、データの解釈に役立つ。

rgl.snapshot("file_name.png")

rgl関数で描かれたグラフをPNG形式で保存する。

主成分分析の実行例

> pet<-read.table("petunia.csv",sep=",",header=T)
> res<-prcomp(pet[,3:25],scale.=T) # pl, spを除く全データを使う。
> summary(res) # 固有値、寄与率、累積寄与率の表示
Importance of components:
                          PC1     PC2     PC3     PC4     PC5     PC6     PC7     PC8     PC9    PC10    PC11    PC12    PC13    PC14
Standard deviation     3.3646 1.49100 1.35508 1.23312 1.07091 0.98562 0.94125 0.75862 0.69786 0.62428 0.58708 0.53921 0.51377 0.49051
Proportion of Variance 0.4922 0.09666 0.07984 0.06611 0.04986 0.04224 0.03852 0.02502 0.02117 0.01694 0.01499 0.01264 0.01148 0.01046
Cumulative Proportion  0.4922 0.58886 0.66869 0.73481 0.78467 0.82691 0.86542 0.89045 0.91162 0.92857 0.94355 0.95619 0.96767 0.97813
                          PC15    PC16   PC17    PC18   PC19    PC20    PC21    PC22    PC23
Standard deviation     0.44034 0.35269 0.2876 0.22273 0.1797 0.11493 0.06514 0.04459 0.02665
Proportion of Variance 0.00843 0.00541 0.0036 0.00216 0.0014 0.00057 0.00018 0.00009 0.00003
Cumulative Proportion  0.98656 0.99197 0.9956 0.99772 0.9991 0.99970 0.99988 0.99997 1.00000
> res$rotation[,1:4] # 因子負荷量を第４主成分まで表示
             PC1          PC2         PC3          PC4
X1  -0.232902265 -0.007737726 -0.11926598  0.084209120
X3  -0.125560614 -0.038512370 -0.43077888 -0.075944713
X4  -0.073251021 -0.228037568 -0.50934208 -0.027266084
X5  -0.283243440  0.115905707  0.01346202 -0.013690070
X6  -0.280968230  0.129511194  0.07787914 -0.016946203
X7  -0.262326416  0.135132295  0.13562767 -0.026904644
X8  -0.279998306  0.128339399  0.03669440 -0.012244306
X9  -0.275328334  0.082197053  0.02629658 -0.044473096
X10 -0.234093251 -0.138156813 -0.15981333  0.217121593
X11 -0.165980897  0.248784151  0.20097755 -0.272031591
X12 -0.138290237 -0.024120199 -0.23978806  0.293020160
X14 -0.175804168 -0.101687140  0.16463744 -0.145212460
X16 -0.197096009 -0.185338117 -0.14846991  0.071400841
X17 -0.206668986 -0.150850302 -0.07571860  0.177249191
X19 -0.251178368  0.213004676  0.14398787 -0.047846778
X21 -0.272395960  0.004539319 -0.05266473  0.023877109
X24 -0.268917541  0.006342497  0.02638788 -0.004010512
X25 -0.259344731 -0.041710520 -0.07050005  0.060615579
X29  0.004998925 -0.041842991 -0.24654798 -0.629364943
X30 -0.053711753 -0.585100094  0.22898313 -0.034536584
X31 -0.061202688 -0.536945309  0.28471550 -0.005553255
X33 -0.193511651 -0.142001416  0.31716958 -0.050094976
X34 -0.086758434 -0.161283144 -0.10831709 -0.558643570

二次元散布図を描く

> my_pch<-c(16,17)
> plot(res$x[,1:2],pch=my_pch[pet$sp])
> abline(h=0,v=0) # x=0, y=0に線を引く

三次元散布図を描き、ファイルに保存する

> library(rgl)
> my_colors<-c("red", "blue")
> rgl.spheres(res$x[,1:3],radius=0.1,col=my_colors[pet$sp]) # 第１、２、４主成分を使う場合などは、座標はres$x[,c(1,2,4)]と指定する。
> rgl.lines(c(-6,6),0,0,col="cyan")
> rgl.lines(0,c(-6,6),0,col="yellow")
> rgl.lines(0,0,c(-6,6),col="white")
> rgl.texts(c(6,0,0), c(0,6,0) c(0,0,6),text=c("PC1", "PC2", "PC3"),col="white")
> rgl.snapshot("3d_scatter.png")

最後の軸を描く部分は関数にしておくと汎用性がある。次のファイルをダウンロードし、「ファイル」メニューの「ソースを読み込む...」 からそのファイルを選択するか、source関数を使って読み込む。
draw_axis.R

使う時は

> source("draw_axis.R") # またはメニューの「ソースを読み込む」から
> draw.axis(2) # -2から2の3本の軸を描き、デフォルトの軸ラベルx, y, zを書く。
> draw.axis(6,c("PC1","PC2","PC3"))
などとする。

# draw_axis.R
draw.axis <- function( x,labels=c("x","y","z"),colors=c("cyan","yellow","white","white") ) {
  rgl.texts(c(x,0,0),c(0,x,0),c(0,0,x),text=labels,col=colors[4])
  rgl.lines(c(-x,x),0,0,col=colors[1])
  rgl.lines(0,c(-x,x),0,col=colors[2])
  rgl.lines(0,0,c(-x,x),col=colors[3])
}