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

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

Rの練習: 因子分析の分析対象から除去すべき項目を割り出すプログラムを書いてみる

 心理学で心理測定尺度の因子分析を行う際に、30項目なら30項目の質問を並べて「とても当てはまる」〜「全く当てはまらない」までの7段階の回答を取り、7点〜1点を割り振って、因子分析を行うことがよくあります。


 で、たとえば仮説によって「3因子」の構造が想定される(あるいはスクリープロットなどから3因子構造が示唆される)として、分析してみたときに因子1と因子3の両方に高い負荷をもつ質問項目があったりすると、尺度の構成上困ってしまいます。
 そういう項目は、「ああ、質問文の作りが良くなかったんだな」と判断して分析対象から除いて、再度因子分析を行ったりします。
 また、各質問項目の「共通性」をみて、著しく低い値が出ている場合は、共通因子によって説明される分散の割合が小さいということなので、その項目も除いてみたりします。


 たとえば、手元にある調査データの値を適当に変換してつくったダミーデータで以下のように分析してみます。

> # dというダミーデータに,Q1〜Q30までの質問に対する,100人分の回答結果が入っている.
> ncol(d)
[1] 30
> nrow(d)
[1] 100
> d[1:5, 1:10]  # 途中まで表示
  Q1 Q2 Q3 Q4 Q5 Q6 Q7 Q8 Q9 Q10
1  5  5  4  4  3  4  5  4  4   2
2  5  1  1  1  1  3  1  1  1   2
3  5  5  5  4  5  6  5  5  5   3
4  6  5  4  4  4  5  4  4  4   4
5  6  6  6  5  5  5  5  6  5   6
> 
>
> # 3因子を仮定し,最尤法,プロマックス回転と設定して因子分析を実行
> fa1 <- fa(d, nfactors=3, rotate="promax", fm="ml", normalize = TRUE)
>
> # 因子負荷行列と共通性を表示
> cbind(fa1$loadings, Communality=fa1$communality)
             ML1          ML2          ML3 Communality
Q1   0.440092919  0.252129385  0.106704526  0.42186495
Q2   0.877347598  0.068420081 -0.075151917  0.79029623
Q3   0.880892066  0.029174292 -0.050631633  0.77340139
Q4   0.906008329 -0.217738063  0.065185081  0.71531032
Q5   0.940186457 -0.245670297  0.083483863  0.77095367
Q6   0.473027956  0.027416496  0.316386553  0.44121911
Q7   0.621496971  0.132040332  0.124220827  0.55831782
Q8   0.857372790 -0.036850618 -0.007835338  0.70224750
Q9   0.456445572  0.283011683 -0.145718584  0.36259502
Q10  0.182259105  0.388044735 -0.418319727  0.27399807
Q11  0.221939958  0.532855190 -0.108466156  0.40518757
Q12 -0.141619263  0.542818149  0.138319104  0.29488552
Q13 -0.004335725  0.621561660 -0.002386915  0.38288265
Q14  0.163568311  0.293328042 -0.342374714  0.17573429
Q15 -0.077084287  0.659286839 -0.171556066  0.35972784
Q16 -0.267276869  0.696001074  0.067660795  0.40130417
Q17 -0.028763427  0.649654359 -0.048204926  0.38877515
Q18  0.089992972  0.401949430  0.152377105  0.27461881
Q19 -0.047480145  0.369660986  0.358691874  0.32262318
Q20 -0.227678321  0.459829590  0.064602554  0.17690950
Q21  0.154349220  0.317294219  0.405582098  0.45739175
Q22  0.171645227  0.370004504  0.143841491  0.29698209
Q23  0.002454247 -0.006396383  0.654273537  0.42654771
Q24 -0.059417465  0.239524253  0.363595043  0.21981782
Q25  0.137740454  0.007175661  0.689052487  0.56054958
Q26 -0.096950193  0.142295568  0.743313088  0.58764531
Q27 -0.079146909 -0.089308471  0.545175097  0.25919401
Q28 -0.176008243  0.574273844  0.407900013  0.53033114
Q29  0.290043810 -0.096281760  0.041824751  0.07411997
Q30  0.234860294 -0.054848044  0.590181381  0.46556020
> 


 因子負荷行列をみると、たとえばQ6なんかは第1因子に0.473、第3因子に0.316という負荷量を持っていて、「どっちやねん」と解釈しづらい項目になっています。
 また、共通性が0.1台などのものもあり、共通因子によって説明できているとはあまり言えない項目もたくさんまじっていますね。


 因子分析には、どういう因子構造として解釈するのが適切かという唯一の正解は存在しませんが*1、説明力の低い結果とか、複数の因子にかかってしまっている項目があるような結果になるのは、なるべく避けたいです。


 たとえば上の結果をみて、共通性0.3未満のやつをとりあえず削ってみるとか、2つ以上の因子に対して絶対値0.3以上の負荷量を持っている項目を削ってみるとかしてみて、もう1回因子分析をやってみたいですね。(最終的にたとえば共通性は0.5以上ないとダメとかは言われますが。)


 しかし、因子分析を1回行って、結果の因子負荷行列などをみて、「あーこの項目削ってみよう」とか言ってエクセルで削ったテーブルをつくり、再度読み込んでやり直す・・・といったことはしたくありません。
 なるべく自動的に、どの項目が削除候補となるかを計算し、それを除いて再度分析するといったところまで自動でやりたいですね。


 というわけで、以下のような関数を作成しました。
 パッケージpsychが読み込まれていることが前提です。(因子分析のfa()関数を使うため。)

############################################
# 因子分析で,除去すべき項目を割り出す関数
############################################

ToBeRem <- function (data, load, com, fact, rot, fm) {
    # data: 因子分析にかけるデータフレーム
    # load: 「複数因子への負荷」の基準を,負荷量絶対値いくら以上にするか
    # com: 「共通性不足」の基準はいくらにするか
    # 因子分析関数FAの引数: fact=因子数,rot=回転法,fm=推定法
    # dataに指定したDFを因子分析にかけ,2個以上の因子に対して負荷量の絶対値がload以上になる項目の番号及び名称と,共通性の値がcom未満になる項目の番号及び名称を,それぞれベクトルで返す.

FA <- fa(data, nfactors=fact, rotate=rot, fm=fm, normalize = TRUE)  # 因子分析を実行
loadings <- FA$loadings  # 因子負荷量の行列を取り出し

# 複数因子に負荷している質問項目を特定
multifactor.number <- c()
multifactor.name <- c()
for (i in 1:nrow(loadings)) {
    # 各共通因子への負荷量を,大きい順に並べたときの2番目がload以上になるか,
    # 小さい順に並べたときの2番目が-load以下になった場合,項目の番号及び名称を返す
x1 <- loadings[i,]
x2 <- -sort(-abs(x1))  #絶対値降順に並べる
if (x2[2] >= load) {
	multifactor.number <- c(multifactor.number, i)
	multifactor.name <- c(multifactor.name, rownames(loadings)[i])
   }
}

# 共通性が低い項目を特定

lowcommunality.number <- c()
lowcommunality.name <- c()
if (min(FA$communality) < com) {
low <- as.vector(which(FA$communality < com))
lowcommunality.number <- low
lowcommunality.name <- rownames(loadings)[low]
}


# 出力結果をリストにまとめる
list <- list(multifactor.number, multifactor.name, lowcommunality.number, lowcommunality.name)
names(list) <- c("multifactor.number", "multifactor.name", "lowcommunality.number", "lowcommunality.name")
return(list)
}


 戻り値の定義の仕方をあまり分かってないのですが、複数あるときはリストにしないとダメなようです。
 ではこれをつかって、最初のデータで分析を実行してみます。

> # 削除候補となる項目を判定!
> remlist <- ToBeRem(data=d, load=0.3, com=0.3, fact=3, rot="promax", fm="ml")
> print(remlist)
$multifactor.number
[1]  6 10 19 21 28

$multifactor.name
[1] "Q6"  "Q10" "Q19" "Q21" "Q28"

$lowcommunality.number
[1] 10 12 14 18 20 22 24 27 29

$lowcommunality.name
[1] "Q10" "Q12" "Q14" "Q18" "Q20" "Q22" "Q24" "Q27" "Q29"

> 


 おおww
 13項目(Q10は両方の基準にひっかかっている)も、削除候補となりました。
 これらを削ったデータを用意してもう1回分析してみます。

> # 判定された項目を除去したデータをつくる
> remnumber <- sort(unique(c(remlist$multifactor.number, remlist$lowcommunality.number)))  # unique()で重複を削る
> d.new <- d[,-remnumber]  # やり直し用のDFを作成
> ncol(d.new) # 生き残り項目数
[1] 17


 17項目残っていますね。

> print(remlist.new)
$multifactor.number
[1] 1 8

$multifactor.name
[1] "Q1" "Q9"

$lowcommunality.number
[1]  9 12 13

$lowcommunality.name
[1] "Q11" "Q16" "Q17"

> 


 また5つも削除候補になってしまいました。
 だめだこの分析は感じですね。(じつはスクリープロットや平行分析の結果からすると3因子ではなく4因子が示唆されている。)


 まぁ、実際の分析はここまで機械的に進めていくわけでもなく、もう少しいろんな要素を考慮して削除項目を決めていきますが、とりあえず判定のコードは書けたということで練習にはなりました。実際、使うと思います。

*1:いわゆる確証的、確認的因子分析の場合はどうかとかありますけど