Rで因子分析やってみた

今更感ありますが、R言語で因子分析を試してみたのでそのメモを書きます。

# 因子分析とは

>因子分析とは、ある観測された変数（質問項 目への回答など）が、どのような潜在的な変 数（観測されない、仮定された変数）から影響 を受けているかを探る手法。 

 [因子分析](cogpsy.educ.kyoto-u.ac.jp/personal/Kusumi/datasem13/masuda.pdf)より引用

ざっくり言うと、複数の観測変数の背後にある潜在的な要因を明らかにするための分析です。
例えば、定期試験で`国語、社会、英語、数学、理科`のテストを行ったとします。
もし、`数学、理科`の点数が高い人がいれば、理系科目が得意そうですよね。
この「`数学、理科`の点数が高い」という背景にある**理系能力**という潜在的な因子を探るのが**因子分析**です。

## 分析用データ

今回分析するデータは、統計教育推進委員会が提供する`国語, 英語, 数学, 物理, 化学`の5科目の成績データをダウンロードして使用します。このデータを使ってR言語で因子分析を行い、5科目の成績に影響を与える潜在因子を探索したいと思います。

[5科目の成績データ](https://www.google.co.jp/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&ved=0ahUKEwiYj_DtuK_QAhUBoZQKHW9RC-cQFggbMAA&url=http%3A%2F%2Festat.sci.kagoshima-u.ac.jp%2Fdata%2Fcgi-bin%2Fdata%2Fwhats_data%2Fdata%2Fimg%2F932722923_9821.xls&usg=AFQjCNGqWn0acSBm0XaieWtnDGrZXi5GrA&sig2=xRSveyRwT6ak2n5spbopbQ&cad=rja)

# データ読み込み

先ほどダウンロードした[5科目の成績データ](https://www.google.co.jp/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&ved=0ahUKEwiYj_DtuK_QAhUBoZQKHW9RC-cQFggbMAA&url=http%3A%2F%2Festat.sci.kagoshima-u.ac.jp%2Fdata%2Fcgi-bin%2Fdata%2Fwhats_data%2Fdata%2Fimg%2F932722923_9821.xls&usg=AFQjCNGqWn0acSBm0XaieWtnDGrZXi5GrA&sig2=xRSveyRwT6ak2n5spbopbQ&cad=rja)をR言語で記述したプログラムから読み込みます。このとき、ダウンロードしたファイルを`csv`形式に変換し、不要な空白のセル等は削除しておくと良いです。
データを読み込むプログラムは以下の通りです。

```

data <- read.csv(file="ファイルパス.csv",
                   header=T, fileEncoding="Shift-JIS", row.names=1)

```

# 因子数の決定

抽出する因子数を決定するために、相関行列を算出した後、その行列の固有値を算出します。


* 相関行列
 * 相関行列は関数`cor()`で算出することができます。
 
```

cor <- cor(data)

#出力
cor

```

相関行列
<img width="400" alt="スクリーンショット 2016-11-17 23.58.02.png" src="https://qiita-image-store.s3.amazonaws.com/0/134515/6dd8255c-e003-4b64-6b3d-7387d6861de7.png">


* 固有値
 * 固有値および固有ベクトルは関数`eigen()`で算出することができます。
  * 固有値は`$values`,固有ベクトルは`$vectors`に格納されます。

```

eigen <- eigen(cor)$values

#出力
eigen

```

固有値
<img width="396" alt="スクリーンショット 2016-11-18 0.11.59.png" src="https://qiita-image-store.s3.amazonaws.com/0/134515/9523ca22-e966-682d-f30e-edec32661a0a.png">

* スクリープロット
 * 先ほど算出した固有値の減衰をスクリープロットで見て因子数を決定します。
 * スクリープロットの描画は関数`plot()`を使用します。
 
```

#type="b"でグラフ上に点と線両方を表示
plot(eigen, type="b")

```

スクリープロット
<img width="454" alt="スクリーンショット 2016-11-18 0.30.17.png" src="https://qiita-image-store.s3.amazonaws.com/0/134515/f84fcf44-9252-a8c7-5f9f-77bf5f88d167.png">

### 因子数

算出された固有値のうち、1より大きい固有値は2つで、スクリープロットの固有値減衰をみると、第3固有値以降は減衰が緩やかになっていますので、今回は因子数を**2**として因子分析を行います。


# 因子分析

R言語には因子分析を行うための標準関数が用意されています。
因子分析は関数`factanal()`で行うことができます。
なお、標準関数では初期解の算出方法は最尤法のみ、因子の回転法はバリマックス回転(直行回転)かプロマックス回転(斜行回転)のみとなっております。

* 直行回転
 * 因子間の相関を仮定しない
* 斜行回転
 * 因子間の相関を仮定する

(psychパッケージ等を利用すれば初期解の算出法や因子の回転法などの選択の幅が広がります)


今回は、バリマックス回転法で因子分析を行いたいと思います。

```

factanal <- factanal(x=data,factors=2,rotation="varimax")
#出力(cutoff=0で全ての因子負荷量を表示)
print(factanal,　cutoff=0)

```

結果
<img width="421" alt="スクリーンショット 2016-12-04 19.57.59.png" src="https://qiita-image-store.s3.amazonaws.com/0/134515/ad1418de-3cad-6545-33f5-4d92e433c941.png">


### 結果の見かた

* **uniquenesses**
    * factanal()では共通性の値が出力されないかわりに独自性(1-共通性)の値が出力されます。
    * 共通性は、各測定値に対して、共通因子の部分がどの程度あるのかを示す指標。

今回の分析結果では、どの項目も共通因子からある程度影響を受けていることがわかります。


* **Loadings**
 * 因子負荷
  * 観測変数に対して共通因子がどの程度影響を与えているかを示したもの。
 
今回の分析結果では、第一因子は`数学、物理、化学`が高い因子負荷量を示しており、第二因子は`国語、英語`が高い因子負荷量を示していることがわかります。

* **SS loadings**
 * 因子寄与
  * その因子が全体分散のうちどの程度を説明しているか。

* **Proportion Var**
 * 因子寄与率
  * その因子が全体分散のうちどの程度の割合を説明しているか。

* **Cumulative Var**
 * 累積寄与率
  * 因子寄与率を大きい順に足していったもの。
  * そこまでの因子で全体の情報量のどの程度が説明されているかを示す。
  * 今回の分析結果では、2つの因子で5項目の全分散を説明する割合が78%であることがわかります。

## 分析結果まとめ

* 固有値の変化(`3.1188278 1.1563143 0.3410577 0.2036203 0.1801800`)より、**2因子構造が妥当**
* 第一因子は`数学、物理、化学`が高い因子負荷量を示しており、第二因子は`国語、英語`が高い因子負荷量を示していることから、第一因子を**理系能力**、第二因子を**文系能力**と命名します。
* ２因子で全体の分散を説明する割合は78%。


以上です。Rを使えば中の数式がわからなくても簡単に因子分析ができます。
今後は`psychパッケージ`などを使ってもっと色々試してみたいと思います。