8  第11章 探索的因子分析

pacman::p_load(
  psych, # 因子分析のためのパッケージ
  GPArotation, # 因子分析の回転方法を提供するパッケージ
  tidyverse, # データ操作と可視化のためのパッケージ群
  readxl, # Excelファイルの読み込みのためのパッケージ
  gt, # 表の作成のためのパッケージ
  gtExtras, # gtパッケージの拡張機能を提供するパッケージ
  dendextend, # 階層的クラスタリングのためのパッケージ
  cluster, # クラスタリングのためのパッケージ
  factoextra, # クラスタリングの可視化のためのパッケージ
  useful, # PCAのためのパッケージ
  ggrepel # ラベルを重ならなくする
  )

The downloaded binary packages are in
    /var/folders/lf/9kvp0m950bx72lzg24t53g280000gn/T//RtmpltLXdD/downloaded_packages

8.1 本性の概要

8.2 因子分析とは

因子分析(factor analysis)は、観測された変数の背後に存在する潜在的な共通因子(common factor)を明らかにするための多変量解析の手法です。 因子分析には2つのアプローチがあります。

  1. 探索的因子分析：データから潜在的な因子構造を探索する方法
  2. 確認的因子分析：因子数と因子負荷量について、事前に仮説を立てて分析する方法

テキストでは、探索的因子分析について説明します。

8.2.1 因子分析のモデル構造

教科書のp.305からの解説を参照

8.2.2 因子分析の特徴

たとえば，お店を訪れたお客さんに10個の質問からなるアンケートに答えてもらい，1000名分のデータが集まったとしましょう。 いま，皆さんの手には以下のような1000人分の10個の回答結果のデータがあります。

ID Q1 Q2 Q3 Q4 Q5 Q6 Q7 Q8 Q9 Q10
1 5 3 4 2 3 4 4 2 4 2
2 4 4 5 3 4 5 5 3 5 3
… … … … … … … … … … …
1000 3 2 3 1 2 3 3 1 3 1

このとき，IDがi = 1のQ1の解答 5 が，1さんが心の中に持っている価値観（例：お店の清潔さ）に従って決められたものであると考えることができます。 このように，回答者の心の中にある価値観や性格などの潜在的な因子が，質問項目への回答を決定していると考えることができます。 この潜在的な因子を特定する統計的手法の1つが因子分析(factor analysis)です。 上の表のように，1000人が答えた10個の質問項目の回答から，少数の潜在的因子を特定することで，少ない数の因子でデータを説明することができます。 1000人の回答は1000通りある，というのではアンケート調査から分かることはありませんが， 少数の共通因子を特定することで，1000人の回答を理解しやすくまとめることができます。 このように，因子分析は，観測された変数の背後に存在する潜在的な共通因子を明らかにするための多変量解析の手法です。

8.3 2因子モデル

因子の数を2と仮定したモデルを2因子モデルといいます。 ここで，j = 1, 2, \dots , J 個の質問項目をもつアンケート調査をi, i = 1, 2, ..., n人の回答者に実施したとします。

  • \boldsymbol{x}_i：観測変数。実際にデータとして得られた質問項目（例：テストの点数）
  • \boldsymbol{f}_i：共通因子。目に見えない潜在的な能力や性質（例：計算力、読解力）
  • \boldsymbol{\alpha}：因子負荷量（Factor Loading）。各共通因子が観測変数に与える影響の強さ
  • \boldsymbol{\varepsilon}_i：独自因子（Unique Factor）。その項目特有の誤差や要因

J個の項目の背後には2つの潜在因子 f_1 と f_2 が存在すると仮定して、以下のモデルを考えます。 iさんのj番目の質問項目への回答をx_{ji} とします。 今手元にあるデータは，x_{ji} だけです。 その質問項目への回答は次のように表すとします。 \begin{aligned} \begin{cases} x_{1i} &= \alpha _{11} f_{1i} + \alpha _{12} f_{2i} + \varepsilon_{1i} \\ \ \vdots \\ x_{Ji} &= \alpha _{J1} f_{1i} + \alpha _{J2} f_{2i} + \varepsilon_{Ji} \end{cases} \end{aligned}

掛けて足す、の構造をもっているので、行列で表すことができます。

\begin{aligned} \underbrace{ \begin{pmatrix} x_{1i} \\ \vdots \\ x_{Ji} \end{pmatrix} }_{J\times1} &= \underbrace{ \begin{pmatrix} \alpha _{11} & \alpha _{12} \\ \vdots & \vdots \\ \alpha _{J1} & \alpha _{J2} \end{pmatrix} }_{J\times2} \underbrace{ \begin{pmatrix} f_{1i} \\ f_{2i} \end{pmatrix} }_{2\times1} + \underbrace{ \begin{pmatrix} \varepsilon_{1i} \\ \vdots \\ \varepsilon_{Ji} \end{pmatrix} }_{J\times1} \end{aligned}

まとめて書くと、次のようになります。

\boldsymbol{x}_i = \boldsymbol{\alpha} \boldsymbol{f}_i + \boldsymbol{\varepsilon}_i この式を2因子モデルといいます。 任意の非特異な2 \times 2 の正則行列\boldsymbol{T}を用いて、次のように変形します。 \begin{aligned} \boldsymbol{x}_i &= \boldsymbol{\alpha} (\boldsymbol{TT^{-1}}) \boldsymbol{f}_i + \boldsymbol{\varepsilon}_i\\ & = \underbrace{(\boldsymbol{\alpha T})}_{J\times 2, 2\times 2} \underbrace{(\boldsymbol{T^{-1}f}_i)}_{2\times 2, 2\times 1} + \boldsymbol{\varepsilon}_i\\ & = \boldsymbol{\alpha}^* \boldsymbol{f}_i^* + \boldsymbol{\varepsilon}_i \end{aligned}

この式の変形を詳しく見ていきましょう。 これは、因子の回転(factor rotation)と呼ばれる操作で、\boldsymbol{\alpha} を \boldsymbol{\alpha}^* に、\boldsymbol{f}_i を \boldsymbol{f}_i^* に変換することを意味します。 この \boldsymbol{T} をどのように選ぶか，つまりどのように回転させるのかによって、\boldsymbol{\alpha}^* と \boldsymbol{f}_i^* の値が変わります。 ではこの回転について、詳しく見ていきましょう。

8.3.1 因子軸の回転

因子分析には「回転の不定性」という性質があります。 行列 \boldsymbol{T} を挟み込んでも、計算結果として得られる観測変数 \boldsymbol{x}_i の値（あるいは分散共分散構造）は全く変わりません。 これは、潜在因子の座標軸，つまり縦軸と横軸を回転させて新しい因子負荷量 \boldsymbol{\alpha}^* と新しい因子スコア \boldsymbol{f}_i^* を作っていることに相当します。 初期の計算で得られる \boldsymbol{\alpha} は、すべての因子がすべての変数に中途半端に影響していることが多く、解釈が困難です。 そこで、適切な \boldsymbol{T} を選ぶことで、特定の変数には強く、別の変数には弱く影響するように「整理」し直すのです。 この \boldsymbol{\alpha}^* を回転後の因子負荷量といい、\boldsymbol{f}_i^* を回転後の因子得点といいます。

データの読み込みと整理
factor_exdata <- readxl::read_excel("data/factor_ex.xlsx", na = ".")
summary(factor_exdata)
       ID              V1              V2            V3            V4     
 Min.   : 1.00   Min.   :1.000   Min.   :2.0   Min.   :1.0   Min.   :2.0  
 1st Qu.: 8.25   1st Qu.:2.000   1st Qu.:3.0   1st Qu.:2.0   1st Qu.:3.0  
 Median :15.50   Median :4.000   Median :4.0   Median :4.0   Median :4.0  
 Mean   :15.50   Mean   :3.933   Mean   :3.9   Mean   :4.1   Mean   :4.1  
 3rd Qu.:22.75   3rd Qu.:6.000   3rd Qu.:5.0   3rd Qu.:6.0   3rd Qu.:5.0  
 Max.   :30.00   Max.   :7.000   Max.   :7.0   Max.   :7.0   Max.   :7.0  
       V5            V6       
 Min.   :1.0   Min.   :2.000  
 1st Qu.:2.0   1st Qu.:3.000  
 Median :3.5   Median :4.000  
 Mean   :3.5   Mean   :4.167  
 3rd Qu.:5.0   3rd Qu.:4.750  
 Max.   :7.0   Max.   :7.000

読み込んだデータは、変数IDと5つの質問項目V1からV6で構成されています。 質問項目は、1から7のリッカート尺度で回答されており、数値が大きいほど肯定的な回答を意味しています。 V5は逆質問であるため、数値が大きいほど重要度が低いことを意味するので、8 - 得点 として逆転させておきます。

factor_exdata$V5 <- 8 - factor_exdata$V5
summary(factor_exdata)
       ID              V1              V2            V3            V4     
 Min.   : 1.00   Min.   :1.000   Min.   :2.0   Min.   :1.0   Min.   :2.0  
 1st Qu.: 8.25   1st Qu.:2.000   1st Qu.:3.0   1st Qu.:2.0   1st Qu.:3.0  
 Median :15.50   Median :4.000   Median :4.0   Median :4.0   Median :4.0  
 Mean   :15.50   Mean   :3.933   Mean   :3.9   Mean   :4.1   Mean   :4.1  
 3rd Qu.:22.75   3rd Qu.:6.000   3rd Qu.:5.0   3rd Qu.:6.0   3rd Qu.:5.0  
 Max.   :30.00   Max.   :7.000   Max.   :7.0   Max.   :7.0   Max.   :7.0  
       V5            V6       
 Min.   :1.0   Min.   :2.000  
 1st Qu.:3.0   1st Qu.:3.000  
 Median :4.5   Median :4.000  
 Mean   :4.5   Mean   :4.167  
 3rd Qu.:6.0   3rd Qu.:4.750  
 Max.   :7.0   Max.   :7.000

次に、IDを行名に設定します。

IDを行名に設定
factor_exdata <- factor_exdata |>
  column_to_rownames(var = "ID")
glimpse(factor_exdata)
Rows: 30
Columns: 6
$ V1 <dbl> 7, 1, 6, 4, 1, 6, 5, 6, 3, 2, 6, 2, 7, 4, 1, 6, 5, 7, 2, 3, 1, 5, 2…
$ V2 <dbl> 3, 3, 2, 5, 2, 3, 3, 4, 4, 6, 4, 3, 2, 6, 3, 4, 3, 3, 4, 5, 3, 4, 2…
$ V3 <dbl> 6, 2, 7, 4, 2, 6, 6, 7, 2, 2, 7, 1, 6, 4, 2, 6, 6, 7, 3, 3, 2, 5, 1…
$ V4 <dbl> 4, 4, 4, 6, 3, 4, 3, 4, 3, 6, 3, 4, 4, 5, 2, 3, 3, 4, 3, 6, 3, 4, 5…
$ V5 <dbl> 6, 3, 7, 6, 2, 6, 4, 7, 2, 1, 6, 3, 7, 5, 2, 5, 5, 7, 2, 4, 3, 6, 4…
$ V6 <dbl> 4, 4, 3, 5, 2, 4, 3, 4, 3, 6, 3, 4, 3, 6, 4, 4, 4, 4, 3, 6, 3, 4, 4…
knitr::kable(summary(factor_exdata))
V1 V2 V3 V4 V5 V6
Min. :1.000 Min. :2.0 Min. :1.0 Min. :2.0 Min. :1.0 Min. :2.000
1st Qu.:2.000 1st Qu.:3.0 1st Qu.:2.0 1st Qu.:3.0 1st Qu.:3.0 1st Qu.:3.000
Median :4.000 Median :4.0 Median :4.0 Median :4.0 Median :4.5 Median :4.000
Mean :3.933 Mean :3.9 Mean :4.1 Mean :4.1 Mean :4.5 Mean :4.167
3rd Qu.:6.000 3rd Qu.:5.0 3rd Qu.:6.0 3rd Qu.:5.0 3rd Qu.:6.0 3rd Qu.:4.750
Max. :7.000 Max. :7.0 Max. :7.0 Max. :7.0 Max. :7.0 Max. :7.000 
# コード11-3
knitr::kable(cor(factor_exdata), caption = "相関行列")
相関行列
V1 V2 V3 V4 V5 V6
V1 1.0000000 -0.0532178 0.8730902 -0.0861622 0.8576366 0.0041681
V2 -0.0532178 1.0000000 -0.1550200 0.5722121 -0.0197456 0.6404649
V3 0.8730902 -0.1550200 1.0000000 -0.2477879 0.7778480 -0.0180688
V4 -0.0861622 0.5722121 -0.2477879 1.0000000 0.0065819 0.6404649
V5 0.8576366 -0.0197456 0.7778480 0.0065819 1.0000000 0.1364029
V6 0.0041681 0.6404649 -0.0180688 0.6404649 0.1364029 1.0000000 
# コード11-4
cor.exdata <- cor(factor_exdata)

VSS.scree(cor.exdata)

# コード11-5
fa <- fa(r = factor_exdata, nfactors = 2,
  rotate = "promax", fm = "ml")
fa
Factor Analysis using method =  ml
Call: fa(r = factor_exdata, nfactors = 2, rotate = "promax", fm = "ml")
Standardized loadings (pattern matrix) based upon correlation matrix
     ML1   ML2   h2    u2 com
V1  0.97  0.00 0.94 0.063   1
V2 -0.02  0.75 0.56 0.437   1
V3  0.89 -0.12 0.83 0.174   1
V4 -0.06  0.78 0.62 0.378   1
V5  0.89  0.11 0.79 0.205   1
V6  0.08  0.83 0.69 0.309   1

                       ML1  ML2
SS loadings           2.54 1.89
Proportion Var        0.42 0.32
Cumulative Var        0.42 0.74
Proportion Explained  0.57 0.43
Cumulative Proportion 0.57 1.00

 With factor correlations of 
      ML1   ML2
ML1  1.00 -0.06
ML2 -0.06  1.00

Mean item complexity =  1
Test of the hypothesis that 2 factors are sufficient.

df null model =  15  with the objective function =  4.25 with Chi Square =  111.31
df of  the model are 4  and the objective function was  0.21 

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

The harmonic n.obs is  30 with the empirical chi square  0.33  with prob <  0.99 
The total n.obs was  30  with Likelihood Chi Square =  5.21  with prob <  0.27 

Tucker Lewis Index of factoring reliability =  0.95
RMSEA index =  0.095  and the 90 % confidence intervals are  0 0.314
BIC =  -8.39
Fit based upon off diagonal values = 1
Measures of factor score adequacy             
                                                   ML1  ML2
Correlation of (regression) scores with factors   0.98 0.92
Multiple R square of scores with factors          0.96 0.84
Minimum correlation of possible factor scores     0.92 0.68
# コード11-6
fa.diagram(fa)

# コード11-7
fs <- data.frame(fa$scores)
factor_exdata$rowname <- rownames(factor_exdata)
fs$rowname <- rownames(fs)
factor_exdata <- left_join(factor_exdata, fs, by = "rowname")
knitr::kable(head(factor_exdata), caption = "結合後データ")
結合後データ
V1 V2 V3 V4 V5 V6 rowname ML1 ML2
7 3 6 4 6 4 1 1.3105846 -0.2899670
1 3 2 4 3 4 2 -1.2878473 -0.2069844
6 2 7 4 7 3 3 1.1827458 -0.7999135
4 5 4 6 6 5 4 0.1475745 1.0035312
1 2 2 3 2 2 5 -1.3909639 -1.3063607
6 3 6 4 6 4 6 0.9927775 -0.2878395 
# コード11-8

ex_cluster <- factor_exdata %>%
  select(ML1, ML2)
Hier1 <- agnes(ex_cluster, metric = "euclidian", stand = TRUE)
pltree(Hier1)

# コード11-9
cl_1 <- kmeans(ex_cluster, 3)
clus_fa <- data.frame(cl_1$cluster)
clus_fa$rowname <- rownames(clus_fa)
factor_exdata <- left_join(factor_exdata, clus_fa, by = "rowname")
factor_exdata$cl_1.cluster <- factor(factor_exdata$cl_1.cluster)

# Visualizing the clusters with 2 factors
p1 <- ggplot(data = factor_exdata,
             mapping = aes(x = ML1, y = ML2, color = cl_1.cluster))
p1 + geom_point() +
  geom_text_repel(mapping = aes(label = rownames(factor_exdata))) +
  labs(x = "Ease of use", y = "Design")

# コード11-10
factor_exdata3 <- factor_exdata %>%
  dplyr::select(V1:V6) # 質問項目だけを厳密に抽出
pca_res <- pca(factor_exdata3, nfactors = 2, rotate = "none")
pca_res
Principal Components Analysis
Call: principal(r = r, nfactors = nfactors, residuals = residuals, 
    rotate = rotate, n.obs = n.obs, covar = covar, scores = scores, 
    missing = missing, impute = impute, oblique.scores = oblique.scores, 
    method = method, use = use, cor = cor, correct = 0.5, weight = NULL, 
    n.rotations = n.rotations)
Standardized loadings (pattern matrix) based upon correlation matrix
     PC1  PC2   h2    u2 com
V1  0.93 0.25 0.93 0.074 1.1
V2 -0.30 0.80 0.72 0.277 1.3
V3  0.94 0.13 0.89 0.106 1.0
V4 -0.34 0.79 0.74 0.261 1.4
V5  0.87 0.35 0.88 0.122 1.3
V6 -0.18 0.87 0.79 0.210 1.1

                       PC1  PC2
SS loadings           2.73 2.22
Proportion Var        0.46 0.37
Cumulative Var        0.46 0.82
Proportion Explained  0.55 0.45
Cumulative Proportion 0.55 1.00

Mean item complexity =  1.2
Test of the hypothesis that 2 components are sufficient.

The root mean square of the residuals (RMSR) is  0.05 
 with the empirical chi square  1.97  with prob <  0.74 

Fit based upon off diagonal values = 0.99
# コード11-11
# データの相関行列から固有値と固有ベクトルを抽出
cor_ex <- cor(factor_exdata3)
eigen_list <- eigen(cor_ex)

# 固有ベクトルについての情報をオブジェクトとして定義
pc1_eig <- eigen_list$vectors[, 1]
pc2_eig <- eigen_list$vectors[, 2]

# 固有ベクトルに，固有値の平方根を乗じる
PC1 <- pc1_eig * sqrt(eigen_list$values[1])
PC2 <- pc2_eig * sqrt(eigen_list$values[2])
knitr::kable(data.frame(cbind(PC1, PC2)), caption = "計算結果")
計算結果
PC1 PC2
0.9283425 -0.2532285
-0.3005297 -0.7952496
0.9361812 -0.1308894
-0.3415817 -0.7889663
0.8687553 -0.3507939
-0.1766389 -0.8711581 
# 結果の表示（行名をV1-V6に合わせて見やすくします）
loadings_manual <- data.frame(PC1, PC2)
rownames(loadings_manual) <- colnames(factor_exdata3)
knitr::kable(loadings_manual, caption = "主成分負荷量の計算結果")
主成分負荷量の計算結果
PC1 PC2
V1 0.9283425 -0.2532285
V2 -0.3005297 -0.7952496
V3 0.9361812 -0.1308894
V4 -0.3415817 -0.7889663
V5 0.8687553 -0.3507939
V6 -0.1766389 -0.8711581
Back to top