Rによる主成分分析

質問

統計解析用ソフトRを用いて主成分分析を行いたいと考えております。
しかし中々うまくいきません。

そこで、詳しい方がいたっしゃいましたら、主成分分析、主成分負荷量、バイプロットまでのプログラムを教えていただけないでしょうか。

解析対象のデータは以下の様なものです。データ数が２１個、変数が１３個です。

　 F G　H ・・・　　R
1　 6 6　　5 ・・・　　2
2　 5 6　　1 ・・・　　6
3　 6 4　　6 ・・・　　5
： ： ：　　： ・・・　　：
21　3 6　　4 ・・・　　6

よろしくお願いいたします。

回答No.1

kamiyasiro · ベストアンサー

データをスクリプト中に入れましたが、
csvファイルから読み込むことも可能です。
それくらいは、ご存知ですよね。

Rのバイ・プロットは、記号表示となり、
私は好きではありませんが、
下記のようにすれば、できます。

コピペして動作を確認してみて下さい。

# テストデータによる主成分分析
dat <- data.frame(
matrix(c(
3,4,1,4,1,4,1,1,2,2,4,2,2,
3,3,3,3,3,3,3,1,3,3,3,3,3,
1,4,5,3,1,1,1,1,1,1,1,1,1,
3,1,3,1,1,2,2,2,4,4,3,4,2,
3,3,3,4,3,3,3,3,3,3,3,3,3,
3,3,3,3,3,3,3,3,3,3,3,3,3,
4,4,4,4,4,4,3,3,2,2,2,2,2,
4,3,4,3,2,2,2,2,3,3,2,2,2,
3,4,4,4,1,3,3,2,2,2,2,3,3,
3,3,3,3,3,3,4,3,2,2,2,3,3,
3,3,2,4,4,3,2,2,3,2,2,2,2,
3,2,3,3,3,2,3,3,4,3,3,2,3,
4,4,4,4,4,3,3,3,2,2,2,2,2,
3,4,4,4,2,3,3,3,3,3,2,2,3,
3,3,3,3,3,3,3,3,3,3,3,3,3,
3,4,4,3,3,3,3,3,3,3,3,3,3,
3,4,4,4,4,4,3,3,3,3,3,3,4,
3,3,3,3,3,3,3,3,3,3,3,3,3,
3,3,4,3,3,3,2,1,2,2,2,2,2,
3,5,1,5,1,5,1,1,1,1,1,1,1,
3,3,2,2,2,2,2,2,2,2,4,2,4)
,ncol=13,nrow=21,byrow=T)
)
names(dat)=c("F","G","H","I","J","K","L","M","N","O","P","Q","R")
dat
result1 <- prcomp(dat,scale=TRUE)
biplot(result1)

お礼

ご回答いただきありがとうございます。

動作確認いたしましたところ、無事にバイプロットを出力することができました。

大変図々しいこととは承知しておりますが２点質問させてください。
１．このプログラムは相関行列を用いた主成分分析なのでしょうか。私の説明が不十分だったのですが、用いるデータファイルは最小値が１、最大値が６であり、単位の揃ったものです。（正確には単位はありません。）なので今回は分散共分散行列を用いるのが好ましいのではないかと考えております。

２．バイプロットと一緒に、各主成分に対する主成分負荷量の図も出力したいのですが、コードを教えていただませんでしょうか。

専攻は建築系なもので、勉強不足な点が多々ございますが、
何卒よろしくお願いいたします。

回答No.4

kamiyasiro

#1です。

データが１～６の非負の離散値ということが少し気になったので、
おせっかいですがコメントさせて下さい。

もしかして、アンケート結果を解析しようとしてみえますか。
すると、結果が思わしくないと思います。
なぜなら、人による採点の「甘・辛」（バイアス）や
同内容の設問間の「自己相関」が出て、
主成分分析には適さないデータになっているからです。

バイアスを取るには、行が「被験者」なら、行間を中心化します。
自己相関を取るには、列が「設問」なら、列間を中心化します。
これを「二重中心化」といいます。

スクリプトをのせておきます。
走らせてみると、因子負荷量の矢印が四方八方に広がるのが分かると思います。
第１、第２主成分を取り除いた感じです。
第１主成分のサイズファクター（甘いか辛いか）と
第２主成分のシェイプファクター（設問に対する反応の凹凸）を
取り除き偏りが取れて広がるのです。
「コレスポンデンス分析」の結果に近いと思って下さい。

ただ、相対評価（主観的回答）以外のデータには適用しないで下さい。

# テストデータによる主成分分析
par(mfrow=c(1,1))
dat <- data.frame(
matrix(c(
3,4,1,4,1,4,1,1,2,2,4,2,2,
3,3,3,3,3,3,3,1,3,3,3,3,3,
1,4,5,3,1,1,1,1,1,1,1,1,1,
3,1,3,1,1,2,2,2,4,4,3,4,2,
3,3,3,4,3,3,3,3,3,3,3,3,3,
3,3,3,3,3,3,3,3,3,3,3,3,3,
4,4,4,4,4,4,3,3,2,2,2,2,2,
4,3,4,3,2,2,2,2,3,3,2,2,2,
3,4,4,4,1,3,3,2,2,2,2,3,3,
3,3,3,3,3,3,4,3,2,2,2,3,3,
3,3,2,4,4,3,2,2,3,2,2,2,2,
3,2,3,3,3,2,3,3,4,3,3,2,3,
4,4,4,4,4,3,3,3,2,2,2,2,2,
3,4,4,4,2,3,3,3,3,3,2,2,3,
3,3,3,3,3,3,3,3,3,3,3,3,3,
3,4,4,3,3,3,3,3,3,3,3,3,3,
3,4,4,4,4,4,3,3,3,3,3,3,4,
3,3,3,3,3,3,3,3,3,3,3,3,3,
3,3,4,3,3,3,2,1,2,2,2,2,2,
3,5,1,5,1,5,1,1,1,1,1,1,1,
3,3,2,2,2,2,2,2,2,2,4,2,4)
,ncol=13,nrow=21,byrow=T)
)
names(dat)=c("F","G","H","I","J","K","L","M","N","O","P","Q","R")
dat
xs <- dat
dat <- sweep(dat,2,apply(xs,2,mean)) # 列平均を引く
dat <- sweep(dat,1,apply(xs,1,mean)) # 行平均を引く
dat <- sweep(dat,c(1,2),mean(colMeans(xs)),FUN="+") # 全平均を足す
#
result1 <- prcomp(dat)
biplot(result1)
par(ask=TRUE)
#
result1$rotation
p <- result1$rotation[,1:2]
p1 <- p[order(-abs(p[,1])),1]
p2 <- p[order(-abs(p[,2])),2]
par(mfrow=c(1,2))
barplot(p1,horiz=T,names.arg=rownames(p1),las=1,cex.names=0.8)
barplot(p2,horiz=T,names.arg=rownames(p2),las=1,cex.names=0.8)
#

お礼

ご丁寧にありがとうございます。
分析に用いるデータは、アンケート結果ではなく健康診断の結果の様なものです。
F～Rのアルファベットは、例えば血圧ですとか血糖値といった項目（６が最も状態の良い評価、１が最も状態の悪い評価）であり、１～２１は検体数といった感じです。
主成分分析を行うことで、多くの診断結果の中から最も体の状態を把握することのできる因子を読み取ることができるのではないか、と考えました。

頂いたプログラムを回してみました。
確かに中心から放射状に矢印が伸びていることが分かりました。

私の勉強不足にお付き合いいただき本当にありがとうございます。

回答No.3

kamiyasiro

#1です。

グラフを昇順にしたいときは、データを降順にする必要があります。

p1 <- p[order(abs(p[,1])),1]
のorderで並べ替えをしていますので、
p1 <- p[order(-abs(p[,1])),1]
と「－」を入れてやれば、絶対値の符号が逆転しますので降順になります。

horiz=T
は、棒グラフを縦棒グラフから横棒グラフに変える指示です。

ここまで来ると表示上の趣味の問題で、本来のご質問の趣旨から外れますね。

回答No.2

kamiyasiro

#1です。

１．出発行列を相関係数行列にするか、分散共分散行列にするかは、
result1 <- prcomp(dat,scale=TRUE)
の、scaleで指定します。scale=TRUEだと相関係数行列を使います。
デフォルトはscale=FALSEなので、分散共分散から出発したければ、
result1 <- prcomp(dat)
として下さい。

２．主成分負荷量（一般的には因子負荷量と言いますが）は、
バイ・プロットでは、赤い矢印で示されています。
これを、棒グラフにでもしたいのですか？ならば、
result1$rotation
に値が入っていますので、
それをソートしてグラフにすればいいです。

スクリプトを追加しました。
グラフに一時停止を入れてありますので、次に進むには、
プロット窓のどこかをクリックして下さい。

# テストデータによる主成分分析
par(mfrow=c(1,1))
dat <- data.frame(
matrix(c(
3,4,1,4,1,4,1,1,2,2,4,2,2,
3,3,3,3,3,3,3,1,3,3,3,3,3,
1,4,5,3,1,1,1,1,1,1,1,1,1,
3,1,3,1,1,2,2,2,4,4,3,4,2,
3,3,3,4,3,3,3,3,3,3,3,3,3,
3,3,3,3,3,3,3,3,3,3,3,3,3,
4,4,4,4,4,4,3,3,2,2,2,2,2,
4,3,4,3,2,2,2,2,3,3,2,2,2,
3,4,4,4,1,3,3,2,2,2,2,3,3,
3,3,3,3,3,3,4,3,2,2,2,3,3,
3,3,2,4,4,3,2,2,3,2,2,2,2,
3,2,3,3,3,2,3,3,4,3,3,2,3,
4,4,4,4,4,3,3,3,2,2,2,2,2,
3,4,4,4,2,3,3,3,3,3,2,2,3,
3,3,3,3,3,3,3,3,3,3,3,3,3,
3,4,4,3,3,3,3,3,3,3,3,3,3,
3,4,4,4,4,4,3,3,3,3,3,3,4,
3,3,3,3,3,3,3,3,3,3,3,3,3,
3,3,4,3,3,3,2,1,2,2,2,2,2,
3,5,1,5,1,5,1,1,1,1,1,1,1,
3,3,2,2,2,2,2,2,2,2,4,2,4)
,ncol=13,nrow=21,byrow=T)
)
names(dat)=c("F","G","H","I","J","K","L","M","N","O","P","Q","R")
dat
result1 <- prcomp(dat)
biplot(result1)
par(ask=TRUE)
#
result1$rotation
p <- result1$rotation[,1:2]
p1 <- p[order(abs(p[,1])),1]
p2 <- p[order(abs(p[,2])),2]
par(mfrow=c(1,2))
barplot(p1,horiz=T,names.arg=rownames(p1),las=1,cex.names=0.8)
barplot(p2,horiz=T,names.arg=rownames(p2),las=1,cex.names=0.8)
#

お礼

ご対応いただきありがとうございます。

1.そういうことだったのですね。
分散共分散行列なのか相関行列なのか、解決しました。

2.早速プログラムを回したところ、因子負荷量のグラフを出力することができました。
細かいお心遣いもありがとうございます。
まさに棒グラフで示したかったのです。
試しに第５主成分の因子負荷量まで求めてみましたが、問題なく出力することができました。
少し気になる点が１つあります。
「horiz=T」の部分が図中の軸を決めるものだと思うのですが、負荷量の少ないものから表記することは可能なのでしょうか。