主成分分析・因子分析
藤野秀則
初稿日:“2024-10-15” 更新日:2025-02-07
取得したデータの変数が多い場合（多変量データと呼ぶ）に，それら1つ1つの変数を分析するのは非常に骨が折れる作業になる．そういう場合には，「似たようなデータ」になっているものを集約して，もとの変数よりも少ない数の集約変数で分析するのが効率的である．
そのような変数の集約をするための分析が主成分分析（PCA: Principal Component Analysis）と因子分析（FA: Factor Analysis）である．
以下ではこれら2つの分析の仕方について説明する．
なお，例によってtidyverseパッケージは予め読み込んでおくこと．
1 主成分分析
主成分分析は多数の変数を持ったデータから，文字通りデータの「主成分」を抽出する分析である．主成分とは，もとのデータ全体のばらつき方の情報を可能な限り反映した変数のことである． 最も強く反映した変数を第1主成分，２番目に反映した変数は第2主成分，３番目を第3主成分…という具合に名付けられる．
では，実際に以下の例題を使って主成分分析を行なっていこう．
こちらのデータは，大学入試共通テストの5教科9科目分の結果である．これを主成分分析せよ．
（あくまで2024年度の結果からシミュレーションしたものであり，実際の得点ではない）
まずはこのデータの概要を把握しておこう．今回のデータは全部で1万点もあるデータなのでprint(data)とするととんでもない表が出力されてしまうので，冒頭6行だけを表示させるようにhead()を使っている．
また，データの概要として標準偏差も表示させたいのでsummary()ではなくpsychパッケージに含まれるdescribe()を使っている．describe()は，こちらで紹介をしたデータの要約統計量を出力する関数である．
data <- read.csv("./practice/example_13_Test.csv")
head(data)## 国語 数学IA 数学IIB 世界史B 地理B 物理 化学 英語R 英語L
## 1 119 58 42 51 54 66 58 47 67
## 2 39 38 20 58 50 46 55 17 60
## 3 122 36 37 47 75 62 36 26 67
## 4 194 63 43 102 74 119 108 83 110
## 5 121 55 35 88 74 69 74 74 81
## 6 128 75 49 83 77 83 57 66 77psych::describe(data)## vars n mean sd median trimmed mad min max range skew
## 国語 1 10000 116.07 35.47 117 116.20 35.58 -23 257 280 -0.04
## 数学IA 2 10000 51.24 21.41 51 51.27 20.76 -21 138 159 -0.01
## 数学IIB 3 10000 57.69 26.25 58 57.62 26.69 -48 151 199 0.01
## 世界史B 4 10000 60.19 23.84 60 60.19 23.72 -40 168 208 0.01
## 地理B 5 10000 65.73 16.65 66 65.77 16.31 5 126 121 -0.04
## 物理 6 10000 62.58 27.81 62 62.53 28.17 -45 175 220 0.03
## 化学 7 10000 54.42 23.29 54 54.34 23.72 -35 143 178 0.04
## 英語R 8 10000 51.33 22.08 51 51.30 22.24 -32 138 170 0.01
## 英語L 9 10000 67.15 22.30 67 67.19 22.24 -10 156 166 0.00
## kurtosis se
## 国語 -0.01 0.35
## 数学IA -0.06 0.21
## 数学IIB 0.02 0.26
## 世界史B 0.08 0.24
## 地理B -0.03 0.17
## 物理 0.03 0.28
## 化学 0.01 0.23
## 英語R 0.01 0.22
## 英語L 0.04 0.22なお，psych::describe()という書き方はpsychパッケージを，describe()関数を使うときだけ適用する場合に用いる書き方である．もちろんlibrary(psych)を実行してもよいが，今回はdescribe()を実行するためだけなので，このような対応をした．
::の詳細
::は正確には「名前空間」（name
space）を指定する演算子（Operator)である．
各パッケージはそれぞれそのパッケージ名と同一の「空間」を持っており，各パッケージに含まれる関数はその空間の中でその機能が定義されている．
library()関数は現在実行中のプロジェクトにパッケージを読み込む関数であるが，内部で行っていることは，そのパッケージの名前空間をプロジェクトと結びつけてる，というものである．そしてRで関数を実行する際には，その都度，そのプロジェクトと結びつけられた名前空間から関数の定義を検索し，その名前空間まで定義を読み込みに行くという処理を内部で行っている．
つまりlibrary()を使うことによって，その名前空間を関数の検察範囲を含めることができるので，その名前空間に含まれる関数を使用できるようになるのである．
::を使うと，関数を呼び出す時点で明示的に名前空間を指定することができる．こうすると，検索を掛ける必要がなく，直接その名前空間の中の定義を読み込みに行ける．このため，PCにそのパッケージが存在していさえすればlibrary()でそのパッケージを読み込んでいなくてもその関数を使うことができるのである．
::で明確にパッケージを指定することによって，期待通りの関数が呼び出されるようになる．
さらに各科目間の相関もついでに見ておこう．
round(cor(data),2)## 国語 数学IA 数学IIB 世界史B 地理B 物理 化学 英語R 英語L
## 国語 1.00 0.13 0.11 0.40 0.31 0.11 0.09 0.33 0.22
## 数学IA 0.13 1.00 0.60 0.08 0.09 0.48 0.36 0.26 0.32
## 数学IIB 0.11 0.60 1.00 0.08 0.10 0.52 0.37 0.20 0.23
## 世界史B 0.40 0.08 0.08 1.00 0.33 0.09 0.07 0.16 0.13
## 地理B 0.31 0.09 0.10 0.33 1.00 0.10 0.07 0.15 0.13
## 物理 0.11 0.48 0.52 0.09 0.10 1.00 0.25 0.16 0.20
## 化学 0.09 0.36 0.37 0.07 0.07 0.25 1.00 0.14 0.15
## 英語R 0.33 0.26 0.20 0.16 0.15 0.16 0.14 1.00 0.60
## 英語L 0.22 0.32 0.23 0.13 0.13 0.20 0.15 0.60 1.00この結果から，国語・世界史B・地理の間，また，数学IA・数学IIB・化学・物理の間，さらに英語R・英語Lの間にはある程度相関があることが分かる．
1.1 主成分分析の実行
では主成分分析を実行してみよう
主成分分析を行う関数はprcomp()である．引数は以下の通り．
- 第1引数: 主成分分析に掛ける多変量が含まれたデータフレーム
- 第2引数
scale.=TRUE: データを標準化するかどうか．省略するとFALSE(標準化しない)が設定される．通常はTRUEにすべきとされる．
第2引数について
第2引数についてはデフォルトでTRUEとなっている．これは例えば多変量データの中に，単位が異なる変数が混在している場合，ばらつき方の情報が単位に依存してしまい，主成分分析の結果が変わってしまうためである．例えば体重[kg]と身長[cm]の場合，体重の分散は身長の分散に比べて小さくなってしまう．これは実際に「ばらつき方」として小さいということではなく，体重が身長に比べて小さな値となる単位を用いているためである．逆に身長を[m]で表記した場合には，体重の方が身長よりも分散が大きくなる．このように，そのままの値ではばらつき方の情報がそれぞれの変数の単位に依存してしまうため，標準化を行って単位系の違いが影響しないようにして，真のばらつき方の情報を利用するのが一般的である．
結果は一旦オブジェクトPCAに収めた上で，summary()で出力させる．
pca <- prcomp(data, scale. = TRUE)
summary(pca)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.6964 1.2735 1.0751 0.88098 0.83799 0.7625 0.71576
## Proportion of Variance 0.3197 0.1802 0.1284 0.08624 0.07803 0.0646 0.05692
## Cumulative Proportion 0.3197 0.4999 0.6284 0.71460 0.79262 0.8572 0.91414
## PC8 PC9
## Standard deviation 0.63100 0.61202
## Proportion of Variance 0.04424 0.04162
## Cumulative Proportion 0.95838 1.00000主成分分析では，与えた変数の数だけ「主成分」が出力され，それらは順に第1主成分，第2主成分，第3主成分，…と呼ばれる．今回の例の場合だと9つの変数を与えているので出力される主成分は全部で9つとなる．
第1主成分が元のデータのばらつきをもっとも大きく反映している変数であり，第2主成分は第1主成分に含まれていないばらつきを反映している変数，第3主成分は第1主成分と第2主成分に含まれていないばらつきを反映している変数となる．こうしたことから，9つ全部使えば，当然ながら元のデータのすべてのばらつきを反映することができる．
summary()関数で得られる結果は，合成された変数の元のばらつきを反映している度合についての情報を，その度合いが大きいものから順に出力する．
- 1段目: 主成分として抽出（正確には合成）されたそれぞれの成分の標準偏差
- 2段目: 分散説明率と呼ぶ．各主成分がそれぞれ元のデータのばらつき情報（分散）をどの程度説明しているか（反映しているか）の割合
- 3段目: 累積説明率と呼ぶ．分散説明率の累積値である．
例えば，上記の結果の場合，第1主成分（PC1）は2段目の情報から元のデータの約32%の情報を説明しているということであり，第1～第3主成分の3つの主成分で元のデータの情報（ばらつき）の約63%を説明（反映）しているということとなる．
1.2 どこまでの主成分を採用するか
主成分分析を使う目的は，多変量データをより少ない数の変数に集約することであり，文字通り，元のデータの「主な成分」を抽出することである．そのため，どこまでの主成分を使うかをそれぞれの目的に応じて選択する必要がある． 選択のパターンとしては主には以下のものがある．
- 単純に第1主成分のみを使う．
- 第1主成分と第2主成分を使う．図として描画する場合にはこのパターンが多い．
- 第1主成分と第2主成分と第3主成分を使う．人が頭の中でイメージできるのは精々ここまで．
- 累積説明率が6割程度になるところまでの主成分を使う．
- fa.parallel()関数を使って，主成分数を決める．
- 次に述べる主成分の解釈可能性から選択する．
fa.parallel()関数は，因子分析の際に因子数を決めるために使われることの多い関数であるが，主成分分析でも同様に使うことができる．詳しくは因子分析を説明する際に説明する．
今回のデータだと，第3主成分で累積説明率が約6割になるので，とりあえず，第3主成分までを仮採用することにしよう．
つまり，もともと9科目分あった変数が3つの主成分に集約された，ということである．
1.3 主成分の解釈
上記の方法で実際に採用された主成分について，その主成分が一体どういう意味合いを持った成分なのかの解釈が求められることは多い．
解釈をする場合には，主成分の各変数に対する主成分負荷量を見ることが重要である． 主成分負荷量とは，主成分がどの程度その変数に依存しているかを示す指標であり，主成分負荷量の絶対値が大きいほどその変数が主成分に影響を与えていることを示す．正負の向きは関係の向きを表す．
主成分負荷量は，pcaオブジェクトのrotationという変数に格納されており，pca$rotationで取得できる．
そのまま表示させると小数点以下の桁数が多くなってしまうので，以下の例では，
round()関数によって四捨五入し，小数点第2位まで表示させている．
round(pca$rotation,2)## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9
## 国語 -0.27 0.46 -0.13 -0.01 -0.44 -0.67 0.04 -0.16 0.14
## 数学IA -0.43 -0.30 -0.04 0.11 -0.02 0.02 0.55 -0.43 -0.48
## 数学IIB -0.42 -0.34 -0.17 0.13 -0.02 -0.02 0.32 0.55 0.51
## 世界史B -0.21 0.46 -0.37 0.00 -0.36 0.69 0.02 0.05 -0.04
## 地理B -0.21 0.39 -0.37 0.04 0.80 -0.11 0.02 0.01 -0.04
## 物理 -0.37 -0.29 -0.19 0.45 -0.05 -0.03 -0.72 -0.10 -0.11
## 化学 -0.30 -0.23 -0.15 -0.88 0.02 -0.01 -0.25 -0.05 -0.01
## 英語R -0.35 0.24 0.55 -0.04 0.01 -0.03 -0.09 0.54 -0.47
## 英語L -0.36 0.15 0.57 0.02 0.17 0.24 -0.06 -0.43 0.51この結果から，PC1は絶対値の大きいものから順に数学IA，数学IIB，物理が反映された主成分であり「理系得点」と解釈できる．
逆にPC2は国語，世界史，地理Bの順で大きく反映されており「文系得点」と解釈できるだろう．
第3主成分は英語R，英語Lが他よりも強く反映されており「英語」と解釈できる．
1.4 主成分得点
さて，採用する主成分が定まったあとは，それぞれのサンプルの主成分得点を得ることが次のステップとなる．主成分得点は，主成分分析の結果を収めたオブジェクトの中にxという変数名で収められており，先の例であればpca$x`で取得できる．
head()で最初の6行だけを表示してみる．
head(pca$x)## PC1 PC2 PC3 PC4 PC5 PC6
## [1,] 0.3002578 -0.4287877 0.323789360 -0.1409791 -0.46338443 -0.22365871
## [2,] 2.5509094 -1.0032696 0.008991016 -0.4938770 0.23399564 1.51781545
## [3,] 1.2418160 0.4365513 -0.366427388 0.5703851 0.58246611 -0.50787269
## [4,] -3.7140030 1.5410311 0.089394589 -1.1549117 -0.93280178 0.01969507
## [5,] -1.0291064 1.1073519 0.246354421 -0.7394713 0.04760494 0.77567274
## [6,] -1.4739507 0.6274427 -0.182148453 0.3101549 0.08064049 0.44733989
## PC7 PC8 PC9
## [1,] -0.14867455 -0.6289950 -0.32396371
## [2,] -0.34169630 -0.8400043 -0.06907116
## [3,] -0.31937709 -0.7483367 0.51226630
## [4,] -2.02032808 -1.1774320 -0.27212998
## [5,] -0.65078259 -0.2936209 -0.77337633
## [6,] -0.09165064 -0.5664564 -0.88519141さらに主成分得点の概要を把握するために，describe()関数を使ってみる．
psych::describe(pca$x)## vars n mean sd median trimmed mad min max range skew kurtosis
## PC1 1 10000 0 1.70 0.00 0.00 1.72 -5.96 6.20 12.16 0.00 -0.04
## PC2 2 10000 0 1.27 0.02 0.00 1.28 -5.20 5.28 10.48 -0.01 0.04
## PC3 3 10000 0 1.08 -0.01 0.00 1.07 -4.49 3.67 8.16 0.00 0.00
## PC4 4 10000 0 0.88 0.00 0.00 0.88 -3.16 3.25 6.41 0.01 -0.06
## PC5 5 10000 0 0.84 0.01 0.01 0.85 -3.55 3.09 6.63 -0.06 -0.05
## PC6 6 10000 0 0.76 0.00 0.00 0.75 -3.17 2.64 5.81 0.01 -0.01
## PC7 7 10000 0 0.72 0.00 0.00 0.71 -2.57 2.90 5.47 0.01 0.06
## PC8 8 10000 0 0.63 0.00 0.00 0.63 -2.40 2.36 4.75 -0.02 -0.07
## PC9 9 10000 0 0.61 0.00 0.00 0.62 -2.36 2.35 4.71 0.00 -0.05
## se
## PC1 0.02
## PC2 0.01
## PC3 0.01
## PC4 0.01
## PC5 0.01
## PC6 0.01
## PC7 0.01
## PC8 0.01
## PC9 0.01今回の主成分分析を実行する際に標準化を行っているため，主成分得点も標準化された得点として抽出されている．このため主成分得点の平均は0となっている．一方標準偏差は元データのばらつきが大きく反映されている主成分だと1を超える（つまり，元の単独の変数が持っていた情報よりも多くの情報を持っている）し，あまり反映されていないものだと1を下回る（つまり元の単独の変数が持っていた情報よりも少ない情報しか持っていない）．
1.5 主成分分析の使いみち
さて，以上が主成分分析の具体的な手法であるが，主成分分析はどのような場面で使われるのだろうか．
主成分分析の目的は，多変量データをより少ない数の変数に集約することである．そのため，主成分分析は「数多くの変数」を扱う場面で用いられる．
例えば，以下のような場面で主成分分析が使われることがある．
マクロ経済指標の分析
- 目的: 多数の経済指標（GDP成長率，失業率，インフレ率，貿易収支など）の関係性を理解し，データを簡略化する．
- 例: 世界銀行やIMFが複数の国の経済指標を主成分分析して，各国の経済状況を総合的に評価するための指数（例えば「経済健全性指標」）を作成する．
株式市場のセクター分析
- 目的: 株式市場における複数の株価指標（売上高，利益率，PER，PBRなど）を主成分分析して，主要な変動要因を特定する．
- 例: 投資ファンドが株式市場のセクター（テクノロジー，金融，エネルギーなど）のパフォーマンスに共通する要因を抽出し，リスクを分散した投資戦略を策定する．
顧客行動データの分析
- 目的: 購買履歴やアンケート結果など，多次元の顧客データから行動パターンを明らかにする．
- 例: 小売業やEコマースで，購買頻度，支出額，購入カテゴリなどを主成分分析して，顧客を類似したグループに分類し，マーケティング施策を最適化する．
もちろん，実際には主成分析だけでこれらの分析ができるわけではなく，他の分析手法（例えば回帰分析やクラスター分析）と組み合わせて使われることが多い．
特に回帰分析の際に，数多くの変数を説明変数（独立変数）として設定してしまうと，個々の説明変数の回帰係数が有意になりにくくなってしまったり，多重共線性が生じやすくなってしまう．そういう場合には，まず主成分分析を実施して，多数の説明変数を小数のものへと集約しておくと，回帰分析の結果として解釈しやすい結果が出力されるようになる．
2 因子分析
因子分析も主成分分析と同様に沢山ある測定項目をそれより少ない項目に集約する分析である．ただ，「計測値データから，データのばらつきをよく表す成分を抽出する」ことを目的とする主成分分析とは異なり，因子分析では「計測値データから，背後に潜在していて結果に影響を与えている要因（因子）を推定する」ことを目的とする（以下の図を参照）．
因子分析は組織行動論や消費者行動論，あるいはミクロ経済学など人の心理や行動を扱う様々な分野で利用されるが，いずれもアンケート結果に対して因子分析を掛けるというものである．そこでを以下ではアンケート調査と因子分析の理論的背景について解説し，そののち，因子分析の実施方法について説明する．
2.1 アンケートと潜在変数（因子）
たとえば，心理アンケートをやっていて「似たような質問項目が繰り返し出てきている」と感じたことはないだろうか．実際に心理アンケートでは意味内容が似通った（ただし文言は異なる）質問を複数並べることが一般的である．なぜそのようなことをするのだろうか．それを理解するためには次のことを理解しておく必要がある．
- 基本的に人の心理はあくまでその人の心の中に潜在的に存在しているものであり，それを物理的に直接計測することはできない．
- そのため，質問項目を示してそれへの回答（心理が反映された行動）を元にその心理を推定する必要がある（ここでは質問への回答は１～７や１～５などの数値で「程度」を回答しているものとする）
- しかし，質問項目文を構成している個々の単語の意味の受け止め方やイメージされることは人それぞれ微妙に異なっているし，その時の気分によっても印象は変わるため，どれだけアンケート作成者が熟慮して作ったとしても，質問項目がその心理を計りとれる保証はない
つまり，一つの質問項目を作成したときに，その回答には，以下の図のようにこちらが計りたいと思っている心理が反映されている部分と，人々の間で生じるわずかな言葉の意味の受け止め方の違いやその時々の気分の違いに起因してランダムに混入する誤差との，2つの部分で構成されていることになる．
質問項目が1項目しかないと，この誤差の部分が相対的に大きくなってしまうため，その項目の回答を「測りたい心理」として用いることが難しくなる．
ではどうすればよいだろうか． 文言が異なるが「似たような意味」を持つ質問項目をもう一つ持ってきて，それの回答と先の回答を足し合わせてみると，「測りたい心理」が反映されている箇所は，そのまま重ね合わされて増幅される．一方で誤差の部分は全くランダムに混入しているものになるので，サンプル単位では増幅されたり，相殺されたり，サンプルごとに異なってくるが，データ全体でみれば結局トータルでは元の誤差とそれほど誤差の大きさは変わらなくなる．
数式的に表現すれば，例えば
\[ 1つ目の項目の回答 = 0.5 \times 測りたい心理 + 0.5 \times 誤差_1\\ 2つ目の項目の回答 = 0.5 \times 測りたい心理 + 0.5 \times 誤差_2 \]
とすれば
\[ 1つ目の回答+2つ目の回答 =1.0 \times 測りたい心理 + 0.5 \times 誤差_3 \]
となる．\(誤差_3\)は\(誤差_1\)と\(誤差_2\)の和であり，サンプル単独では増幅されることもあれば相殺されることもあるが，そうした増幅と相殺が組み合わさって，データ全体では結局\(誤差_3\)は\(誤差_1\)や\(誤差_2\)と同程度の大きさになる．
つまり単独の項目では，測りたい心理と誤差とが1:1であったのが，2つを足し合わせることによって2:1とできるようになる． 質問項目が3つになれば3:1となるし，4つになれば4:1となる．
このような原理から，一見すると「似たような意味」となる質問項目文を複数作成し，それらを重ね合わせることによって，誤差の部分を相殺させるとともに，計りたい心理が反映されている部分を増幅させているのが心理アンケートの項目なのである．
そして，以上のような考えをベースとして，実際に測定した複数の項目から潜在している（直接計測ができない）因子（の値）を推定しようというのが因子分析である．因子分析の構図を図で示すと以下のようになる．
主成分分析はあくまで得られた結果からばらつき情報をうまく反映した成分を「抽出」する結果駆動型であるのに対して，因子分析ではあくまで理論として潜在的な因子を仮定し，その因子の表れが測定結果であるという前提の下に，結果から因子を逆算する形で推定する理論駆動型である，ということが両者の違いとなる．
2.2 基本的な因子分析
では実際に以下の例題を使って，因子分析を実行してみよう．
学生のSNSの利用状況と学生の性格（好奇心，社交性，計画性）との関連を調査するために，これら3つの性格特性を測定するためのアンケートとして以下に示す9つの項目のアンケートを作成した．順にQ1, Q4, Q7は「好奇心」に関する項目， Q2, Q5，Q8は「社交性」に関する項目， Q3, Q6, Q9は「計画性」に関する項目である．また(-)をつけているものは逆転項目である．
- Q1: 新しいものに興味がある
- Q2: 他人とのコミュニケーションが好きである
- Q3: 物事をじっくり考える方である
- Q4: 知らないものにはとりあえず近づいてみる
- Q5: 人の意見よりも自分の意見を重視する（－）
- Q6: 思いついたらすぐに行動する方である（－）
- Q7: 同じことを続けるのは苦手である
- Q8: 誰とでも仲良くなれる方である
- Q9: 計画通りに物事が進まないとイライラする方である
学生400名に対して，これらのアンケートに回答（1：全く当てはまらない～5：すごく当てはまる，の5段階尺度）してもらうとともに，各自のスマートフォンに記録されているSNSアプリの1日に当たりの平均利用時間を回答してもらった．その結果，このようなデータが得られた．
このデータから，性格とSNSの利用時間との関係を明らかにせよ．
このような課題の場合，手続として，各性格特性の値を因子分析によって求め，得られた因子得点を説明変数として，SNS利用時間を目的変数とした回帰分析を行うことになる．
そこで，以下で因子分析を行なっていく．
2.2.1 データの確認
因子分析を行う前に，まずはどういうデータを確認してみよう．
# データの読み込み
data <- read.csv("./practice/example_13_personality.csv")
# データの内容確認
head(data) # 冒頭6行表示## Q1 Q2 Q3 Q4 Q5 Q6 Q7 Q8 Q9 SNS
## 1 2 3 3 3 3 3 3 3 3 3.6
## 2 3 2 2 2 3 3 2 2 3 4.5
## 3 4 3 3 4 3 3 4 2 3 4.1
## 4 3 3 4 3 3 2 2 3 5 2.8
## 5 3 3 3 3 3 3 4 3 2 4.4
## 6 4 2 3 4 4 3 4 2 3 4.8summary(data) # 記述統計量の表示## Q1 Q2 Q3 Q4 Q5
## Min. :1.000 Min. :1.000 Min. :1.000 Min. :1.000 Min. :1.0
## 1st Qu.:2.000 1st Qu.:2.000 1st Qu.:2.000 1st Qu.:2.000 1st Qu.:2.0
## Median :3.000 Median :3.000 Median :3.000 Median :3.000 Median :3.0
## Mean :2.728 Mean :2.905 Mean :2.772 Mean :2.888 Mean :2.9
## 3rd Qu.:3.000 3rd Qu.:3.000 3rd Qu.:3.000 3rd Qu.:3.000 3rd Qu.:3.0
## Max. :5.000 Max. :5.000 Max. :5.000 Max. :5.000 Max. :5.0
## Q6 Q7 Q8 Q9
## Min. :1.000 Min. :1.000 Min. :1.000 Min. :1.000
## 1st Qu.:2.000 1st Qu.:2.000 1st Qu.:3.000 1st Qu.:3.000
## Median :3.000 Median :3.000 Median :3.000 Median :3.000
## Mean :2.855 Mean :2.947 Mean :3.013 Mean :3.195
## 3rd Qu.:3.000 3rd Qu.:3.000 3rd Qu.:3.000 3rd Qu.:4.000
## Max. :5.000 Max. :5.000 Max. :5.000 Max. :5.000
## SNS
## Min. :1.200
## 1st Qu.:3.275
## Median :3.800
## Mean :3.811
## 3rd Qu.:4.500
## Max. :6.200round(cor(data),2) # 相関行列の表示.round関数を使って2桁までにしている## Q1 Q2 Q3 Q4 Q5 Q6 Q7 Q8 Q9 SNS
## Q1 1.00 0.28 -0.01 0.56 -0.21 0.03 0.45 0.15 -0.04 0.26
## Q2 0.28 1.00 0.24 0.23 -0.59 -0.17 0.21 0.46 0.15 0.26
## Q3 -0.01 0.24 1.00 0.07 -0.12 -0.69 0.09 0.07 0.52 -0.14
## Q4 0.56 0.23 0.07 1.00 -0.20 -0.03 0.68 0.19 0.06 0.22
## Q5 -0.21 -0.59 -0.12 -0.20 1.00 0.03 -0.11 -0.68 -0.01 -0.28
## Q6 0.03 -0.17 -0.69 -0.03 0.03 1.00 -0.03 -0.01 -0.70 0.19
## Q7 0.45 0.21 0.09 0.68 -0.11 -0.03 1.00 0.13 0.02 0.26
## Q8 0.15 0.46 0.07 0.19 -0.68 -0.01 0.13 1.00 0.00 0.29
## Q9 -0.04 0.15 0.52 0.06 -0.01 -0.70 0.02 0.00 1.00 -0.24
## SNS 0.26 0.26 -0.14 0.22 -0.28 0.19 0.26 0.29 -0.24 1.00相関を見ると基本的に想定していた通りの３項目ずつで相関値が心理学的に高いとされる値になっていることや，逆転項目になっている項目については負の値になっていることがわかる．
加えて，一部の項目では0.2〜0.3ほどの相関（心理学的には弱いながらも相関があるとされる）も見られる．
2.2.2 因子分析の実行
Rで因子分析を行う関数はいくつかあるが，もっとも使いやすいのはpsychパッケージに含まれるfa()関数である．引数は以下の通り．
- 第1引数: 因子分析に掛ける多変量が含まれたデータフレーム
- 第2引数: 因子数． 何個の因子を推定するのかを指定する．
- 第3引数 rotate: 回転方法（省略可）． 因子をどのように回転させるかを指定する．
- 第4引数 fm: 推定方法（省略可）． 因子の推定方法を指定する．
第3引数や第4引数については後で説明する．デフォルト値が設定されているので，基本的には省略してもよい．
結果は一旦オブジェクトに収めた上で，print()で出力させる．その際には以下の引数とする．
- 第1引数: faの出力結果を収めたオブジェクト
- 第2引数 sort（省略可）: TRUEにすると因子負荷量の絶対値の大きい順に表示される．省略するとデータ上の項目の並びの順に出力される．TRUEにするのが一般的．．
- 第3引数 digits（省略可）: 表記の小数点以下の桁数を指定する．省略すると2桁表記となる．
#psychパッケージの読み込み
library(psych)##
## Attaching package: 'psych'## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha# 因子分析の実行
fa_result <- fa(data[c(1:9)], 3) #dataには最後の列にSNSの利用時間が含まれているので，それを除いたものを指定する
print(fa_result, sort=TRUE) #digits=3とすると結果が3桁表記になる## Factor Analysis using method = minres
## Call: fa(r = data[c(1:9)], nfactors = 3)
## Standardized loadings (pattern matrix) based upon correlation matrix
## item MR2 MR1 MR3 h2 u2 com
## Q6 6 -0.96 0.02 0.02 0.91 0.092 1.0
## Q9 9 0.74 -0.03 0.01 0.54 0.459 1.0
## Q3 3 0.72 0.09 0.02 0.54 0.462 1.0
## Q5 5 0.02 -0.95 0.03 0.89 0.114 1.0
## Q8 8 -0.03 0.72 0.01 0.52 0.481 1.0
## Q2 2 0.16 0.59 0.14 0.46 0.545 1.3
## Q4 4 0.01 -0.01 0.90 0.81 0.187 1.0
## Q7 7 0.01 -0.04 0.76 0.56 0.438 1.0
## Q1 1 -0.07 0.10 0.60 0.40 0.603 1.1
##
## MR2 MR1 MR3
## SS loadings 2.01 1.82 1.79
## Proportion Var 0.22 0.20 0.20
## Cumulative Var 0.22 0.43 0.62
## Proportion Explained 0.36 0.32 0.32
## Cumulative Proportion 0.36 0.68 1.00
##
## With factor correlations of
## MR2 MR1 MR3
## MR2 1.00 0.08 0.06
## MR1 0.08 1.00 0.26
## MR3 0.06 0.26 1.00
##
## Mean item complexity = 1
## Test of the hypothesis that 3 factors are sufficient.
##
## df null model = 36 with the objective function = 3.63 with Chi Square = 1434.4
## df of the model are 12 and the objective function was 0.06
##
## The root mean square of the residuals (RMSR) is 0.02
## The df corrected root mean square of the residuals is 0.03
##
## The harmonic n.obs is 400 with the empirical chi square 7.73 with prob < 0.81
## The total n.obs was 400 with Likelihood Chi Square = 22.27 with prob < 0.035
##
## Tucker Lewis Index of factoring reliability = 0.978
## RMSEA index = 0.046 and the 90 % confidence intervals are 0.012 0.076
## BIC = -49.63
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy
## MR2 MR1 MR3
## Correlation of (regression) scores with factors 0.96 0.95 0.93
## Multiple R square of scores with factors 0.93 0.91 0.87
## Minimum correlation of possible factor scores 0.85 0.81 0.732.2.3 結果の見方
様々な結果が表記されるがとりあえず確認すべき箇所は以下の箇所である．
## item MR2 MR1 MR3
## Q6 6 -0.96 0.02 0.02
## Q9 9 0.74 -0.03 0.01
## Q3 3 0.72 0.09 0.02
## Q5 5 0.02 -0.95 0.03
## Q8 8 -0.03 0.72 0.01
## Q2 2 0.16 0.59 0.14
## Q4 4 0.01 -0.01 0.90
## Q7 7 0.01 -0.04 0.76
## Q1 1 -0.07 0.10 0.60この部分は推定された3つの因子が各質問項目にどの程度結びついているか（因子を説明変数，質問項目を目的変数としたときの回帰係数）を表している．この値のことを因子負荷量と呼ぶ．因子負荷量は−1から1の値を取り，絶対値が1に近いほど強い結びつきを，0に近いほど弱い結びつきを意味する．正負は関係性の向きを表している
因子負荷量は絶対値が0.4を超えている場合に，その項目はその因子に結びついていると判断するのが一般的である． 今回の結果だと，MR2はQ3, Q6, Q9と，MR1はQ2, Q5, Q8と，MR3はQ1, Q4, Q7とそれぞれ結びついている，と判断される． つまりこの結果から，当初想定していた通りの因子構造が推定された，ということがわかる．
2.2.4 結果の図示
因子負荷量を図示することで，因子と項目の関係をより直感的に理解することができる．
psychパッケージにはfa()関数の結果を描画してくれる便利な関数として，fa.diagram()関数が用意されている．
引数は以下の通り．
- 第1引数: faの出力結果を収めたオブジェクト
- digits: 表記の小数点以下の桁数を指定する．省略すると2桁表記となる．
- cut: 因子負荷量の絶対値がこの値を超える場合に線を引く．省略すると0.3が設定される．
- simple: TRUEにすると各項目につき因子負荷量が最も高い因子との結びつきのみが表示される．FALSEにすると全ての因子との結びつきが表示する．ただしいずれもcutを下回るものについては表示されない．
図の中で赤い破線の矢印は因子負荷量が負の値となっているものを示している．
# 図示
fa.diagram(fa_result, digits=2, cut=0.4, simple = TRUE)
試しに，cutを0.05まで下げて，simpleをFALSEにしてみると，上では表示されたなかった結びつきや，因子同士の相関が表示される．
# SimpleとCutをいじる
fa.diagram(fa_result, digits=2, cut=0.05, simple = FALSE)
2.2.5 因子得点
推定された因子得点はfa()の出力結果のscoresに格納されている． 以下では，推定された因子得点の冒頭6行を表示している．
head(fa_result$scores)## MR2 MR1 MR3
## [1,] -0.1062418 -0.11977804 0.005913257
## [2,] -0.3743964 -0.45588883 -1.088542941
## [3,] -0.1291166 -0.21248013 1.512048370
## [4,] 1.2922966 -0.07725064 -0.140764867
## [5,] -0.2703695 -0.05647853 0.493830786
## [6,] -0.1774206 -1.43937200 1.440993591因子得点は平均が0，標準偏差が1になるように標準化されている． （実際には，計算誤差等により標準偏差が完全に1となるわけではないが1に近い値になっている）
次のようcbind()関数を用いると，元のデータに因子得点が格納される．
data <- cbind(data, fa_result$scores)
head(data)## Q1 Q2 Q3 Q4 Q5 Q6 Q7 Q8 Q9 SNS MR2 MR1 MR3
## 1 2 3 3 3 3 3 3 3 3 3.6 -0.1062418 -0.11977804 0.005913257
## 2 3 2 2 2 3 3 2 2 3 4.5 -0.3743964 -0.45588883 -1.088542941
## 3 4 3 3 4 3 3 4 2 3 4.1 -0.1291166 -0.21248013 1.512048370
## 4 3 3 4 3 3 2 2 3 5 2.8 1.2922966 -0.07725064 -0.140764867
## 5 3 3 3 3 3 3 4 3 2 4.4 -0.2703695 -0.05647853 0.493830786
## 6 4 2 3 4 4 3 4 2 3 4.8 -0.1774206 -1.43937200 1.4409935912.2.6 参考：例題の遂行
因子分析としては以上となるが，さらに例題での課題となっている学生のSNS利用時間との関連を調べてみよう．
まずはこれらの因子得点をdataに格納した上で，SNS利用時間との相関を見てみよう．
cor(data[c("SNS","MR1","MR2","MR3")])## SNS MR1 MR2 MR3
## SNS 1.0000000 0.30875867 -0.19158604 0.27155917
## MR1 0.3087587 1.00000000 0.08352857 0.28194345
## MR2 -0.1915860 0.08352857 1.00000000 0.06512139
## MR3 0.2715592 0.28194345 0.06512139 1.00000000SNSの列をみると，相関係数は0.2〜0.3程度であり，心理学的には弱いながらも相関があると言える．
さらに，「社交性（MR1）や好奇心（MR3）が高く，計画性（MR2）が低い学生ほどSNSを多く利用する」という仮説を検証するために，標準化重回帰分析を行ってみよう．
#標準化
data[c("SNS","MR1","MR2","MR3")] <- scale(data[c("SNS","MR1","MR2","MR3")])
lm_result <- lm(SNS ~ MR1+MR2+MR3, data)
summary(lm_result)##
## Call:
## lm(formula = SNS ~ MR1 + MR2 + MR3, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.79377 -0.56073 0.03097 0.58107 2.53823
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.392e-16 4.535e-02 0.000 1
## MR1 2.684e-01 4.743e-02 5.658 2.94e-08 ***
## MR2 -2.277e-01 4.560e-02 -4.994 8.91e-07 ***
## MR3 2.107e-01 4.737e-02 4.449 1.12e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9069 on 396 degrees of freedom
## Multiple R-squared: 0.1837, Adjusted R-squared: 0.1775
## F-statistic: 29.71 on 3 and 396 DF, p-value: < 2.2e-16この結果から，いずれも有意な回帰係数が得られており，「社交性や好奇心が高く，計画性が低い学生ほどSNSを多く利用する」という仮説が支持されたことがわかる．
因子得点か尺度化か
実際に心理学や行動科学の論文では，実は因子得点を使うことはあまりなく，多くの場合は因子分析の結果として各因子と項目の関係が見出された後は，関係のある項目同士の単純な平均値や合計値を各因子の「尺度」とすることが多い．
例えば，上記の例だと，以下のようにする．なお，Q5とQ6は逆転項目であるため，数値をひっくり返す必要があるため，6-をつけている．6から引いているのは5段階尺度のためであり，7段階尺度の場合には8から引くことになる．
# 因子得点を尺度化する
社交性<- (data$Q2 + 6-data$Q5 + data$Q8)/3
好奇心<- (data$Q1 + data$Q4 + data$Q7)/3
計画性<- (data$Q3 + 6-data$Q6 + data$Q9)/3実際にこうして得られた尺度得点を用いて，SNS利用時間との関係を見てみよう．
#作成した尺度をデータフレームに格納
data <- cbind(data, 社交性, 好奇心, 計画性)
#相関
cor(data[c("SNS","社交性","好奇心","計画性")])## SNS 社交性 好奇心 計画性
## SNS 1.0000000 0.3264385 0.29234403 -0.21677103
## 社交性 0.3264385 1.0000000 0.26870709 0.12575885
## 好奇心 0.2923440 0.2687071 1.00000000 0.03368339
## 計画性 -0.2167710 0.1257588 0.03368339 1.00000000#回帰分析
data[c("SNS","社交性","好奇心","計画性")] <- scale(data[c("SNS","社交性","好奇心","計画性")])
lm_result <- lm(SNS ~ 社交性+好奇心+計画性, data)
summary(lm_result)##
## Call:
## lm(formula = SNS ~ 社交性 + 好奇心 + 計画性, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.49443 -0.57877 0.05327 0.59414 2.73793
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.098e-16 4.435e-02 0.000 1
## 社交性 3.001e-01 4.644e-02 6.462 3.03e-10 ***
## 好奇心 2.205e-01 4.610e-02 4.784 2.43e-06 ***
## 計画性 -2.619e-01 4.476e-02 -5.852 1.02e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.887 on 396 degrees of freedom
## Multiple R-squared: 0.2192, Adjusted R-squared: 0.2133
## F-statistic: 37.06 on 3 and 396 DF, p-value: < 2.2e-16先ほどの因子得点の場合と似たような結果となっている．
では，因子得点と尺度得点のいずれを使うかは，研究目的によって異なるが，一般に尺度得点を用いることが多いのは，因子得点はあくまでその時のデータでのみ有効な算出方法によって因子得点を算出しているのであり，他のデータを用いて算出した因子得点との比較が難しいためである．
一方，尺度得点は，得点の算出方法が明確であり，他のデータでも同じ方法で算出できるため，他の研究との比較が容易である．
したがって，自身の研究結果と他の研究結果とを比較したいような場合には，尺度得点を用いることが一般的である．
ただし尺度得点を用いる場合には，因子分析だけでなく，それらを合算することの妥当性（内的整合性：つまり互いに十分な相関を持っているかどうか）をそれぞれの尺度得点を求める前に，それぞれごとに評価する必要もある．一般にはクロンバックの\(\alpha\)係数や，\(\omega\)係数を用いて評価することが多い．以下にそれぞれ求めている例を紹介する．
alpha()関数，omega()関数はいずれもpsychパッケージに含まれている．いずれも引数は以下の通り．
- 第1引数: 分析対象となるデータの列を指定する．
- check.keys: TRUEにすると，逆転項目と思われるものが含まれている場合に，自動で値を逆転させて評価する（ただし，必ずしも本当に逆転項目であるとは限らないため，注意が必要）．
#クロンバックのα係数
alpha(data[c("Q1", "Q4", "Q7")], check.keys = TRUE)##
## Reliability analysis
## Call: alpha(x = data[c("Q1", "Q4", "Q7")], check.keys = TRUE)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.79 0.79 0.74 0.56 3.9 0.018 2.9 0.64 0.56
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.76 0.79 0.83
## Duhachek 0.76 0.79 0.83
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## Q1 0.81 0.81 0.68 0.68 4.3 0.019 NA 0.68
## Q4 0.62 0.62 0.45 0.45 1.6 0.038 NA 0.45
## Q7 0.72 0.72 0.56 0.56 2.5 0.028 NA 0.56
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## Q1 400 0.79 0.79 0.61 0.55 2.7 0.74
## Q4 400 0.89 0.89 0.82 0.73 2.9 0.77
## Q7 400 0.85 0.84 0.74 0.64 2.9 0.77
##
## Non missing response frequency for each item
## 1 2 3 4 5 miss
## Q1 0.03 0.36 0.48 0.13 0.00 0
## Q4 0.03 0.25 0.53 0.17 0.01 0
## Q7 0.03 0.23 0.53 0.20 0.02 0alpha(data[c("Q2", "Q5", "Q8")], check.keys = TRUE)## Warning in alpha(data[c("Q2", "Q5", "Q8")], check.keys = TRUE): Some items were negatively correlated with the first principal component and were automatically reversed.
## This is indicated by a negative sign for the variable name.##
## Reliability analysis
## Call: alpha(x = data[c("Q2", "Q5", "Q8")], check.keys = TRUE)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.8 0.8 0.75 0.58 4.1 0.018 3 0.62 0.59
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.76 0.8 0.83
## Duhachek 0.76 0.8 0.83
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## Q2 0.81 0.81 0.68 0.68 4.3 0.019 NA 0.68
## Q5- 0.62 0.63 0.46 0.46 1.7 0.037 NA 0.46
## Q8 0.74 0.74 0.59 0.59 2.9 0.026 NA 0.59
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## Q2 400 0.82 0.80 0.63 0.57 2.9 0.80
## Q5- 400 0.89 0.89 0.83 0.74 3.1 0.71
## Q8 400 0.83 0.84 0.73 0.63 3.0 0.69
##
## Non missing response frequency for each item
## 1 2 3 4 5 miss
## Q2 0.04 0.25 0.51 0.19 0.01 0
## Q5 0.01 0.26 0.54 0.17 0.01 0
## Q8 0.01 0.20 0.58 0.20 0.01 0alpha(data[c("Q3", "Q6", "Q9")], check.keys = TRUE)## Warning in alpha(data[c("Q3", "Q6", "Q9")], check.keys = TRUE): Some items were negatively correlated with the first principal component and were automatically reversed.
## This is indicated by a negative sign for the variable name.##
## Reliability analysis
## Call: alpha(x = data[c("Q3", "Q6", "Q9")], check.keys = TRUE)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.84 0.84 0.8 0.64 5.3 0.014 3 0.68 0.69
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.81 0.84 0.87
## Duhachek 0.81 0.84 0.87
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## Q3 0.83 0.83 0.70 0.70 4.8 0.017 NA 0.70
## Q6- 0.69 0.69 0.52 0.52 2.2 0.031 NA 0.52
## Q9 0.82 0.82 0.69 0.69 4.5 0.018 NA 0.69
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## Q3 400 0.85 0.85 0.72 0.66 2.8 0.79
## Q6- 400 0.92 0.92 0.87 0.80 3.1 0.77
## Q9 400 0.85 0.85 0.74 0.67 3.2 0.79
##
## Non missing response frequency for each item
## 1 2 3 4 5 miss
## Q3 0.04 0.32 0.48 0.15 0.01 0
## Q6 0.03 0.28 0.50 0.18 0.01 0
## Q9 0.00 0.18 0.47 0.31 0.04 0#ω係数
omega(data[c("Q1", "Q4", "Q7")], check.keys = TRUE)## Warning in cov2cor(t(w) %*% r %*% w): diag(V) had non-positive or NA entries;
## the non-finite result may be dubious
## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar, check.keys = TRUE)
## Alpha: 0.79
## G.6: 0.74
## Omega Hierarchical: 0.02
## Omega H asymptotic: 0.02
## Omega Total 0.81
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* F2* F3* h2 h2 u2 p2 com
## Q1 0.61 0.40 0.40 0.60 0.03 1.18
## Q4 0.89 0.81 0.81 0.19 0.02 1.04
## Q7 0.76 0.60 0.60 0.40 0.01 1.09
##
## With Sums of squares of:
## g F1* F2* F3* h2
## 0.03 1.74 0.04 0.00 1.18
##
## general/max 0.02 max/min = Inf
## mean percent general = 0.02 with sd = 0.01 and cv of 0.5
## Explained Common Variance of the general factor = 0.02
##
## The degrees of freedom are -3 and the fit is 0
## The number of observations was 400 with Chi Square = 0 with prob < NA
## The root mean square of the residuals is 0
## The df corrected root mean square of the residuals is NA
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 0 and the fit is 0.98
## The number of observations was 400 with Chi Square = 386.7 with prob < NA
## The root mean square of the residuals is 0.56
## The df corrected root mean square of the residuals is NA
##
## Measures of factor score adequacy
## g F1* F2* F3*
## Correlation of scores with factors 0.13 0.92 0.28 0
## Multiple R square of scores with factors 0.02 0.85 0.08 0
## Minimum correlation of factor score estimates -0.97 0.70 -0.84 -1
##
## Total, General and Subset omega for each subset
## g F1* F2* F3*
## Omega total for total scores and subscales 0.81 0.81 NA NA
## Omega general for total scores and subscales 0.02 0.02 NA NA
## Omega group for total scores and subscales 0.80 0.80 NA NAomega(data[c("Q2", "Q5", "Q8")], check.keys = TRUE)## Warning in cov2cor(t(w) %*% r %*% w): diag(V) had non-positive or NA entries;
## the non-finite result may be dubious
## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar, check.keys = TRUE)
## Alpha: 0.8
## G.6: 0.75
## Omega Hierarchical: 0
## Omega H asymptotic: 0
## Omega Total 0.82
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* F2* F3* h2 h2 u2 p2 com
## Q2 0.63 0.43 0.43 0.57 0.01 1.14
## Q5- 0.91 0.84 0.84 0.16 0.00 1.01
## Q8 0.75 0.58 0.58 0.42 0.00 1.08
##
## With Sums of squares of:
## g F1* F2* F3* h2
## 0.01 1.80 0.05 0.00 1.23
##
## general/max 0 max/min = Inf
## mean percent general = 0 with sd = 0.01 and cv of 1.73
## Explained Common Variance of the general factor = 0
##
## The degrees of freedom are -3 and the fit is 0
## The number of observations was 400 with Chi Square = 0 with prob < NA
## The root mean square of the residuals is 0
## The df corrected root mean square of the residuals is NA
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 0 and the fit is 1.05
## The number of observations was 400 with Chi Square = 417.32 with prob < NA
## The root mean square of the residuals is 0.58
## The df corrected root mean square of the residuals is NA
##
## Measures of factor score adequacy
## g F1* F2* F3*
## Correlation of scores with factors 0.07 0.94 0.30 0
## Multiple R square of scores with factors 0.00 0.88 0.09 0
## Minimum correlation of factor score estimates -0.99 0.75 -0.82 -1
##
## Total, General and Subset omega for each subset
## g F1* F2* F3*
## Omega total for total scores and subscales 0.82 0.82 NA NA
## Omega general for total scores and subscales 0.00 0.00 NA NA
## Omega group for total scores and subscales 0.82 0.82 NA NAomega(data[c("Q3", "Q6", "Q9")], check.keys = TRUE)## Warning in cov2cor(t(w) %*% r %*% w): diag(V) had non-positive or NA entries;
## the non-finite result may be dubious
## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar, check.keys = TRUE)
## Alpha: 0.84
## G.6: 0.8
## Omega Hierarchical: 0.83
## Omega H asymptotic: 0.96
## Omega Total 0.86
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* F2* F3* h2 h2 u2 p2 com
## Q3 0.72 0.20 0.55 0.55 0.45 0.93 1.15
## Q6- 0.93 0.90 0.90 0.10 0.96 1.08
## Q9 0.73 0.20 0.57 0.57 0.43 0.93 1.15
##
## With Sums of squares of:
## g F1* F2* F3* h2
## 1.91 0.06 0.06 0.00 1.44
##
## general/max 1.33 max/min = Inf
## mean percent general = 0.94 with sd = 0.02 and cv of 0.02
## Explained Common Variance of the general factor = 0.94
##
## The degrees of freedom are -3 and the fit is 0
## The number of observations was 400 with Chi Square = 0 with prob < NA
## The root mean square of the residuals is 0
## The df corrected root mean square of the residuals is NA
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 0 and the fit is 0.01
## The number of observations was 400 with Chi Square = 3.55 with prob < NA
## The root mean square of the residuals is 0.02
## The df corrected root mean square of the residuals is NA
##
## Measures of factor score adequacy
## g F1* F2* F3*
## Correlation of scores with factors 0.94 0.24 0.24 0
## Multiple R square of scores with factors 0.89 0.06 0.06 0
## Minimum correlation of factor score estimates 0.77 -0.88 -0.88 -1
##
## Total, General and Subset omega for each subset
## g F1* F2* F3*
## Omega total for total scores and subscales 0.86 0.57 0.83 NA
## Omega general for total scores and subscales 0.83 0.53 0.80 NA
## Omega group for total scores and subscales 0.02 0.04 0.03 NAいずれもさまざまな結果が表記されるが，alpha()関数の場合にはraw_alphaの値を，omega()関数の場合にはOmega Totalの値を見るとよい．一般的にはこれらが0.7を上回っていれば合算や平均にすることに十分な妥当性があると評価される．一方，0.6〜0.7だと妥当性が十分とはいえず注意しながら使うことになる．0.6を下回ってしまった場合には妥当性でないと判断され，使うべきではない，とされる．
2.3 探索的因子分析と検証的因子分析
ところで，心理アンケートは上記に示したロジックから，アンケート項目を多く作れば，それだけ誤差は相殺し合う一方で，測りたい心理が反映されている部分は重ね合わさって増幅していく．こうしたことから，似たような意味内容の質問項目を数多く設ければ，測定したい心理をより純粋に測定することができるようになると考えられる．
しかしながら，実際にそうしてみると，もともとは同じ意味内容だと思って設けた（＝ある共通した1つの因子が反映されるだろうと思っていた）質問項目であっても，微妙な違いがそこに存在していて，その違いが増幅されることによって，複数の因子が潜在していたことが発見される，といったこともしばしばある．
そこで自分が作成して実施した質問項目から，いくつの因子を抽出できるかを探索的に見つけ出すために因子分析を行うこともある．こうした因子分析を探索的因子分析（Exploratory Factor Analysis: EFA）と呼ぶ．通常，因子分析という場合には，この探索的因子分析のことを指す．
それに対して，「これらの項目への回答はこの因子から生じるものである」という仮説を自分の中でより強く持っていて，その仮説を検証するために行われる因子分析を検証的因子分析もしくは確認的因子分析（Confirmatory Factor Analysis: CFA）と呼ぶ．特に既存の研究で十分に検討されたうえで作りこまれたアンケートを自身の研究に利用するような場合には，探索的に因子を推定するというより，既存研究で言われている通りに因子が推定できるかの検証が重要となるためCFAが用いられることが多い．
CFAを行うには，lavaanパッケージという別のパッケージ（共分散構造分析に用いるパッケージ）を使うことが一般的であるが，そちらはまた別の機会に紹介することとして，ここではEFAのみを扱うこととする．
2.4 探索的因子分析の実施方法
探索的因子分析を実施する際には，以下の手順で行うことが一般的である．
- 推定する因子数の模索
- 因子分析の実施と解釈
- 因子負荷量・適合度指標の確認
これらを順に次の例を使って説明していく．
次のデータは実際に著者が自身の研究の中で実施した病院看護師に対するアンケートの回答データである．このアンケートでは，看護師らに対して自分の上司（つまりは看護師長）のリーダーシップを評価してもらうことを目的に22の項目を設けた．Q1～7までと，19，20はどちらかと言うと「良い上司のリーダーシップ」が反映されるだろうと思って設けたのに対して，Q8～18と21，22は「嫌な上司のリーダーシップ」が反映されるだろうと思って設けた項目である．これら22個の質問項目から得られた結果から，いくつの因子が潜在しているかを調べるとともに，推定された因子の内容を解釈せよ．
具体的な質問項目：
- Q1 上司は私が落ち込んでる時にはそのことに気づくことができる
- Q2 上司は自職場だけでなく病院全体や地域全体に貢献することの重要性を強調する
- Q3 上司は何か悪いことが起ころうとしていると，そのことに気づき，伝えてくれる
- Q4 上司は，私が困難な状況に対応する際には，「私自身が最良と感じる方法で対応すれば良い」と言ってくれる
- Q5 上司は私のキャリアの発達を優先してくれている
- Q6 上司 は上司自身の利益よりも私自身の利益を優先する
- Q7 上司はいつも誠実である
- Q8 上司は私に対して指示に完全に従うように言ってくる
- Q9 上司はそれが重要なものかどうかにかかわらず，組織におけるあらゆることを自分で決めようとする
- Q10 上司はいつも会議で最後に一言発言をする
- Q11 上司はいつも従業員の前では威圧的なふるまいをする
- Q12 私は上司とともに仕事をするときにはプレッシャーを感じる
- Q13 上司は部下を厳しく指導している
- Q14 上司はタスクが達成できないと厳しく叱責してくる
- Q15 上司は「我々の部署は組織の全部署のなかで，一番の業績を上げなければならない」と強調してくる
- Q16 我々はモノゴトを進める際には上司のルールに従わなければならないさもないと厳しく罰せられる
- Q17 上司は私の（私たちの）業務の実態を分かっていないと感じる
- Q18 上司は私の（私たちの）業務を数字でしか捉えていないと感じる
- Q19 上司は私の（私たちの）業務を遂行している現場によく顔を出す
- Q20 上司は私の（私たちの）業務の状況をよく聞きにくる
- Q21 上司は私の（私たちの）ことよりも，上役の方を見ていると感じる
- Q22 上司に私の（私たちの）業務実態について話ができる機会・場が設けられていない
まずはpsychパッケージを読み込みと，データの読み込みをしておく．
library(psych)
data <- read.csv("practice/example_13_EFA.csv")また読み込んだデータの内容も一応確認してみよう．
なお[-1]をつけているのは，今回のデータの場合，冒頭の列はサンプルNo.であるためである．このように-列番号とすることで指定した列を除外することができる（参考）．
## SampleNo. Q1 Q2 Q3 Q4 Q5 Q6 Q7 Q8 Q9 Q10 Q11 Q12 Q13 Q14 Q15 Q16 Q17 Q18 Q19
## 1 1 4 4 4 4 4 4 4 2 2 3 2 2 1 2 1 1 1 1 4
## 2 2 2 2 4 3 3 3 4 2 2 4 1 2 2 2 2 4 2 2 4
## 3 3 5 2 5 5 4 4 4 2 1 1 1 1 4 2 1 1 2 1 5
## 4 4 2 4 3 4 4 3 4 2 3 3 1 3 2 1 1 1 2 2 4
## 5 5 3 3 3 4 4 3 5 2 2 3 1 1 1 1 1 3 3 3 3
## 6 6 1 1 1 1 4 2 2 1 1 2 1 1 1 1 1 1 2 2 2
## Q20 Q21 Q22
## 1 4 2 2
## 2 3 2 2
## 3 5 2 2
## 4 4 2 2
## 5 3 3 3
## 6 2 3 1## vars n mean sd median trimmed mad min max range skew kurtosis se
## Q1 1 986 3.00 1.03 3.0 3.08 1.48 1 5 4 -0.32 -0.68 0.03
## Q2 2 986 3.37 0.89 3.0 3.42 1.48 1 5 4 -0.52 0.06 0.03
## Q3 3 986 3.47 0.90 4.0 3.56 0.00 1 5 4 -0.89 0.55 0.03
## Q4 4 986 3.15 0.91 3.0 3.21 1.48 1 5 4 -0.43 -0.13 0.03
## Q5 5 986 3.34 0.88 3.0 3.41 1.48 1 5 4 -0.59 0.31 0.03
## Q6 6 986 2.97 0.92 3.0 3.03 1.48 1 5 4 -0.36 -0.06 0.03
## Q7 7 986 3.48 0.97 4.0 3.52 1.48 1 5 4 -0.54 -0.04 0.03
## Q8 8 986 2.57 1.02 2.5 2.54 0.74 1 5 4 0.29 -0.47 0.03
## Q9 9 986 2.58 0.96 2.5 2.56 0.74 1 5 4 0.33 -0.29 0.03
## Q10 10 986 3.18 0.99 3.0 3.20 1.48 1 5 4 -0.27 -0.25 0.03
## Q11 11 986 2.22 1.06 2.0 2.12 1.48 1 5 4 0.59 -0.39 0.03
## Q12 12 986 2.75 1.13 3.0 2.73 1.48 1 5 4 0.17 -0.85 0.04
## Q13 13 986 2.60 1.01 3.0 2.59 1.48 1 5 4 0.22 -0.54 0.03
## Q14 14 986 2.25 0.98 2.0 2.16 1.48 1 5 4 0.56 -0.07 0.03
## Q15 15 986 2.09 0.95 2.0 2.00 1.48 1 5 4 0.67 0.08 0.03
## Q16 16 986 2.05 1.01 2.0 1.93 1.48 1 5 4 0.70 -0.28 0.03
## Q17 17 986 2.81 1.20 3.0 2.79 1.48 1 5 4 0.13 -0.97 0.04
## Q18 18 986 2.68 1.20 3.0 2.62 1.48 1 5 4 0.29 -0.86 0.04
## Q19 19 986 3.09 1.15 3.0 3.12 1.48 1 5 4 -0.23 -0.82 0.04
## Q20 20 986 3.11 1.06 3.0 3.16 1.48 1 5 4 -0.31 -0.61 0.03
## Q21 21 986 2.75 1.08 3.0 2.72 1.48 1 5 4 0.25 -0.43 0.03
## Q22 22 986 2.66 1.02 3.0 2.64 1.48 1 5 4 0.36 -0.28 0.032.4.1 因子数の模索
質問紙をつくる際には当然ながら「推定したい因子」があるはずであり，その因子が反映されるであろう質問項目が作られているので，「推定したい因子」の数がここでいう因子数となるのが基本である．
しかしながら，実際には自分が思い描いた通りの因子構成にならず，想定よりも少ない因子構成になったり，逆に想定よりも多い因子構成となったりすることも多い．
そこで分析にかけるにあたって，まずはfa.parallel()関数や，vss()関数を使って，どの程度の因子数となるのかの目処を得ようとすることは多い．
実際に先のデータに対してこれらの関数を実行してみよう．
fa.parallel()
fa.parallel()関数の引数は以下の通りである．
- 第1引数: データフレーム
- 第2引数
fa（省略可）:"fa"，"pc","both"のいずれかを入力する．これは出力結果を因子分析用(fa)のみ表記するのか，主成分分析用(pc)のみ表記するのか，両方の結果を表記させるのかをコントロールするものである． 省略するとデフォルトは"both"となる．
fa.parallel(data[-1], fa="fa")
## Parallel analysis suggests that the number of factors = 5 and the number of components = NAこのようにfa.parallel()関数を実行すると，スクリープロットと呼ばれる図が表示されるとともに，メッセージとして因子数の候補が出力される．
図には\(\triangle\)印の青い折れ線グラフと赤点線，赤破線のグラフが表示されている（赤点線はほぼ赤破線に重なるのでわかりにくくなっている）．このグラフが表記しているものを理解するには高度な数学の知識が必要となるが，見方としては，青い折れ線グラフと赤いグラフとが交差するあたりが因子数候補となる．
実際にメッセージで表示されているものは，このようにして求めた因子数候補である．
ちなみに，こちらで述べたようにfa.parallel()関数は主成分分析の時の採用する主成分数を決めるときにも使える．その場合には，fa=“pc”と設定すればよい．
vss()
vss()関数は簡易的に因子分析を様々な因子数で行い，その結果として得られる各種適合度指標をグラフや表で出力するものである．
vss()関数の引数は以下の通りである．
- 第1引数: データフレーム
- 第2引数
n(省略可):シミュレートするときの最大因子数．省略時のデフォルトは8．
vss(data[-1])
##
## Very Simple Structure
## Call: vss(x = data[-1])
## VSS complexity 1 achieves a maximimum of 0.83 with 1 factors
## VSS complexity 2 achieves a maximimum of 0.91 with 2 factors
##
## The Velicer MAP achieves a minimum of 0.02 with 2 factors
## BIC achieves a minimum of -391.09 with 6 factors
## Sample Size adjusted BIC achieves a minimum of -85.29 with 8 factors
##
## Statistics by number of factors
## vss1 vss2 map dof chisq prob sqresid fit RMSEA BIC SABIC complex
## 1 0.83 0.00 0.035 209 4146 0.0e+00 14.8 0.83 0.138 2705 3369 1.0
## 2 0.64 0.91 0.018 188 2180 0.0e+00 8.2 0.91 0.104 884 1481 1.3
## 3 0.52 0.82 0.019 168 1465 1.1e-205 6.4 0.93 0.088 307 841 1.6
## 4 0.51 0.79 0.019 149 844 1.4e-97 5.5 0.94 0.069 -183 290 1.6
## 5 0.45 0.73 0.022 131 605 5.2e-62 4.8 0.95 0.061 -298 118 1.7
## 6 0.45 0.71 0.027 114 395 1.3e-32 4.3 0.95 0.050 -391 -29 1.8
## 7 0.45 0.71 0.031 98 315 1.3e-24 4.0 0.95 0.047 -361 -49 1.9
## 8 0.43 0.69 0.036 83 223 8.5e-15 3.7 0.96 0.041 -349 -85 2.0
## eChisq SRMR eCRMS eBIC
## 1 5421 0.109 0.115 3980
## 2 1469 0.057 0.063 173
## 3 737 0.040 0.047 -421
## 4 396 0.029 0.037 -631
## 5 249 0.023 0.031 -654
## 6 141 0.018 0.025 -645
## 7 102 0.015 0.023 -573
## 8 71 0.012 0.021 -501グラフの見方
まずはグラフの見方であるが，これはVSS（Very Simple Structure）と呼ばれる方法で因子数を決定する際に用いられるグラフである，横軸が設定する因子数，縦軸が評価指標となるVSSの値，そしてグラフ中の1から4までがそれぞれVSS1，VSS2，VSS3，VSS4の値を示している．
ここでVSS1～VSS4とは，各項目に反映されると仮定する因子の数である．つまり，VSS1は1つの項目に反映されているのはあくまで1つの因子だけという仮定の下（もっとも単純な因子構造！）で，因子数を順に増やしたときにVSSの値がどう推移するかを示している．VSS2の場合は，1つの項目に対して2つの因子が反映されているという仮定の下でのVSSの値の推移である．同様にVSS3は1つの項目に3つの因子が，VSS4の場合1つの項目に4つの因子が反映されている，という仮定の下でのVSSの値の推移である[参考]．
質問項目を作成するときの理想的な状況を考えるなら，1つの質問項目には1つのみ因子が反映されているということが理想である．そのため，VSS1の値が最も高いところが因子数として適切であると考えられる．すなわち，今回だとVSS1だと1因子が示唆される．
ただ実際にはそこまで理想的な状況が得られることは少ないため，VSS2やVSS3，VSS4の値も参考にしながら因子数を決定することが一般的である．例えば，VSS1よりもVSS2が全般に高く，VSS2よりVSS3が全般に高い，さらにVSS4はわずかだがVSS3より高い．そしてVSS4のなかでも特に5因子の時が一番高い値となっている．
こうしたことから，5を因子数の候補として採用することが適切であると考えられる．
メッセージと表の見方
一方メッセージの方を見ると， まずはVSS1とVSS2の結果が表記されている．それぞれ1と2が示唆されるという結果である．これについては上記の説明の通りである．
続いて，MAP (Minimum average partial; 最小平均偏相関)の結果が表記されている．MAPという指標が最も小さくなる因子数を候補として選択するというもので，今回だとこれは2を示唆している．
さらにBIC (Bayesian Information Criterion; ベイズ情報量規準)の結果が表記されている．BICはモデルの複雑さとデータの適合度をバランスさせた指標で，値が小さいほど良いモデルであるとされる．今回だとこれは6を示唆している．
他にもサンプル数の影響を調整したBICも表示されており，この場合には8が示唆されている．ただ，ここでの8は今回の場合だとシミュレーションの最大因子数として設定されている値でもあるので，実際にはこれよりも多い因子数の方がマッチしている可能性もある．
表はこれらの実際の指標値を示したものである．表ではほかにRMSEAなどの様々な適合度指標も表示されているが，それらまで見るのはあまり一般的ではないので，ここでは省略する．
2.4.2 因子分析の実施
とりあえずfa.parallel()関数やvss()関数の結果から5を最初の候補として因子分析を行なってみよう．
因子分析の実施方法は先ほどと同じく，psychパッケージのfa()関数を使う．
5因子
とりあえず因子構造だけ把握したいので，細かな結果は表示させずに，因子負荷量だけを図で表示させる．
fa_result_5 <- fa(data[-1], 5)
fa.diagram(fa_result_5,cut=.4, digits=2, simple = FALSE)
Q10はどの因子に対しても0.4の基準を満たさなかった．
6因子
今度は6因子構造で実施してみよう．
fa_result_6 <- fa(data[-1], 6)
fa.diagram(fa_result_6,cut=.4, digits=2, simple = FALSE)
すると，このように6つ目の因子はいずれの質問項目にも0.4を超える因子負荷量が得られないという結果となった．
多すぎる因子数を設定すると，このようにいずれの質問項目にも因子負荷量が得られない因子が出てくることがある．したがって今回は6以上の因子数というのは候補から外れることになる．
4因子
今度は，因子数を4としてみよう．
fa_result_4 <- fa(data[-1], 4)
fa.diagram(fa_result_4,cut=.4, digits=2, simple = FALSE)
どうやら5因子の時と比べて，Q8,
Q9が1つ目に表示されている因子に吸収されたようである．
Q10は今度もいずれにも因子負荷量が0.4を超えない結果となった．
3因子
さらに，因子数を3にしてみる．
fa_result_3 <- fa(data[-1], 3)
fa.diagram(fa_result_3,cut=.4, digits=2, simple = FALSE)
こうすると，Q19, Q20が3つ目の因子に吸収されたようである．
今度もQ10はいずれの因子とも結び付かなかった．
まとめ
この辺りで，これらの結果を総合的に判断すると，どうやらQ10はそもそもどの因子にも関連していない項目だった，と考える方が良さそうである．
2.4.3 回転と推定方法
因子分析を実行するfa()関数には第3引数，第4引数を設定することが可能である．
これらが一体どういうものかについてはやや高度な内容となるため，ここではこれらの引数の設定の仕方のみを説明する．
第3引数 回転方法（rotate）
| 回転方法 | タイプ | 特徴 | 使用目的 |
|---|---|---|---|
varimax |
直交回転 | 解釈しやすく最も広く使われる | 因子の独立性を保ちたい場合 |
quartimax |
直交回転 | 単一因子の影響を強調 | 単純構造を強調したい場合 |
promax |
斜交回転 | 計算効率が良く，大規模データに適している | 因子の相関を許容したい場合 |
oblimin |
斜交回転 | 因子の相関を許容しつつ適応的に単純化 | 現実のデータに近いモデルを探る（デフォルト） |
simplimax |
直交・斜交両方 | 特定の変数を特定の因子に結びつける | 解釈を明確化したい場合 |
直交回転と斜交回転の違いは，因子間の相関を許容するかどうかである．直交回転は因子間の相関を許容しない(因子間の相関は0であるという前提おく)ため，因子間の独立性を保つことができる．一方，斜交回転は因子間の相関を許容するため，因子間の関係性を考慮した因子構造を推定することができる．
第4引数 推定方法（fm）
| 手法 | 名称 | 説明 |
|---|---|---|
ml |
最尤法（Maximum Likelihood） | 尤度を最大化する形で因子を推定する．データが正規分布に従う場合に有効 |
minres |
最小二乗法（Minimum Residual） | 誤差分散を最小化する形で因子を推定する．因子構造のシンプル化に適している（デフォルト） |
pa |
主軸法（Principal Axis Factoring） | 相関行列の固有値分解に基づき，共通因子の抽出を行う |
uls |
最小二乗法（Unweighted Least Squares） | 標準誤差を考慮せず，単純な回帰モデルに基づいて因子を推定する |
wls |
重み付き最小二乗法（Weighted Least Squares） | 誤差分散を考慮した最小二乗推定法で因子を推定する．信頼性の高い推定が可能 |
gls |
一般化最小二乗法（Generalized Least Squares） | 共分散構造に基づく推定方法で，観測変数の分散の違いを考慮 |
pc |
主成分法（Principal Components） | 因子分析とは異なり，データの分散を最大限に説明する成分を抽出する手法 |
これらから適宜選んで実行することになる．どれを選べば良いか正直悩ましいものがあるが，一般的には，以下の組み合わせがよく用いられる．
oblimin回転とminres推定（デフォルトの設定）varimax回転とpa推定promax回転とml推定
実際にvarimax回転とpa推定と，promax回転とml推定で，5因子構造での分析してみよう．
まずはvarimaxとpaで．
fa_result_v_pa <- fa(data[-1], 5, rotate="varimax", fm="pa")
fa.diagram(fa_result_v_pa,cut=.4, digits=2, simple = FALSE)
続いてpromaxとmlで．
fa_result_p_ml<- fa(data[-1], 5, rotate="promax", fm="ml")
fa.diagram(fa_result_p_ml,cut=.4, digits=2, simple = FALSE)
varimax回転は因子間に相関を仮定しないため，それぞれの因子間に相関係数は表示されない．一方で，promax回転は因子間に相関を仮定するため，それぞれの因子間の相関係数も表示される．
またpromaxの場合，因子負荷量が1を超える値になることもある．実際に上記の図でもML3とQ18や，ML1とQ19は1を超えている．
それぞれ因子負荷量が微妙に変わっているのかわかるだろう．このように回転や推定方法の設定を変えることによって，因子負荷量が変わるので，場合によってはどの因子にも負荷量が0.4を超えなかった項目で0.4を超える負荷量が出てくることもある．
したがって，デフォルトの設定（Oblimin回転,minres推定）でわずかに0.4を超えなかったものがあった場合には，他の回転方法を試してみると良い．
今回の場合，結局回転や推定方法を変えてもQ10はどの因子にも結びつかなかった．したがって，Q10はやはりどの因子にも関連していないと考えるのが妥当である．
2.4.4 適合度指標
因子分析の結果をprint()すると，前半は因子負荷量が表示されるが，後半では推定された因子構造がどの程度データと適合しているかを評価するための指標（適合度指標）が出力される．
様々な適合度指標が表記されているが，代表的なものは以下のものになる．
TLI（Tucker Lewis Index）は，1に近いほど適合度が高いことを示す．一般には0.95以上であれば適合度が高いと解釈され，逆に0.90を下回る場合，適合度が低いとされることが多い．
RMSEA（Root Mean Square Error of Approximation）は，0に近いほど適合度が高いことを示す．一般には0.05以下であれば適合度が高いとされ，逆に0.10を超える場合，適合度が低いとされることが多い．
BIC (Bayesian Information Criterion)は絶対的な評価基準はなく，比較のために使われることが多い．BICが小さいほど適合度が高いとされる．
これらの値を使って妥当な因子構造を探ることもできる．先ほどの３〜6因子構造のそれぞれで適合度指標を見てみよう．
3因子構造
##
## Factor analysis with Call: fa(r = data[-1], nfactors = 3)
##
## Test of the hypothesis that 3 factors are sufficient.
## The degrees of freedom for the model is 168 and the objective function was 1.5
## The number of observations was 986 with Chi Square = 1465.34 with prob < 1.1e-205
##
## The root mean square of the residuals (RMSA) is 0.04
## The df corrected root mean square of the residuals is 0.05
##
## Tucker Lewis Index of factoring reliability = 0.846
## RMSEA index = 0.088 and the 10 % confidence intervals are 0.084 0.093
## BIC = 307.2
## With factor correlations of
## MR1 MR2 MR3
## MR1 1.00 -0.44 0.42
## MR2 -0.44 1.00 -0.58
## MR3 0.42 -0.58 1.004因子構造
##
## Factor analysis with Call: fa(r = data[-1], nfactors = 4)
##
## Test of the hypothesis that 4 factors are sufficient.
## The degrees of freedom for the model is 149 and the objective function was 0.87
## The number of observations was 986 with Chi Square = 844.22 with prob < 1.4e-97
##
## The root mean square of the residuals (RMSA) is 0.03
## The df corrected root mean square of the residuals is 0.04
##
## Tucker Lewis Index of factoring reliability = 0.907
## RMSEA index = 0.069 and the 10 % confidence intervals are 0.064 0.073
## BIC = -182.94
## With factor correlations of
## MR1 MR2 MR4 MR3
## MR1 1.00 -0.42 0.60 -0.19
## MR2 -0.42 1.00 -0.59 0.42
## MR4 0.60 -0.59 1.00 -0.47
## MR3 -0.19 0.42 -0.47 1.005因子構造
##
## Factor analysis with Call: fa(r = data[-1], nfactors = 5)
##
## Test of the hypothesis that 5 factors are sufficient.
## The degrees of freedom for the model is 131 and the objective function was 0.62
## The number of observations was 986 with Chi Square = 605.07 with prob < 5.2e-62
##
## The root mean square of the residuals (RMSA) is 0.02
## The df corrected root mean square of the residuals is 0.03
##
## Tucker Lewis Index of factoring reliability = 0.928
## RMSEA index = 0.061 and the 10 % confidence intervals are 0.056 0.066
## BIC = -298
## With factor correlations of
## MR2 MR4 MR1 MR5 MR3
## MR2 1.00 -0.38 -0.59 -0.39 0.42
## MR4 -0.38 1.00 0.53 0.67 -0.15
## MR1 -0.59 0.53 1.00 0.55 -0.46
## MR5 -0.39 0.67 0.55 1.00 -0.21
## MR3 0.42 -0.15 -0.46 -0.21 1.006因子構造
##
## Factor analysis with Call: fa(r = data[-1], nfactors = 6)
##
## Test of the hypothesis that 6 factors are sufficient.
## The degrees of freedom for the model is 114 and the objective function was 0.41
## The number of observations was 986 with Chi Square = 394.79 with prob < 1.3e-32
##
## The root mean square of the residuals (RMSA) is 0.02
## The df corrected root mean square of the residuals is 0.03
##
## Tucker Lewis Index of factoring reliability = 0.951
## RMSEA index = 0.05 and the 10 % confidence intervals are 0.045 0.055
## BIC = -391.09
## With factor correlations of
## MR2 MR4 MR1 MR5 MR3 MR6
## MR2 1.00 -0.34 -0.56 -0.37 0.40 -0.08
## MR4 -0.34 1.00 0.52 0.66 -0.13 -0.05
## MR1 -0.56 0.52 1.00 0.56 -0.47 0.07
## MR5 -0.37 0.66 0.56 1.00 -0.22 0.01
## MR3 0.40 -0.13 -0.47 -0.22 1.00 -0.06
## MR6 -0.08 -0.05 0.07 0.01 -0.06 1.00結果をまとめると，以下の表の通りとなる．
| 因子数 | TLI | RMSEA | BIC |
|---|---|---|---|
| 3 | 0.85 | 0.09 | 307.2 |
| 4 | 0.91 | 0.07 | -182.94 |
| 5 | 0.93 | 0.06 | -298 |
| 6 | 0.95 | 0.05 | -391.09 |
これらの結果から，適合度指標としては6因子構造が最も適合している，という結果となる． ただ，先に述べた通り，6因子構造は第6因子がいずれの項目とも0.4以上の負荷量を持たなかったため，除外される． 5因子構造では，TLIにせよRMSEAにせよ，いずれも適合度指標としては完全に十分とはいえないが，低いとされる基準よりは上であるため，5因子構造で良さそうである．
2.4.5 因子の解釈
最後は因子の解釈である．因子分析の結果から得られた因子をどのように解釈するかは，その因子と結びついている質問項目を見ていくことが基本である．ここで解釈が難しい結果となるならば，因子構造を見直すことになる．
改めて今回の5因子構造を整理すると以下の通りとなる．
- Q5 上司は私のキャリアの発達を優先してくれている
- Q3 上司は何か悪いことが起ころうとしていると，そのことに気づき，伝えてくれる
- Q4 上司は，私が困難な状況に対応する際には，「私自身が最良と感じる方法で対応すれば良い」と言ってくれる
- Q2 上司は自職場だけでなく病院全体や地域全体に貢献することの重要性を強調する
- Q1 上司は私が落ち込んでる時にはそのことに気づくことができる
- Q6 上司 は上司自身の利益よりも私自身の利益を優先する
- Q7 上司はいつも誠実である
これらの結果から第1因子は「支援型リーダーシップ性」と整理できるだろう(実はこの項目は実際には「サーバントリーダーシップ」としてLiden et al. (2008)によって開発された項目である）．
- Q14 上司はタスクが達成できないと厳しく叱責してくる
- Q13 上司は部下を厳しく指導している
- Q15 上司は「我々の部署は組織の全部署のなかで，一番の業績を上げなければならない」と強調してくる
- Q12 私は上司とともに仕事をするときにはプレッシャーを感じる
- Q16 我々はモノゴトを進める際には上司のルールに従わなければならないさもないと厳しく罰せられる
- Q11 上司はいつも従業員の前では威圧的なふるまいをする
これらの結果から第2因子は「圧迫型リーダーシップ性」と整理できるだろう．
- Q18 上司は私の（私たちの）業務を数字でしか捉えていないと感じる
- Q17 上司は私の（私たちの）業務の実態を分かっていないと感じる
- Q22 上司に私の（私たちの）業務実態について話ができる機会・場が設けられていない
- Q21 上司は私の（私たちの）ことよりも，上役の方を見ていると感じる
これらの結果から第3因子は「現場軽視型リーダーシップ性」と整理できるだろう．
- Q8 上司は私に対して指示に完全に従うように言ってくる
- Q9 上司はそれが重要なものかどうかにかかわらず，組織におけるあらゆることを自分で決めようとする
これらの結果からは第4因子は「自己中心型リーダーシップ性」と整理できるだろう．
- Q19 上司は私の（私たちの）業務を遂行している現場によく顔を出す
- Q20 上司は私の（私たちの）業務の状況をよく聞きにくる
最後とにこれら2つから第5因子は「現場重視型リーダーシップ性」と整理できるだろう．
このように，因子分析の結果から得られた因子を質問項目と結びつけて解釈することが重要である． 今回は上記の通りに解釈することができたが，場合によっては解釈に悩むような結果が出てくることもある． そのような場合には因子数を見直したり，回転や推定方法を設定し直したり，あるいは除外した方が良い項目を除外して，因子分析を再実行してみることが必要になる．
逆に，もし解釈がしやすい結果なのであれば，適合度がそれよりも高い因子構造があったとしても，解釈がしやすいものを選択することもよく行われる．例えば今回の結果では5因子構造よりも4因子構造の方が解釈しやすいならば，4因子構造を選択することもあり得る．
ちなみに4因子構造を仮定した場合，「自己中心型リーダーシップ」と「圧迫型リーダーシップ」が一つの因子としてまとまる． 実は，これら2つは元々「権威主義的リーダーシップ」という一つの概念を測定する質問項目としてChen et al.(2004)において開発されたものである．したがって，あくまでもともとは1つの因子であるという点を重視するならば，適合度としては5因子構造よりも劣るものの，4因子構造を採用することも十分にあり得ることである．
2.4.6 探索的因子分析のまとめ
以上，具体例を用いて探索的因子分析の進め方を紹介した． 実際には探索的という言葉の通り，手順が明確に定まっているわけではなく，因子数を変えたり，回転方法を変えたり，推定方法を変えたりしつつ，適合度指標を見たり，結果の解釈可能性を考えたりを繰り返しながら試行錯誤していくことになる．正解というものが存在しないため，「より良い結果」を求めて延々と分析を繰り返してしまいがちであるが，ある程度自分として納得できる結果が得られたら，そこで割り切ることも重要である．
次のデータは実際に著者が指導したゼミ生が2021年度に行った優柔不断度を測るアンケート結果である． アンケート項目は以下の通りである．
- Q1: 何かを選ぶ時は些細なことは気にしない
- Q2: 重要なことを決めるとき，些細な問題にこだわり多くの時間を費やす
- Q3: 決断を先送りにしようとする
- Q4: 重要な問題は，直前まで決定を延ばす
- Q5: 他の人が何を選ぶのか気になる
- Q6: 自分が選んだものが他の人と違うと不安になる
- Q7: 決断した後，その決断に後悔しない
- Q8: 決めた後で，間違ってしまったと思うことがよくある
このデータを使って，以下の2つのことを行え．
- 主成分分析を実施し，累積説明率が60％を超えたところまでの主成分について，その内容を質問項目と照らしてを解釈せよ．
- 探索的因子分析を実施し，適切だと思う因子構造を各自で見つけ，推定された因子の内容を質問項目に照らして解釈せよ．
data <- read.csv("./practice/13_work_fa.csv")
head(data)## time_stamp Q1 Q2 Q3 Q4 Q5 Q6 Q7 Q8
## 1 2021/11/2 11:19 5 5 2 5 1 2 4 4
## 2 2021/11/2 11:21 2 4 2 4 4 3 4 2
## 3 2021/11/2 11:21 4 2 4 5 3 4 3 5
## 4 2021/11/2 11:22 2 2 4 5 5 4 4 4
## 5 2021/11/2 11:23 2 5 3 5 4 4 2 4
## 6 2021/11/2 11:25 1 4 4 3 5 5 2 4# colnames(data) <- c("time_stamp","Q1","Q2","Q3","Q4","Q5","Q6","Q7","Q8")
#
# head(data[-c(10,11)])
#
# data <- data[-c(10,11)]
# write_csv(data,"./practice/13_work_fa.csv")
pc_res<-prcomp(data[-1], scale. = TRUE)
summary(pc_res)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.7219 1.0947 0.9793 0.9677 0.79811 0.73633 0.63820
## Proportion of Variance 0.3706 0.1498 0.1199 0.1171 0.07962 0.06777 0.05091
## Cumulative Proportion 0.3706 0.5204 0.6403 0.7574 0.83700 0.90478 0.95569
## PC8
## Standard deviation 0.59540
## Proportion of Variance 0.04431
## Cumulative Proportion 1.00000biplot(pc_res)
round(pc_res$rotation,2)## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## Q1 -0.31 -0.05 -0.61 -0.34 0.47 -0.43 -0.08 -0.03
## Q2 0.35 0.20 0.45 0.27 0.58 -0.48 0.00 0.05
## Q3 0.41 0.35 -0.17 -0.33 -0.24 -0.12 -0.14 0.69
## Q4 0.33 0.55 -0.07 -0.38 0.16 0.32 0.00 -0.56
## Q5 0.38 -0.18 -0.41 0.46 0.08 0.17 -0.64 -0.09
## Q6 0.41 -0.06 -0.42 0.27 -0.23 -0.27 0.65 -0.18
## Q7 -0.31 0.47 -0.23 0.45 0.30 0.40 0.26 0.32
## Q8 0.32 -0.53 0.01 -0.26 0.46 0.45 0.28 0.24fa.parallel(data[-1], fa="fa")
## Parallel analysis suggests that the number of factors = 4 and the number of components = NAvss(data[-1])
##
## Very Simple Structure
## Call: vss(x = data[-1])
## VSS complexity 1 achieves a maximimum of 0.64 with 1 factors
## VSS complexity 2 achieves a maximimum of 0.75 with 3 factors
##
## The Velicer MAP achieves a minimum of 0.05 with 1 factors
## BIC achieves a minimum of -9.52 with 4 factors
## Sample Size adjusted BIC achieves a minimum of -3.17 with 4 factors
##
## Statistics by number of factors
## vss1 vss2 map dof chisq prob sqresid fit RMSEA BIC SABIC complex
## 1 0.64 0.00 0.046 20 1.7e+02 3.9e-26 4.7 0.64 0.16 57.2 120.6 1.0
## 2 0.54 0.73 0.072 13 9.3e+01 3.4e-14 3.5 0.73 0.14 18.7 59.9 1.4
## 3 0.55 0.75 0.114 7 3.2e+01 3.8e-05 2.6 0.80 0.11 -8.0 14.2 1.6
## 4 0.57 0.75 0.150 2 1.9e+00 3.8e-01 1.9 0.85 0.00 -9.5 -3.2 1.5
## 5 0.56 0.73 0.278 -2 8.2e-08 NA 1.5 0.89 NA NA NA 1.5
## 6 0.59 0.74 0.455 -5 0.0e+00 NA 1.5 0.88 NA NA NA 1.5
## 7 0.59 0.73 1.000 -7 3.3e-11 NA 1.5 0.88 NA NA NA 1.5
## 8 0.59 0.73 NA -8 3.3e-11 NA 1.5 0.88 NA NA NA 1.5
## eChisq SRMR eCRMS eBIC
## 1 1.7e+02 1.0e-01 0.118 56.6
## 2 8.1e+01 6.9e-02 0.101 6.9
## 3 3.1e+01 4.3e-02 0.085 -8.7
## 4 9.9e-01 7.6e-03 0.028 -10.5
## 5 4.0e-08 1.5e-06 NA NA
## 6 2.3e-15 3.6e-10 NA NA
## 7 1.3e-11 2.7e-08 NA NA
## 8 1.3e-11 2.7e-08 NA NAfa_result <- fa(data[-1], 4, rotate="varimax", fm="pa")## maximum iteration exceededfa.diagram(fa_result,cut=.4, digits=2, simple = FALSE)