Rで因子分析 基礎

この原稿は統計学勉強会のためにかかれたもので、学校で習うあたりまでの簡単な解説である。

因子分析とは

(探索的)因子分析とは、観測できる変数から観察できない因子(潜在変数)を見出す探索し発見する方法である。このように説明してもよくわからないので、実際に実例を使いながら説明をした方がよさそうである。

Rで因子分析をする際にはpsychパッケージを使用する。psychパッケージにはサンプルデータも入っているため、入門にもちょうどいいと思う。

library("psych")
library("GPArotation")

GPArotationパッケージは回転をするためのパッケージ。 分析の際に使用するので一緒に読み込んでおく。

psychパッケージにはパーソナリティの5因子モデル(IPIP-NEO)のデータが含まれている。 パーソナリティの5因子モデルのデータは"bfi"という名前である。

data(bfi)

5因子モデルの議論にはあまり詳しくないが、IPIP-NEOというのは、人のパーソナリティを5つの因子で分ける試みの一つである(https://sapa-project.org/)。パーソナリティの5因子モデルの試みはいくつもあり、精神医学ではNEO-PI-Rを見かけることがあるが、サンプルデータなので細かいところは措いておこう。

調べたところIPIP-NEOには日本語の定訳がないようだ。変に日本語に翻訳してしまうと、正式な日本語訳だと誤解されると迷惑をかけるので、質問文は原版通りに、それぞれの因子の邦訳だけしておく。

Agreeableness 協調性・調和性
A1 Am indifferent to the feelings of others. (q_146)
A2 Inquire about others’ well-being. (q_1162)
A3 Know how to comfort others. (q_1206)
A4 Love children. (q_1364)
A5 Make people feel at ease. (q_1419)

Conscientiousness 勤勉性・誠実性
C1 Am exacting in my work. (q_124)
C2 Continue until everything is perfect. (q_530)
C3 Do things according to a plan. (q_619)
C4 Do things in a half-way manner. (q_626)
C5 Waste my time. (q_1949)

Extraversion 外向性
E1 Don’t talk a lot. (q_712)
E2 Find it difficult to approach others. (q_901)
E3 Know how to captivate people. (q_1205)
E4 Make friends easily. (q_1410)
E5 Take charge. (q_1768)

Neuroticism 神経症傾向・情緒〔不〕安定性
N1 Get angry easily. (q_952)
N2 Get irritated easily. (q_974)
N3 Have frequent mood swings. (q_1099)
N4 Often feel blue. (q_1479)
N5 Panic easily. (q_1505)

Openness 開放性・経験への開放性
O1 Am full of ideas. (q_128)
O2 Avoid difficult reading material.(q_316)
O3 Carry the conversation to a higher level. (q_492)
O4 Spend time reflecting on things. (q_1738)
O5 Will not probe deeply into a subject. (q_1964)

gender 性別 Males = 1, Females =2
education 教育 1 = HS 高校入学, 2 = finished HS 高卒, 3 = some college 大学入学, 4 = college graduate 大卒 5 = graduate degree 院卒
age 年齢 age in years

この質問文の中から5つの「まとまり」(因子)を導くのか因子分析の目的である。 データは確証的因子分析のものだが、今回はそのデータを探索的因子分析を用いる。

d1 <- bfi[1:25]

26列目以降は性別・教育・性別のデータが入っているので、1-25列だけを使用する。

res01 <- fa(d1, nfactors = 5, fm = "ml", rotate = "promax", scores=TRUE)
print(res01, digits = 3, sort=TRUE)

データ: d1、因子数: nfactors、推定方法: fm、回転方法: rotate、因子得点を表示scores=TRUEで表記する。データを最初に持ってくる以外は順番が入れ替わっても問題ない。

printのオプションは小数点を第3位まで: digits = 3、ソートを指定する。 以下が因子負荷(パターン行列)である。

ここではプロマックス回転、最尤法を使用する。最もよく使われる回転方法と推定法のペアである。最尤法には問題はない。データが正規分布から大きく逸脱している場合には、最尤法がむしろ望ましい。問題は、プロマックス回転の方だろう。プロマックス回転は簡便法であり、必ずしも良い選択肢とはならない。代わりにジオミン回転や独立クラスタ回転(ハリス・カイザー回転)などを使うのが良いだろう。

Standardized loadings (pattern matrix) based upon correlation matrix
  item    ML2    ML1    ML3    ML5    ML4    h2    u2  com
N1   16  0.862  0.224  0.041 -0.266 -0.097 0.705 0.295 1.37
N2   17  0.834  0.168  0.054 -0.247 -0.036 0.657 0.343 1.28
N3   18  0.721 -0.003 -0.012 -0.019 -0.007 0.525 0.475 1.00
N5   20  0.502 -0.159  0.026  0.147 -0.150 0.338 0.662 1.59
N4   19  0.492 -0.321 -0.102  0.054  0.084 0.478 0.522 1.93
E2   12  0.119 -0.709  0.046 -0.019 -0.005 0.548 0.452 1.07
E1   11 -0.055 -0.647  0.159 -0.015 -0.044 0.369 0.631 1.15
E4   14 -0.014  0.608 -0.043  0.274 -0.113 0.519 0.481 1.48
E5   15  0.158  0.480  0.237 -0.008  0.167 0.405 0.595 2.00
E3   13  0.075  0.480 -0.067  0.238  0.251 0.441 0.559 2.15
C2    7  0.141 -0.114  0.662  0.076  0.062 0.427 0.573 1.20
C4    9  0.152  0.029 -0.644  0.053 -0.039 0.465 0.535 1.14
C3    8  0.039 -0.098  0.585  0.088 -0.058 0.317 0.683 1.13
C5   10  0.187 -0.090 -0.562  0.015  0.093 0.435 0.565 1.34
C1    6  0.058 -0.050  0.542 -0.014  0.158 0.321 0.679 1.21
A3    3  0.036  0.203  0.011  0.621 -0.008 0.511 0.489 1.22
A2    2  0.048  0.099  0.082  0.567 -0.005 0.401 0.599 1.12
A5    5 -0.089  0.278 -0.039  0.546  0.013 0.483 0.517 1.56
A4    4 -0.024  0.075  0.196  0.425 -0.169 0.286 0.714 1.84
A1    1  0.151  0.106  0.044 -0.394 -0.036 0.150 0.850 1.50
O3   23  0.037  0.239 -0.037  0.080  0.594 0.473 0.527 1.37
O5   25  0.096  0.025 -0.025  0.052 -0.522 0.274 0.726 1.10
O1   21  0.015  0.155  0.028  0.025  0.508 0.324 0.676 1.20
O2   22  0.174  0.022 -0.075  0.153 -0.450 0.244 0.756 1.62
O4   24  0.158 -0.257 -0.020  0.170  0.381 0.257 0.743 2.627

表の見る部分は、数字の大小である。 因子負荷 は1~-1までの値をとり、1もしくは-1に近ければ非常に関連が高いという意味で、0が関連が低いという意味である。相関係数と同様の解釈の方法でよい。ちなみに、直交回転の場合因子負荷と相関係数は一致するが、斜交回転では一致しないので、斜交回転の際には、値が異なることに注意が必要である。

項目N1のML2(第2因子)は0.862であり、関連が高いことがわかる。他のN項目も同様にML2と正の関連があることがわかる。ML3は0.041と関連が低くなっており、N項目と第3因子(ML3)には関連がないことが確認できる。

printコマンドで"sort=TRUE"にすると自動で並び変えてくれるため、解釈が非常にやりやすくなる。ただし、試行錯誤の段階では、この例のように明確に分かれないため、ソートをしない方が解釈しやすいということもあるので、適宜使い分けが必要だ。

原因はわからないが、RStudioのRMarkdownからでは因子得点は出力できなかった。RStudioのコンソールから直接コマンドを入れると出力される。ノーマルなRであれば、このような問題は起きないので、RStudioを使用する際にだけ注意が必要だ。

次のブロックに表示されるのは統計量である。

                      ML2   ML1   ML3   ML5   ML4
SS loadings           2.651 2.389 1.944 1.875 1.493
Proportion Var        0.106 0.096 0.078 0.075 0.060
Cumulative Var        0.106 0.202 0.279 0.354 0.414
Proportion Explained  0.256 0.231 0.188 0.181 0.144
Cumulative Proportion 0.256 0.487 0.675 0.856 1.000

SS loadings: 因子負荷量平方和
Proportion Var: 寄与率
Cumulative Var: 累積寄与率
Proportion Explained: 説明率
Cumulative Proportion: 累積説明率

論文に書くときにはは、因子負荷量平方和、寄与率、累積寄与率の3つがあればよい。

次のブロックは因子間の相関が表示される。

With factor correlations of
      ML2    ML1    ML3   ML5   ML4
ML2  1.000 -0.272 -0.231 0.014 0.025
ML1 -0.272  1.000  0.375 0.335 0.153
ML3 -0.231  0.375  1.000 0.243 0.211
ML5  0.014  0.335  0.243 1.000 0.187
ML4  0.025  0.153  0.211 0.187 1.000

ほとんどの相関係数が絶対値0.35以下であり、5つの因子間の相関がそれほど高くないことがわかる。

次のブロックが最後で適合度などが表記される。

The degrees of freedom for the null model are  300  and the objective function was  7.228 with Chi Square of  20163.79
The degrees of freedom for the model are 185  and the objective function was  0.628

The root mean square of the residuals (RMSR) is  0.03
The df corrected root mean square of the residuals is  0.038

The harmonic number of observations is  2762 with the empirical chi square  1474.696  with prob <  1.29e-199
The total number of observations was  2800  with Likelihood Chi Square =  1749.883  with prob <  1.39e-252

Tucker Lewis Index of factoring reliability =  0.8721
RMSEA index =  0.0551  and the 90 % confidence intervals are  0.0526 0.0573
BIC =  281.469
Fit based upon off diagonal values = 0.979
Measures of factor score adequacy             
                                                    ML2   ML1   ML3   ML5   ML4
Correlation of (regression) scores with factors   0.929 0.904 0.881 0.868 0.835
Multiple R square of scores with factors          0.863 0.816 0.777 0.754 0.698
Minimum correlation of possible factor scores     0.727 0.633 0.553 0.508 0.395

スクリープロット法とカイザーガットマン基準

因子数を決定するために一般的に使われている方法。この方法は他の方法で代用した方がよいが、普及率が高いので記載しておく。

scree(d1,factors=TRUE,pc=TRUE,main="Scree plot",hline=NULL,add=FALSE)
VSS.scree(d1, main = "scree plot")  

f:id:iDES:20190319093630p:plain

スクリープロット法は図にしてガクッと落ちているところを探して、落ちる前の数を因子数として採用する方法である。 この図では、明確にわからない。2~3でガクッと落ちているものの、固有値はそこまで低くなく、ここで切っていいかは判断しづらい。

次に一般的によく使われるのが、カイザーガットマン基準である。固有値を「1」と決めて、それ以上のものを採用するという基準である。

図の1のところに線が横に引いてあるのが、カイザーガットマン基準である。カイザーガットマン基準に従えば、6因子構造を採用することになる。

ただ、このデータは5因子モデルなので、カイザーガットマン基準の6因子という結果とは異なる。また、スクリープロット法もあまりよくわからないので、探索的因子分析で、このデータをこの2つの方法で、5因子と決定することは難しい。

MAP基準、平行分析、BICの使用が現在では推奨されており、上の2つの方法は参考程度にするのがよいだろう。

クロンバッハのα係数

クロンバッハのαとは内的妥当性の係数の一つである。信頼性係数とも呼ぶ。外的妥当性とは計測したいものをどのくらい計測できているか、他の方法と比較してどのくらい妥当なのかをみるが、内的妥当性は質問項目が因子を導くのにどのくらい妥当かというものである。質問項目にノイズが多くあれば、クロンバッハのαは低くなる。

クロンバッハのαはpsychパッケージに含まれており、以下のコマンドで出力される。

alpha(d1)

各項目の統計量まで出力されるが、論文で使用する際には、以下値さえあればよい。

lower alpha upper   95% confidence boundaries
0.51  0.53  0.56

使用するのは中央の0.53である。95%信頼区間が表示されているが、この値の表記を求められることはないと思う。

0.53という数字だが、これは高いとは言えない数字だと思う。表記が主観的になっているので、いくつであれば良いかというのが明確に決まっていないためである。因子分析を使用している人たちにも意見のばらつきはあるだろうが、0.7以上くらいは求められるのではないだろうか。

クロンバッハのαず最も使用頻度の高い方法だが、もう少し良い指標である「ω係数」などもあるため、そちらの表記の方が望ましい。ただ、クロンバッハのαがまだ主流なので、ω係数とは何か?と質問がでる所(例えば学会など)では、両方の併記が望ましいだろう。

この記事はRでの基本的な因子分析の方法しか扱っていない。 残したお題は、次回研究会で扱うので、その際にまた記事を書くつもりである。