Rで統計学を学ぶ(9)

この講義では、教科書の第16章「因子分析」をとりあげます。 因子分析とは、複数の変数の値において共通して影響する要因を仮定し、 その要因によって変数間の相関関係を説明しようとする手法です。

このウェブページと合わせて教科書を読み進めてください。

学習項目です:

因子分析
Rを用いて説明する。例では、200名のテスト得点(国語、社会、数学、理科、英語)のデータを取り上げる。それぞれの科目のテスト得点が観測変数である。 そして、人によって理系の科目が得意なもの(理系因子)と文系の科目が得意なもの(文系因子)があると想定し、その2つの因子によって成績データを説明しようと試みます。
[成績データ例の作成方法について]
教科書p.300の成績データは、p.308にある次のコードにより生成されます: 
set.seed(9999)
n <- 200
因子負荷行列 <- matrix(c(0.09884, 0.17545, 0.52720, 0.73462, 0.45620, 0.72141, 0.47258,
	 0.17901, 0.07984, 0.37204), nrow=5)
独自性 <- diag(sqrt(c(0.530201, 0.254119, 0.309986, 0.546036, 0.346539)))
因子得点 <- matrix(rnorm(2*n),nrow=2)
独自因子 <- matrix(rnorm(5*n),nrow=5)
五教科 <- round(t(因子負荷行列 %*% 因子得点 + 独自性%*%独自因子)*10+50)
colnames(五教科) <- c("国語","社会","数学","理科","英語")


> 相関行列 <- cor(五教科)
> 相関行列
          国語      社会      数学      理科      英語
国語 1.0000000 0.5500166 0.1953799 0.1630274 0.4275140	# 国語と社会の相関係数が0.55と高い
社会 0.5500166 1.0000000 0.3317530 0.2944938 0.5178159
数学 0.1953799 0.3317530 1.0000000 0.5301135 0.4575891	# 数学と理科の相関係数が0.53と高い
理科 0.1630274 0.2944938 0.5301135 1.0000000 0.3876493
英語 0.4275140 0.5178159 0.4575891 0.3876493 1.0000000	# 英語は他の教科との相関係数が0.39～0.52と全般的に高い
相関行列を求めてみると、国語と社会、数学と理科の相関係数が高いことがわかる(英語はどれとも相関する)。これにより文系、理系という因子の仮定は、それほど間違いではなさそう。。。そこで、因子の数を２（理系因子と文系因子）として因子分析を行う。 (因子分析を実際に使う場面では、因子の数をいろいろ変えてみて、最も妥当な解釈(説明)が可能な因子数を採用します)
[因子の数を決める手がかり: 固有値]
固有値を参考にして因子数を決めるのもよい方法です。 Rではeigen関数を用いて簡単に固有値と固有ベクトルを求められます: 
>  eigen(相関行列)
$values
[1] 2.5573782 1.0656034 0.5058941 0.4461472 0.4249772

$vectors
           [,1]        [,2]       [,3]        [,4]        [,5]
[1,] -0.4040025  0.57915506 -0.3526733 -0.30988872  0.53004894
[2,] -0.4791362  0.36314955 -0.1046695  0.07323541 -0.78881668
[3,] -0.4380493 -0.48395278  0.2494573 -0.71311577 -0.05603054
[4,] -0.4064807 -0.54402387 -0.6062633  0.39786745  0.11383232
[5,] -0.5000967  0.05029462  0.6594556  0.48142803  0.28411115
ここで values　が固有値、 vectorsがそれぞれの固有値に対する固有ベクトルです。 ちなみに、そうであることは次の計算から確かめられます: (Rでは一次元ベクトルは一行で表示されます) 
>  result <- eigen(相関行列)
>  (相関行列 %*% result$vectors[,1])[,1] 
     国語      社会      数学      理科      英語 
-1.033187 -1.225333 -1.120258 -1.039525 -1.278936 
>  result$values[1] * result$vectors[,1] 
[1] -1.033187 -1.225333 -1.120258 -1.039525 -1.278936
>  (相関行列 %*% result$vectors[,2])[,1] 
     国語      社会      数学      理科      英語 
 0.61714961  0.38697340 -0.51570173 -0.57971368  0.05359412 
>  result$values[2] * result$vectors[,2] 
[1]  0.61714961  0.38697340 -0.51570173 -0.57971368  0.05359412
ところで、固有ベクトルとか固有値とはどういう意味でしたでしょうか。 定義は線形代数の講義を思い出してもらうとして、 固有ベクトルと固有値についてここでの重要なポイントは次のことです:
行列Aの固有ベクトルをx、固有値をλとし、 行列Aをベクトルxに作用させる（行列の意味で掛け算する）とは、 幾何的には、ベクトルxに対し線形変換（拡大縮小や回転） を施していているとみなせる
つまり、固有ベクトルとは、行列Aを作用しても向きが変わらないもの、といえます。 因子分析の『文脈に即して言えば』、 行列Aを特徴づける成分の一つがxであるとみなせます。 また固有値は行列Aを作用した時の拡大率になりますが、 因子分析の『文脈に即して言えば』、 固有ベクトルxが行列Aを特徴づける成分での割合を表すものとみなせます。
そういう見方からeigen(相関行列)の結果を見ると 1以上の固有値は2個なので、この例では因子数を2としたと言えます。 

> 五教科因子分析 <- factanal(五教科, factors=2)
> 五教科因子分析

Call:
factanal(x = 五教科, factors = 2)	# factanalが因子分析の関数、factorsの値が因子の数

Uniquenesses:				# 独自性
 国語  社会  数学  理科  英語 
0.471 0.395 0.379 0.548 0.491 

Loadings:				# 因子負荷
     Factor1 Factor2			# 第1因子、第2因子
国語 0.722          
社会 0.730   0.268  
数学 0.177   0.768  
理科 0.156   0.654  
英語 0.537   0.470  

               Factor1 Factor2			# 第1因子、第2因子
SS loadings      1.399   1.317			# 因子負荷の平方和(因子寄与)
Proportion Var   0.280   0.263			# 分散説明率(因子寄与率)
Cumulative Var   0.280   0.543			# 累積因子寄与率
Test of the hypothesis that 2 factors are sufficient.		# 適合度判定
The chi square statistic is 0.08 on 1 degree of freedom.
The p-value is 0.782 			# p値が小さい場合は「不適合」
この出力の最後に表示されたp値は、「適合度の検定」の結果を表しています。 この値は、いま試した「因子数=2」というモデルの「適合度」の値で、 （今までの場合と異なり)p値が大きいほど「適合している」、 p値が小さいと「適合していない」ということを表します。 そして有意水準より小さい場合は「適合していない」として棄却されます。

共通性とは「1から独自性を引いた値」で、因子により説明できる部分を表します。 つまり、共通性の値が大きければ大きいほど、因子によって説明できる部分が大きいことを表します。 

> 共通性 <-  ( 1-五教科因子分析$uniquenesses )		# uniqunesses は小文字から始まることに注意
> 共通性
     国語      社会      数学      理科      英語 
0.5288198 0.6049597 0.6213085 0.4523362 0.5087881
この例では、社会と数学が(それぞれ第1因子と第2因子に対応して)共通性の値が大きく、理科がほかと比べて小さいという結果が得られました。

また『因子負荷』とは、それぞれの観測変数による因子の貢献度を表す数値です。 因子1に対して国語と社会が大きいアタになっているのに対し、 因子2に対しては数学と理科が大きな値になっています。 これらにより「観測変数は因子2個のモデルで説明できる。因子1は国語と社会、 因子2は数学と理科に共通する要因とかんがえられる(英語はその両方）」 と解釈されます（これらに対して「文系因子」とか「理系因子」と名付けるのはセンスが必要でしょう）。

ここで、先の結果で、国語の行においてFactor2の値が表示されていないことについて補足しておきます。 値が表示されなかったのは、その値が小さかったからです。これを表示させるには、つぎのように cutoff=0 を指定します: 

> print(五教科因子分析,cutoff=0)

Call:
factanal(x = 五教科, factors = 2)

Uniquenesses:
 国語  社会  数学  理科  英語 
0.471 0.395 0.379 0.548 0.491 

Loadings:
     Factor1 Factor2
国語 0.722   0.084  
社会 0.730   0.268  
数学 0.177   0.768  
理科 0.156   0.654  
英語 0.537   0.470  

               Factor1 Factor2
SS loadings      1.399   1.317
Proportion Var   0.280   0.263
Cumulative Var   0.280   0.543

Test of the hypothesis that 2 factors are sufficient.
The chi square statistic is 0.08 on 1 degree of freedom.
The p-value is 0.782
この結果の解釈については教科書を読んでください。

斜交回転

因子分析のためのRの関数 factanalでは、指定をしないと、最尤推定法により因子分析が行われ、直交回転によるバリマックス(Varimax)回転が行われます。「直交」とは、因子が共通の要素を持たない(ベクトルとして直交する)ことを意味します。

しかし時には、共通の要素を持つことを許すことで、より少ない因子により観測変数をよりよく説明することができる場合があります。そのために「斜交」回転をします。プロマックス(promax)回転は、斜交回転の有名なものの一つです。プロマックス回転を実行するには以下のようにします(rotationというオプションに注意): 

> 五教科因子分析斜交 <- factanal(五教科,factors=2, rotation="promax")
> print(五教科因子分析斜交,cutoff=0)

Call:
factanal(x = 五教科, factors = 2, rotation = "promax")

Uniquenesses:					# 独自性
 国語  社会  数学  理科  英語 
0.471 0.395 0.379 0.548 0.491 

Loadings:					# 因子負荷
     Factor1 Factor2				# 第1因子、第2因子
国語  0.801  -0.156 
社会  0.749   0.050 
数学 -0.051   0.815 
理科 -0.038   0.692 
英語  0.461   0.348 

               Factor1 Factor2			# 第1因子、第2因子
SS loadings      1.419   1.291			# 因子負荷の平方和(因子寄与)
Proportion Var   0.284   0.258			# 分散説明率(因子寄与率)
Cumulative Var   0.284   0.542			# 累積因子寄与率

Factor Correlations:				# 因子間相関
        Factor1 Factor2
Factor1   1.000   0.546
Factor2   0.546   1.000

Test of the hypothesis that 2 factors are sufficient.		# 適合度判定
The chi square statistic is 0.08 on 1 degree of freedom.
The p-value is 0.782
[回転なしの因子分析と斜交回転]
本文で述べたように、因子分析のためのRの関数 factanalでは、指定をしないと、最尤推定法による、また直交回転によるバリマックス回転による因子分析が行われます。そこで、rotation="promax"とするとプロマックス回転による因子分析になるのでした。
ここで、rotation="none"と指定すれば、「回転なし」の因子分析になります: 
> 五教科因子分析回転なし <- factanal(五教科, factors=2, rotation="none")
> print(五教科因子分析回転なし, cutoff=0)

Call:
factanal(x = 五教科, factors = 2, rotation = "none")

Uniquenesses:					# 独自性
 国語  社会  数学  理科  英語 
0.471 0.395 0.379 0.548 0.491 

Loadings:					# 因子負荷
     Factor1 Factor2				# 第1因子、第2因子
国語  0.583  -0.435 
社会  0.714  -0.307 
数学  0.657   0.436 
理科  0.563   0.368 
英語  0.713  -0.028 

               Factor1 Factor2			# 第1因子、第2因子
SS loadings      2.106   0.611			# 因子負荷の平方和(因子寄与)
Proportion Var   0.421   0.122			# 分散説明率(因子寄与率)
Cumulative Var   0.421   0.543			# 累積因子寄与率

Test of the hypothesis that 2 factors are sufficient.		# 適合度判定
The chi square statistic is 0.08 on 1 degree of freedom.
The p-value is 0.782
「回転なし」の因子分析の結果を初期解と呼びます。この初期解を元に、いろいろな回転をさせて最適な因子分析を求めることができます(探索的因子分析)。例えばプロマックス回転による因子分析を初期解から得るには次のようにpromax関数を使います: 
> 初期解を元にプロマックス回転 <- promax(loadings(五教科因子分析回転なし),m=4)
> print(初期解を元にプロマックス回転)
$loadings

Loadings:					# 因子負荷
     Factor1 Factor2				# 第1因子、第2因子
国語  0.801  -0.156 
社会  0.749         
数学          0.815 
理科          0.692 
英語  0.461   0.348 

               Factor1 Factor2			# 第1因子、第2因子
SS loadings      1.419   1.291			# 因子負荷の平方和(因子寄与)
Proportion Var   0.284   0.258			# 分散説明率(因子寄与率)
Cumulative Var   0.284   0.542			# 累積因子寄与率

$rotmat
           [,1]      [,2]
[1,]  0.6062459 0.5303245
[2,] -1.0284319 1.0695617
ただ、これだと、factanal関数を用いた時に表示されていた「因子間相関」が表示されません。それを得るには、 次のようにします: 
> 因子間相関 <- solve(t(初期解を元にプロマックス回転$rotmat)%*%初期解を元にプロマックス回転$rotmat)
> 因子間相関
          [,1]      [,2]
[1,] 1.0000000 0.5462117
[2,] 0.5462117 1.0000000

参考: 主成分分析

因子分析と「似た」分析に、主成分分析がある。主成分分析の基本的な狙いは、 「多変量データの持つ情報を少数の総合特性値（合成変数）に要約する」 (青木(2009)『Rによる統計解析』オーム社)である。 （さらに言えば、この合成変数の性質として、 相関が０，また分散の大きいものが情報量が多く、 小さいものは情報量が少ない‐したがって無視できる、と考える）
次も参考にすること:中部大学小塩先生による「因子分析と主成分分析の違い」
本項で取り上げた成績データに対して主成分分析をしてみよう: 
> analysis <- prcomp(五教科, scale=TRUE)
> print(analysis)
Standard deviations:              # 固有値の平方根
[1] 1.5991805 1.0322807 0.7112623 0.6679425 0.6519027

Rotation:                         # 固有ベクトル
            PC1         PC2        PC3         PC4         PC5
国語 -0.4040025  0.57915506 -0.3526733  0.30988872  0.53004894
社会 -0.4791362  0.36314955 -0.1046695 -0.07323541 -0.78881668
数学 -0.4380493 -0.48395278  0.2494573  0.71311577 -0.05603054
理科 -0.4064807 -0.54402387 -0.6062633 -0.39786745  0.11383232
英語 -0.5000967  0.05029462  0.6594556 -0.48142803  0.28411115

> summary(analysis)
Importance of components:
                          PC1    PC2    PC3     PC4    PC5
Standard deviation     1.5992 1.0323 0.7113 0.66794 0.6519
Proportion of Variance 0.5115 0.2131 0.1012 0.08923 0.0850
Cumulative Proportion  0.5115 0.7246 0.8258 0.91500 1.0000
主成分分析は、元データ（行列）から、固有値と固有ベクトルを求めることにより計算される。固有ベクトルPC1, ..., PC5は互いに相関のない変数ベクトルであり、 これらの固有ベクトルによって元のデータを記述することができる。 また、この固有ベクトルPC1からPC5は、その記述に対する貢献度によって順番がついており、それぞれの貢献度はsummary関数の結果のProportion of Varianceによって得られる。

上位2つの固有ベクトルによって散布図を書くと右図が得られる。 この図によって、五教科間の関係を考えることができなくはない。

また、上位2つの固有ベクトルの空間において個々のデータがどのように位置するかを図に書くこともできる(左図)。この図をバイプロット(biplot)といい、この図によって個々のデータの位置関係を考えることができる(番号は五教科データの番号＝学生番号に相当)。

各成分による主成分得点は次のようにして得られる: 

> analysis$x
            PC1        PC2         PC3        PC4         PC5
[1,] -2.1311183  0.6511116 -0.30343340  0.1437321 -0.79789374
[2,]  1.1396340 -0.9566057 -0.46724133 -0.5109947  0.03715720
[3,] -3.2538608 -0.2209607  1.09671671 -0.2844395 -0.48947287
[4,]  0.1904725  0.4888520 -0.74063544 -0.3089838  1.12269123
[5,] -2.5159078 -0.8648665 -0.38302419  0.1108199  0.14328172
[6,]  0.1511904  1.5852820 -0.01486356 -0.3087898 -0.01072046
…
ここで第1主成分(PC1)に注目すると、点数の高い3番と5番の学生の主成分得点がそれぞれ－3.25と-2.52、点数の低い2番の学生の主成分得点が1.14と表されている（ここでは値がマイナスの方が成績が良い）。

[上位2つの固有ベクトルの散布図とバイプロットの書き方]
散布図は次のコードによって得られた: 
plot(analysis$rotation[,1],analysis$rotation[,2],type="n")
text(analysis$rotation[,1],analysis$rotation[,2],colnames(五教科))
またバイプロットは次のコードによって得られた: 
biplot(analysis)

演習問題

問題9-1.

五教科の例題にたいし、初期解を求め、初期解から関数varimaxを用いて、バリマックス回転による因子分析の結果を求めよ(初期解とは、「回転なし」の因子分析の結果のこと)。
またこの結果と、promax関数において初期解を用い、mの値を2～4と変えた結果とを比較せよ。
注意: promax関数のmパラメタは「乗数」を表す。詳しくは教科書p.307を見る こと。（「回転なしの因子分析と斜交回転」の項も参照のこと)

問題 9-2.

東京図書のウェブサイトから 豊田(2012)「因子分析入門――Rで学ぶ最新データ解析」のデータである DL02126.zipをダウンロードせよ。その中の YG3.csv ファイルの内容に対して因子分析を行え(プログラムの結果だけではなく「観測変数の関係の説明」も試みること)。ここで因子数は7までの適当な数を試すこと。
なお、これはY-G性格検査(矢田部・ギルフォード性格検査」による980人、12項目(本来は120の項目)の調査結果である。記号の意味は次の通り: D: 抑うつ性（悲観的気分)、C: 回帰性傾向(気分の変りやすさ)、I:劣等感(自己の過小評価)、N: 神経質(敏感さ、心配性、傷つきやすさ)、O:客観性の欠如(空想性、主観的態度)、Co: 協調性の欠如(疑い深さ、不信感)、Ag: 愛想の悪さ(衝動性、攻撃性)、G:一般的活動性(心身両面での活発さ)、R:のんきさ(気軽な性質)、T: 思考的外向(瞑想的、反省的態度とは反対の傾向)、A:心配性(社会的リーダーシップ、服従性)、S:社会的外向(社会的接触を好む)
詳細については 豊田(2012)「因子分析入門――Rで学ぶ最新データ解析」p.32を参照のこと。
Rで統計学(8)　　　　　　　　　トップページに戻る