ねこすたっと

ねこの気持ちと統計について悩む筆者の備忘録的ページ。

データ構造を要約・説明する(2):検証的因子分析(lavaanパッケージなど)[R]

因子分析とは

因子分析には以下の2種類がある。

  • 探索的因子分析(exploratory FA, EFA):何も仮説を考えないまま、データの背景にありそうな因子構造を推定するもの
  • 検証的因子分析(confirmatory FA, CFA):先行する理論をもとに因子・観察変数の関係性をあらかじめ設定して、データがその因子構造に当てはまるのかを検証するもの。確認的因子分析とも呼ばれる。

検証的因子分析の手順

  1. データの確認
  2. 因子構造の指定
  3. モデルの当てはめ
  4. 推定結果の確認
  5. モデルの改善

使用するパッケージとデータセット

CFAは主にlavaan*1パッケージを使うが、他に関連するパッケージもまとめて読み込んでおく。

  • lavaanパッケージ:Latent variable analysis = 潜在変数を扱う解析、つまりCFAや構造方程式モデリング(Structural Equation Modeling, SEM)を行える
  • MVNパッケージ:正規性(後述)の確認
  • semToolsパッケージ:信頼性(後述)の確認
  • semPlot:パス図を描く
library(lavaan)
library(MVN)
library(semTools)
library(semPlot)

lavaanパッケージのHolzingerSwineford1939データセットを用いる。

> head(HolzingerSwineford1939)
  id sex ageyr agemo  school grade       x1   x2    x3       x4   x5        x6       x7   x8       x9
1  1   1    13     1 Pasteur     7 3.333333 7.75 0.375 2.333333 5.75 1.2857143 3.391304 5.75 6.361111
2  2   2    13     7 Pasteur     7 5.333333 5.25 2.125 1.666667 3.00 1.2857143 3.782609 6.25 7.916667
3  3   2    13     1 Pasteur     7 4.500000 5.25 1.875 1.000000 1.75 0.4285714 3.260870 3.90 4.416667
4  4   1    13     2 Pasteur     7 5.333333 7.75 3.000 2.666667 4.50 2.4285714 3.000000 5.30 4.861111
5  5   2    12     2 Pasteur     7 4.833333 4.75 0.875 2.666667 4.00 2.5714286 3.695652 6.30 5.916667
6  6   2    14     1 Pasteur     7 5.333333 5.00 2.250 1.000000 3.00 0.8571429 4.347826 6.65 7.500000

含まれる変数は以下のとおり。

  • id:ID
  • sex:性別
  • ageyr:年齢の年の部分
  • agemo:年齢の月の部分
  • school:児童の所属する学校
  • grade:学年
  • x1:visual perception(視覚)
  • x2:cubes(立方体?)
  • x3:lozenges(菱形?)
  • x4:paragraph comprehension(パラグラフ理解)
  • x5:sentence completion(文章補完)
  • x6:word meaning(語彙)
  • x7:speeded addition(速く足す)
  • x8:speeded counting of dots(速く数える)
  • x9:speeded discrimination straight and curved capitals(速く区別する)

x1x9だけ使う。

data <- HolzingerSwineford1939[-c(1:6)]

1. データの確認

まずはデータ(観察変数)の正規性を確認する。MVNパッケージのmvn( )で一気に確認できる。

> mvn(data)
$multivariateNormality
             Test        Statistic              p value Result
1 Mardia Skewness 344.905276936947 8.86473544030093e-15     NO
2 Mardia Kurtosis 2.83021589156651  0.00465166039260279     NO
3             MVN             <NA>                 <NA>     NO

$univariateNormality
          Test  Variable Statistic   p value Normality
1 Shapiro-Wilk    x1        0.9928  0.1582      YES   
2 Shapiro-Wilk    x2        0.9697  <0.001      NO    
3 Shapiro-Wilk    x3        0.9523  <0.001      NO    
~~~~~(省略)~~~~~

$Descriptives
     n     Mean  Std.Dev   Median       Min       Max     25th     75th       Skew    Kurtosis
x1 301 4.935770 1.167432 5.000000 0.6666667  8.500000 4.166667 5.666667 -0.2543455  0.30753382
x2 301 6.088040 1.177451 6.000000 2.2500000  9.250000 5.250000 6.750000  0.4700766  0.33239397
x3 301 2.250415 1.130979 2.125000 0.2500000  4.500000 1.375000 3.125000  0.3834294 -0.90752645
~~~~~(省略)~~~~~

$multivariateNormalityで多変量正規分布と見なせるかどうかを検定している。結果は"NO"なので正規分布とは見なせない。 多変量正規分布と見なせるならば最尤法を使うが、見なせないならロバスト最尤法を用いる。推定方法は、モデルの当てはめのときに指定する。

2. 因子構造の指定

モデル式を書いて因子構造を指定する。SEMで指定するような複雑なモデルを記述できるように構文が用意されている。4種類の式のタイプがあるみたいですが、詳細はまた勉強します。

  • =~:潜在変数の定義
  • ~:回帰
  • ~~:残差の分散・共分散
  • y ~ 1:切片

ここでは3つの潜在因子F1F3に対して、観察変数を以下のように割り当てた。

model = "
F1 =~ x1 + x2 + x3
F2 =~ x4 + x5 + x6
F3 =~ x7 + x8 + x9
"

図1:想定したモデル(3因子9項目)

ちなみに因子間に独立を想定する場合は、以下のようにすると因子間相関係数=0になる。

model = "
F1 =~ x1 + x2 + x3
F2 =~ x4 + x5 + x6
F3 =~ x7 + x8 + x9
F1 ~~ 0*F2
F1 ~~ 0*F3
F2 ~~ 0*F3
"

さらに、

model = "
F1 =~ x1 + x2 + x3
F2 =~ x4 + x5 + x6
F3 =~ x7 + x8 + x9
F1 ~~ 1*F1
F2 ~~ 1*F2
F3 ~~ 1*F3
"

とすると、各因子の分散を1に固定する(cfa( )でstd.lv = TRUEを指定するのと同じだと思う)。

3. モデルの当てはめ

モデル式とデータをcfa( )に渡して、因子負荷や相関を推定する。 指定する主な引数は、

  • estimator:推定方法。最尤法なら"ML"、ロバスト最尤法なら"MLR"を指定。デフォルトは"ML"みたい。
  • std.lv:TRUEを指定すると因子分散が1となるように標準化される。TRUEを指定しないとデフォルトのauto.fix.first = TRUEが指定され、各因子の中で1つ目に指定した変数のパラメータが1に固定される。因子分散を1に固定すると因子のバラツキ具合について比較できないが、1に固定して表示するのが普通なのかな。
fit <- cfa(model = model, data = data, estimator = "MLR", std.lv = TRUE)

4. 推定結果の確認

結果は以下のコードで表示。

  • fit.measures:適合度指標を表示したいときはTRUEにする
  • standardized:結果でStd.lvStd.allを表示したいときはTRUEにする。
summary(fit, fit.measures=FALSE, standardized=TRUE)

結果は長いので分割しながら見ていく。

4-1. 推定の概要

使用されたlavaanのバージョン、収束までの回数、パラメータ数、観測数などが表示される。

lavaan 0.6-8 ended normally after 20 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                        21
                                                      
  Number of observations                           301  

"ended normally"とあるので、ここでは推定は正常に収束したみたい。

4-2. モデルの適合度

ここではRobustの方を見る。

Model Test User Model:
                                               Standard      Robust
  Test Statistic                                 85.306      87.132
  Degrees of freedom                                 24          24
  P-value (Chi-square)                            0.000       0.000 #
  Scaling correction factor                                   0.979
       Yuan-Bentler correction (Mplus variant)                     

Model Test Baseline Model:

  Test statistic                               918.852     880.082
  Degrees of freedom                                36          36
  P-value                                        0.000       0.000
  Scaling correction factor                                  1.044

User Model versus Baseline Model:

  Comparative Fit Index (CFI)                    0.931       0.925 #
  Tucker-Lewis Index (TLI)                       0.896       0.888 #
                                                                  
  Robust Comparative Fit Index (CFI)                         0.930 #
  Robust Tucker-Lewis Index (TLI)                            0.895 #

Loglikelihood and Information Criteria:

  Loglikelihood user model (H0)              -3737.745   -3737.745
  Scaling correction factor                                  1.133
      for the MLR correction                                      
  Loglikelihood unrestricted model (H1)      -3695.092   -3695.092
  Scaling correction factor                                  1.051
      for the MLR correction                                      
                                                                  
  Akaike (AIC)                                7517.490    7517.490
  Bayesian (BIC)                              7595.339    7595.339
  Sample-size adjusted Bayesian (BIC)         7528.739    7528.739

Root Mean Square Error of Approximation:

  RMSEA                                          0.092       0.093 #
  90 Percent confidence interval - lower         0.071       0.073
  90 Percent confidence interval - upper         0.114       0.115
  P-value RMSEA <= 0.05                          0.001       0.001 
                                                                  
  Robust RMSEA                                               0.092 #
  90 Percent confidence interval - lower                     0.072
  90 Percent confidence interval - upper                     0.114

Standardized Root Mean Square Residual:

  SRMR                                           0.065       0.065 #
  • P-value (Chi-square):最尤統計量から得たカイ2乗分布に基づくP値。Model Test User Modelのところを見る。>0.05ならOKと言われるが、サンプルサイズが大きいと感度が良すぎる(75-200がちょうどよく、400以上ではいつも有意になってしまう*2)。この「感度が良すぎる問題」があるので、普通は以下の指標が使われる。
  • 相対的な適合との指標:ベースラインモデル(最も簡単なモデル)との相対的な比較を行う。Model Test User ModelModel Test Baseline Modelの情報を使っている。
    • Comparative Fit Index (CFI) :0〜1の値を取る。大きいほど良い。>0.95(甘くしても>0.9)ならOK。多分ここでもRobust CFIの方を見るのかな、と。
    • Tucker-Lewis Index (TLI):0〜1の値を取る(1以上になるときは1に丸める)。>0.90なら適合度OK。
  • 絶対的な適合度の指標:
    • SRMR:Standard root mean squareで、≦0.08ならOK
    • RMSEA:root mean square error of approximationの略。RMSEAが大きいほど、適合度が悪い。≦0.05ならぴったり適合、≧0.10なら適合良くない。ロバスト最尤法を使ったのなら、Robustの方を見るのかな。

4-3. 推定値

Parameter Estimates:

  Standard errors                             Sandwich
  Information bread                           Observed
  Observed information based on                Hessian

Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  F1 =~                                                                 
    x1                0.900    0.100    8.973    0.000    0.900    0.772
    x2                0.498    0.088    5.681    0.000    0.498    0.424
    x3                0.656    0.080    8.151    0.000    0.656    0.581
  F2 =~                                                                 
    x4                0.990    0.061   16.150    0.000    0.990    0.852
    x5                1.102    0.055   20.146    0.000    1.102    0.855
    x6                0.917    0.058   15.767    0.000    0.917    0.838
  F3 =~                                                                 
    x7                0.619    0.086    7.193    0.000    0.619    0.570
    x8                0.731    0.093    7.875    0.000    0.731    0.723
    x9                0.670    0.099    6.761    0.000    0.670    0.665

Covariances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  F1 ~~                                                                 
    F2                0.459    0.073    6.258    0.000    0.459    0.459
    F3                0.471    0.119    3.954    0.000    0.471    0.471
  F2 ~~                                                                 
    F3                0.283    0.085    3.311    0.001    0.283    0.283

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .x1                0.549    0.156    3.509    0.000    0.549    0.404
   .x2                1.134    0.112   10.135    0.000    1.134    0.821
   .x3                0.844    0.100    8.419    0.000    0.844    0.662
   .x4                0.371    0.050    7.382    0.000    0.371    0.275
   .x5                0.446    0.057    7.870    0.000    0.446    0.269
   .x6                0.356    0.047    7.658    0.000    0.356    0.298
   .x7                0.799    0.097    8.222    0.000    0.799    0.676
   .x8                0.488    0.120    4.080    0.000    0.488    0.477
   .x9                0.566    0.119    4.768    0.000    0.566    0.558
    F1                1.000                               1.000    1.000
    F2                1.000                               1.000    1.000
    F3                1.000                               1.000    1.000

Std.lvは潜在変数だけが標準化されていて、Std.allは潜在変数と観測変数の両方が標準化されている。

因子負荷量(factor loading, FL)

  • Latent VariableStd.allを見る
  • FLは0〜1の範囲に収まっているはずで、>1の場合はHeywood caseと呼ばれる*3
  • 0.5なら意味のある関連と考えていい

  • 全ての因子負荷量のP値は統計学的に有意でなくてはならない。

因子間相関(factor correlation)

  • CovariancesStd.allを見る
  • <0.85であれば因子が明確に区別されていることを示す。反対に>0.85だと多重共線性の問題ありかも。

因子分散および独自因子の分散

  • VariancesEstimateを見る
  • ここでは因子分散を1に固定しているので上記の結果。

4-4. パス図を描く

semPlotパッケージのsemPaths( )を使ってパス図を描く。

第1引数はcfa( )で作成したオブジェクト。その後、"path", "std"を指定。前者だけだと矢印しか描かれない。 styleは分散を表現するときに、双頭矢印(両方とも同じノードに向かう)なら"ram", "mx", "OpenMx"、一方向矢印なら"lisrel"を指定する。 デフォルトだと矢印の色がグレーなので、黒にしたいときはedge.colorで指定。

semPaths(fit, "path", "std", 
         style = "lisrel", edge.color = "black", intercepts = FALSE)

図2:semPaths( )を使ってCFAのパス図を描いた

先程のパラメータ推定値のStd.allが書き込まれていることを確認。

5. モデルの改善

5-1. 当てはまりの悪い部分を探す

モデルの当てはまりが悪かった場合、どこを変えたらいいのかを見るために、当てはまりの悪い部分を探す。

1つ目はmodification index(MI)で、>3.84を超えるところは改善の余地ありと判断。

> mis_fit <-  modificationIndices(fit)
> subset(mis_fit, mi>3.84)
   lhs op rhs     mi    epc sepc.lv sepc.all sepc.nox
10  F1 ~~  F2 26.552  0.388   0.388    0.388    0.388
11  F1 ~~  F3 22.633  0.393   0.393    0.393    0.393
12  F2 ~~  F3 12.152  0.248   0.248    0.248    0.248
25  F1 =~  x4  4.653  0.112   0.112    0.097    0.097
27  F1 =~  x6  9.912  0.157   0.157    0.144    0.144
28  F1 =~  x7  4.725 -0.145  -0.145   -0.133   -0.133
30  F1 =~  x9 40.044  0.399   0.399    0.396    0.396
31  F2 =~  x1 31.475  0.354   0.354    0.304    0.304
36  F2 =~  x9 10.298  0.175   0.175    0.173    0.173
37  F3 =~  x1 10.504  0.224   0.224    0.192    0.192
39  F3 =~  x3  4.612  0.142   0.142    0.125    0.125
45  x1 ~~  x4  6.547  0.111   0.111    0.197    0.197
50  x1 ~~  x9 14.396  0.200   0.200    0.263    0.263
55  x2 ~~  x7  6.279 -0.148  -0.148   -0.166   -0.166
59  x3 ~~  x5  4.487 -0.095  -0.095   -0.186   -0.186
63  x3 ~~  x9  4.651  0.109   0.109    0.164    0.164
66  x4 ~~  x7  7.095  0.106   0.106    0.198    0.198

もう1つはstandardized residuals(SR)で、|SR|>2.58は改善の余地ありと判断。

> sr <- residuals(fit, type = "standardized")
> sr
$type
[1] "standardized"

$cov
   x1     x2     x3     x4     x5     x6     x7     x8     x9    
x1  0.000                                                        
x2  0.000  0.000                                                 
x3  0.000  0.000  0.000                                          
x4  5.497  2.517  2.784  0.000                                   
x5  4.688  2.394  1.309  0.000  0.000                            
x6  5.119  2.996  3.126  0.000  0.000  0.000                     
x7  1.084 -1.430  1.247  2.928  1.662  2.031  0.000              
x8  3.531  1.682  3.081  1.798  2.228  2.358  0.000  0.000       
x9  6.325  3.613  5.462  3.138  3.849  3.025  0.000  0.000  0.000

5-2. 因子の信頼性を評価する

semToolsパッケージのreliablity( )でRaykov’s rhoを計算する。

> rel <- reliability(fit) 
> print(rel, digits = 3)
          F1    F2    F3
alpha  0.626 0.883 0.688
omega  0.633 0.886 0.696
omega2 0.633 0.886 0.696
omega3 0.633 0.886 0.696
avevar 0.369 0.723 0.438

omegaがRaykov’s rhoに相当。≧0.70なら信頼性あり。 alphaはCronbachのα係数。

おわりに

  • モデル式を使ってパラメータとか分散を自在に固定できるようになるには、まだまだ訓練が必要だが...私には使う機会がなさそう。
  • しばらくウェットの餌が続くと、ドライを出したときに「これじゃないでしょ?」って顔で見てきます。

参考資料

  • lavaanパッケージのチュートリアルです。

lavaan.ugent.be

  • 適合度指標の計算方法については、以下の資料を参考にしました。

stats.oarc.ucla.edu