今回の勉強会
こないだの日曜日は友人とやっている週次の統計勉強会で、金明哲『Rによるデータサイエンス』の主成分分析の章と因子分析の章を扱いました。
Rによるデータサイエンス - データ解析の基礎から最新手法まで
- 作者: 金明哲
- 出版社/メーカー: 森北出版
- 発売日: 2007/10/13
- メディア: 単行本（ソフトカバー）
- 購入: 36人 クリック: 694回
- この商品を含むブログ (59件) を見る
多変量解析の教科書で理論的な学習は終えた上で、演習としてRでの分析を実行してみるという位置づけです。主成分分析と因子分析についての理論的な説明も、一応この教科書に少し書いてありますが、基本的にはRの演習がメインの教科書なので、考え方や数式の展開などの解説は別書で読むべきでしょう。（我々の勉強会では、永田・棟近『多変量解析入門』を使用しました。過去のエントリ参照。）
なお、次回の勉強会では永田・棟近『多変量解析入門』のクラスター分析を扱う予定。
主成分分析
この『Rによるデータサイエンス』の主成分分析の演習のデータは面白くて、↓の図のとおりなんですが、1から16がサンプル（ケース）で、それぞれのA〜Eとの距離を変数としたデータを分析にかけます。
いわば、A〜Eに対応する５本の座標軸で位置関係が表現されている1〜16のサンプルについて（この図上でA〜Eを何か座標軸の矢印のように見なせるわけではないですけど）、A〜Eの座標軸を主成分分析で２座標（第１主成分と第２主成分）に縮約した上で、第１主成分、第２主成分とそれぞれに対する主成分得点から、もとの位置関係どれぐらい再現されるかを見てみるという感じですね。
やることは簡単で、主成分分析の関数にデータを入れて終了です……。
まず教科書のデータを入力します。めんどくさいですが、ScanSnapで教科書を自炊してOCRをかけていたので、ほぼコピペで再現できましたｗ
# データの入力 temp <- c(50,57,74,94,112,128,140,147,150,147,140,128,112,94,74,57,57,50,57,74,94,112,128,140,147,150,147,140,128,112,94,74,74,57,50,57,74,94,112,128,140,147,150,147,140,128,112,94,94,74,57,50,57,74,94,112,128,140,147,150,147,140,128,112,112,94,74,57,50,57,74,94,112,128,140,147,150,147,140,128) # 分析用の行列に変形 okamoto <- matrix(temp, 16, 5, byrow=F) colnames(okamoto) <- c("A", "B", "C","D", "E")
これを、主成分分析にかけて結果を表示させます。
> oka.pc <- princomp(okamoto) # 主成分分析を実行 > > summary(oka.pc) Importance of components: Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Standard deviation 65.8425 39.5194 3.714616 1.2617966 0.29607143 Proportion of Variance 0.7332 0.2641 0.002334 0.0002693 0.00001483 Cumulative Proportion 0.7332 0.9974 0.999716 0.9999852 1.00000000 > > print(oka.pc) Call: princomp(x = okamoto) Standard deviations: Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 65.8425 39.5194 3.7146 1.2618 0.2961 5 variables and 16 observations.
表示される標準偏差は、主成分得点の標準偏差であり、これは固有値の正の平方根に一致します。
さてこの結果から、第2主成分までで累積寄与率が99.7%になっているので、第1主成分と第2主成分だけでデータを表現すればだいたいOKということがわかります。
第1主成分と第2主成分をヨコ・タテの軸にもつ座標軸上に、元の変数A〜Eの負荷量をプロットすると、なんとなく、元のデータにおけるA〜Eの位置関係が再現されたような感じになる。
plot(oka.pc$loadings[,1:2], type="n”) # 作図デバイスを開き、点のプロットは”n”で非表示にする text(oka.pc$loadings, colnames(okamoto)) # 空白の作図デバイス上に、点ではなく対応するテキストをプロットする
さて、主成分得点は$scoresの中に入っています。
# 因子得点を取り出し > oka.pc$scores Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 [1,] 65.349 -51.01938458029395917 -2.880 -1.8768050106654071 -0.13974 [2,] 91.202 -31.09924073244312837 3.874 -2.3102436809885090 -0.07364 [3,] 100.752 0.00000000000005684 7.365 -0.0000000000002469 0.25640 [4,] 91.202 31.09924073244322784 3.874 2.3102436809883180 -0.07364 [5,] 65.349 51.01938458029403023 -2.880 1.8768050106655920 -0.13974 [6,] 29.709 56.94018795323957960 -5.485 -0.5612449106682362 -0.25085 [7,] -7.558 52.48044511377999299 -3.722 -1.3427139157409904 0.52232 [8,] -40.168 42.38021509369466600 -1.114 -0.6460407202128757 0.02797 [9,] -65.140 30.02037386793610807 1.198 -1.1301118658509794 -0.21794 [10,] -80.742 15.18008481607138371 2.724 -0.2549999546623445 0.29651 [11,] -86.053 -0.00000000000003553 3.444 0.0000000000000373 -0.58567 [12,] -80.742 -15.18008481607146720 2.724 0.2549999546621162 0.29651 [13,] -65.140 -30.02037386793616491 1.198 1.1301118658510223 -0.21794 [14,] -40.168 -42.38021509369470863 -1.114 0.6460407202129304 0.02797 [15,] -7.558 -52.48044511378001431 -3.722 1.3427139157409851 0.52232 [16,] 29.709 -56.94018795323955828 -5.485 0.5612449106685977 -0.25085
なぜか、options(digits=4)などと表示桁数を指定していても、第2主成分と第4主成分だけ、すごい桁数で表示される。。。原因不明。
ちなみに主成分得点の値は（負荷量もだが）、計算のアルゴリズムによって、符号の正負が異なる箇所が出てくる場合がある。教科書では、
主成分分析で求まる主成分および主成分得点の正負の符号は、固有値および固有ベクトルを求めるアルゴリズムに依存する。異なるアルゴリズムの間には正・負の符号が逆になる場合があるため、主成分および主成分得点の散布図の上下が逆になったり、左右が逆になったりすることがある。主成分分析では、変数間、個体間の絶対的関係ではなく、相対的関係を分析するので問題がない。
と解説されている。
こんどは、主成分得点を、第1主成分と第2主成分の軸にプロットしてみます。
# まずplotで作図デバイスを開きつつ、点のプロットは”n”で非表示にしておく。 # いわば良い感じの軸だけを表示するということ。 # その上で、空白の作図デバイス上に、点ではなく対応するテキストをプロットする。 plot(oka.pc$scores[,1], oka.pc$scores[,2], type="n”) text(oka.pc$scores[,1], oka.pc$scores[,2], 1:16)
すこしゆがんでますが、だいたい、丸い形に並びました！
次に、別の面から主成分分析の性能を確かめるために、16個あるサンプルのうち10個だけを使って分析を行い、そこで出てきた負荷量から残りの6個の因子得点の計算を行い、プロットします。
oka.pc1 <- princomp(okamoto[1:10,]) # 1〜10で主成分分析 # さっきと同じように結果をプロット plot(oka.pc1$scores[,1], oka.pc1$scores[,2], type="n”) text(oka.pc1$scores[,1], oka.pc1$scores[,2], 1:10)
oka.pre <- predict(oka.pc1, okamoto[11:16,]) # 上記で1〜10のサンプルによりパラメータ推定されたモデルに、11〜16のデータを代入 # さっきと同じように結果をプロット plot(oka.pre[,1], oka.pre[,2], type="n”) text(oka.pre[,1], oka.pre[,2], 11:16)
predict()関数は、モデルとパラメータを決めた（推定した）後に、そのモデルに説明変数の値を新たに代入して目的変数を計算し、予測を行うためのものですね。回帰分析の結果に新しい値を代入して予測するときにも使った気がする。
さて次に、バイプロットを表示してみます。
biplot(oka.pc)
Biplotとはなんなのかというと、Wikipediaの解説によれば、二次元の座標軸上に、変数と値をいっしょに表示するグラフのこと。
この図をみてると、主成分分析というのが何をやっている分析なのかがよく分かりますね。
赤い矢印で表されているA〜Eは、もとの変数（標準化されてる）の主成分負荷量を表しており、左と下の座標軸をみればいい。つまりたとえばCというのは第1主成分とほぼ並行な変数なので、第1主成分に対する負荷量はほぼゼロとなっている。さっきの「loadings」の数字をみると実際そうなっています。
1〜16の点は主成分得点でプロットされており、上と右の座標軸をみればいい。さっきのscoresの値をみると、実際そうなっています。
続けて、Rの練習用データirisの４つの変数を２つに圧縮する主成分分析も行う。
ちなみにRを使ってる人にはおなじみですが、irisというのは、フィッシャーがあやめの花の研究で使ったデータらしい。（参考）
最初のほうの中身を表示してみる。
> head(iris) Sepal.Length Sepal.Width Petal.Length Petal.Width Species 1 5.1 3.5 1.4 0.2 setosa 2 4.9 3.0 1.4 0.2 setosa 3 4.7 3.2 1.3 0.2 setosa 4 4.6 3.1 1.5 0.2 setosa 5 5.0 3.6 1.4 0.2 setosa 6 5.4 3.9 1.7 0.4 setosa
１〜４列目の変数を２変数に縮約して、主成分得点をプロットしますが、元データの５列目に入っているsetosa、versicolor、virginicaという文字列を、１〜３の数字に変えてラベルとして用います。
iris.pca <- princomp(iris[,-5]) plot(iris.pca$scores, type="n") lab <- as.numeric(iris[,5]) text(iris.pca$scores, label=lab, col=lab)
バイプロットも出してみる。
biplot(iris.pca)
ごちゃごちゃして何のことかわかりませんね。
次に、ade4というパッケージも試してみます。
> install.packages("ade4") > library(ade4) > iris.dudi <- dudi.pca(iris[,1:4]) Select the number of axes: 2 # 軸の数を入力 > s.arrow(iris.dudi$c1, lab=colnames(iris[1:4]))
これは、主成分負荷量を軸にして元の変数の位置関係を表現している。
> s.class(iris.dudi$li, iris[,5], cpoint=1)
５列目の変数でまとめた上で、主成分得点の分布を見ている図ですね。何が中心になってるのかは知らない。
因子分析
次は因子分析です。これも教科書の演習をそのまま実行してみる。
まずデータをつくる。
seiseki <- matrix(c(89,90,67,46,50,57,70,80,85,90,80,90,35,40,50,40,60,50,45,55,78,85,45,55,60,55,65,80,75,85,90,85,88,92,95), 7, 5, byrow = TRUE) colnames(seiseki) <- c("算数", "理科", "国語", "英語", "社会") rownames(seiseki) <- c("田中", "佐藤", "鈴木", "本田", "川端", "吉野", "斉藤")
factanal()関数で分析します。
> seiseki.fac <- factanal(seiseki, factors=2) > print(seiseki.fac) Call: factanal(x = seiseki, factors = 2) Uniquenesses: 算数 理科 国語 英語 社会 0.005 0.029 0.241 0.005 0.006 Loadings: Factor1 Factor2 算数 0.997 理科 -0.188 0.967 国語 0.871 英語 0.997 社会 0.989 -0.128 Factor1 Factor2 SS loadings 2.768 1.946 Proportion Var 0.554 0.389 Cumulative Var 0.554 0.943 Test of the hypothesis that 2 factors are sufficient. The chi square statistic is 1.73 on 1 degree of freedom. The p-value is 0.188 >
教科書には書いてませんが、負荷量が小さい値を非表示にされるのは困るので、
> print(seiseki.fac, cutoff=FALSE) （省略） Loadings: Factor1 Factor2 算数 0.043 0.997 理科 -0.188 0.967 国語 0.871 -0.018 英語 0.997 -0.028 社会 0.989 -0.128 （省略）
とやることが多いと思います。負荷量が大きい順番に並べたいときは、
> print(seiseki.fac, cutoff=FALSE, sort=TRUE) （略） Loadings: Factor1 Factor2 国語 0.871 -0.018 英語 0.997 -0.028 社会 0.989 -0.128 算数 0.043 0.997 理科 -0.188 0.967 （略）
とします。
で、どの因子がどの項目に高い負荷量を持っているかを可視化するために、この教科書では、棒グラフを書けといってます。
par(mfrow=c(1,2)) # 教科書には書いてないけど２つのグラフを横に並べるときはこうする barplot(seiseki.fac$loadings[,1], col="lightblue") #第１因子 barplot(seiseki.fac$loadings[,2], col="lightblue") #第２因子
理系科目に反応している因子と、文系科目に反応している因子であることが分かりますね。
しかし個人的には、結局、
write.csv(seiseki.fac$loadings, "seiseki.csv")
とやって負荷量の行列をCSVに出して、エクセルでセルに色を付けるとかしたほうがみやすい気がします。まぁ、グラフの方が比率は直感的に分かるのでいいんですが、項目が増えてくると、因子負荷行列をエクセルのフィルターで並べ替えたりしながら見るのがいいんじゃないかなと。心理尺度なんかだと、どの項目を分析から削るという判断にも使えるし。
エクセルだと文字化けしますがめんどくさいので文字コードの変換とかしてません。
最後に、バイプロットを出してみる。
par(mfrow=c(1,2), mar = c(1, 5, 1, 4)) # marは単なる余白調整 biplot(seiseki.fac2$scores, seiseki.fac2$loadings) biplot(seiseki.fac3$scores, seiseki.fac3$loadings)
