StatsBeginner: 初学者の統計学習ノート

統計学およびR、Pythonでのプログラミングの勉強の過程をメモっていくノート。たまにMacの話題。

勉強会:主成分分析・因子分析(金明哲『Rによるデータサイエンス』)

今回の勉強会

 こないだの日曜日は友人とやっている週次の統計勉強会で、金明哲『Rによるデータサイエンス』の主成分分析の章と因子分析の章を扱いました。


Rによるデータサイエンス - データ解析の基礎から最新手法まで

Rによるデータサイエンス - データ解析の基礎から最新手法まで


 多変量解析の教科書で理論的な学習は終えた上で、演習としてRでの分析を実行してみるという位置づけです。主成分分析と因子分析についての理論的な説明も、一応この教科書に少し書いてありますが、基本的にはRの演習がメインの教科書なので、考え方や数式の展開などの解説は別書で読むべきでしょう。(我々の勉強会では、永田・棟近『多変量解析入門』を使用しました。過去のエントリ参照。)
 なお、次回の勉強会では永田・棟近『多変量解析入門』のクラスター分析を扱う予定。
 
 

主成分分析

 この『Rによるデータサイエンス』の主成分分析の演習のデータは面白くて、↓の図のとおりなんですが、1から16がサンプル(ケース)で、それぞれのA〜Eとの距離を変数としたデータを分析にかけます。


f:id:midnightseminar:20140828233121p:plain


 いわば、A〜Eに対応する5本の座標軸で位置関係が表現されている1〜16のサンプルについて(この図上でA〜Eを何か座標軸の矢印のように見なせるわけではないですけど)、A〜Eの座標軸を主成分分析で2座標(第1主成分と第2主成分)に縮約した上で、第1主成分、第2主成分とそれぞれに対する主成分得点から、もとの位置関係どれぐらい再現されるかを見てみるという感じですね。


 やることは簡単で、主成分分析の関数にデータを入れて終了です……。
 まず教科書のデータを入力します。めんどくさいですが、ScanSnapで教科書を自炊してOCRをかけていたので、ほぼコピペで再現できましたw

# データの入力

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))  # 空白の作図デバイス上に、点ではなく対応するテキストをプロットする


f:id:midnightseminar:20140828233219p:plain


 さて、主成分得点は$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)


f:id:midnightseminar:20140828233309p:plain


 すこしゆがんでますが、だいたい、丸い形に並びました!


 次に、別の面から主成分分析の性能を確かめるために、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)


f:id:midnightseminar:20140828233326p:plain

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) 


f:id:midnightseminar:20140828233357p:plain


 predict()関数は、モデルとパラメータを決めた(推定した)後に、そのモデルに説明変数の値を新たに代入して目的変数を計算し、予測を行うためのものですね。回帰分析の結果に新しい値を代入して予測するときにも使った気がする。


 さて次に、バイプロットを表示してみます。

biplot(oka.pc)


f:id:midnightseminar:20140829070725p:plain


 Biplotとはなんなのかというと、Wikipediaの解説によれば、二次元の座標軸上に、変数と値をいっしょに表示するグラフのこと。
 この図をみてると、主成分分析というのが何をやっている分析なのかがよく分かりますね。
 赤い矢印で表されているA〜Eは、もとの変数(標準化されてる)の主成分負荷量を表しており、左と下の座標軸をみればいい。つまりたとえばCというのは第1主成分とほぼ並行な変数なので、第1主成分に対する負荷量はほぼゼロとなっている。さっきの「loadings」の数字をみると実際そうなっています。
 1〜16の点は主成分得点でプロットされており、上と右の座標軸をみればいい。さっきのscoresの値をみると、実際そうなっています。


 続けて、Rの練習用データirisの4つの変数を2つに圧縮する主成分分析も行う。
 ちなみに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


 1〜4列目の変数を2変数に縮約して、主成分得点をプロットしますが、元データの5列目に入っているsetosa、versicolor、virginicaという文字列を、1〜3の数字に変えてラベルとして用います。

iris.pca <- princomp(iris[,-5])
plot(iris.pca$scores, type="n")
lab <- as.numeric(iris[,5])
text(iris.pca$scores, label=lab, col=lab)


f:id:midnightseminar:20140828233717p:plain


 バイプロットも出してみる。

biplot(iris.pca)


f:id:midnightseminar:20140828233733p:plain


 ごちゃごちゃして何のことかわかりませんね。
 次に、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]))


f:id:midnightseminar:20140829070704p:plain


 これは、主成分負荷量を軸にして元の変数の位置関係を表現している。

> s.class(iris.dudi$li, iris[,5], cpoint=1)


f:id:midnightseminar:20140829070650p:plain


 5列目の変数でまとめた上で、主成分得点の分布を見ている図ですね。何が中心になってるのかは知らない。

因子分析

 次は因子分析です。これも教科書の演習をそのまま実行してみる。
 まずデータをつくる。

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))  # 教科書には書いてないけど2つのグラフを横に並べるときはこうする
barplot(seiseki.fac$loadings[,1], col="lightblue")  #第1因子
barplot(seiseki.fac$loadings[,2], col="lightblue")  #第2因子


f:id:midnightseminar:20140829070628p:plain


 理系科目に反応している因子と、文系科目に反応している因子であることが分かりますね。
 しかし個人的には、結局、

write.csv(seiseki.fac$loadings, "seiseki.csv")


 とやって負荷量の行列をCSVに出して、エクセルでセルに色を付けるとかしたほうがみやすい気がします。まぁ、グラフの方が比率は直感的に分かるのでいいんですが、項目が増えてくると、因子負荷行列をエクセルのフィルターで並べ替えたりしながら見るのがいいんじゃないかなと。心理尺度なんかだと、どの項目を分析から削るという判断にも使えるし。
 エクセルだと文字化けしますがめんどくさいので文字コードの変換とかしてません。


f:id:midnightseminar:20140828233932p:plain


 最後に、バイプロットを出してみる。

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)


f:id:midnightseminar:20140828233950p:plain