いきなり因子分析（その6）：lavaanを使って確証的因子分析

シリーズものの記事。過去の記事はこちら。

その1（とりあえず試した編）
その2（因子数の決定編）
その3（様々な推定法編）
その4（因子軸の回転編）
その5（高次因子モデルと双因子モデル編）

前回までは基本的に探索的因子分析（Exploratory Factor Analysis）に焦点を当てていた。これは観測されたデータを説明するモデルを探索的に調べるときに使われる手法である。

一方で、観測データと因子との間の関係についての仮説的なモデルを立てて、そのモデルをデータに基づき評価する際に使われるのが、確証的因子分析（Confirmatory Factor Analysis）である。（あるいは確認的因子分析や検証的因子分析とも）

前回までのデータセットを用いてやってみる。

lavaanのパッケージ

Rで確認的因子分析をするときはlavaanというパッケージに含まれるcfa()を使えばできるようだ。さっそくやってみる。

まず最初にやるべきことは、モデルを表現する式を書くことらしい。式の記法は決まっており, lavaanのパッケージの説明pdfの日本語版には次のように書いてある。

モデルのタイプ	オペレータ	意味
潜在変数の定義	=∼	~により測定される
回帰	∼	~への回帰
(残差) (共)分散	∼∼	〜と相関を持つ
切片	∼1	切片

これらを使って、潜在変数・観測変数の関係性を記述するのである。

前回までの記事では、観測変数が9つあり、それらが「視空間能力」「言語」「処理速度」と3つの因子によって説明されるとしていたのであった。そこで、潜在変数を定義するオペレータを用いて次のようにモデルを書く。

library(lavaan)
dat <- HolzingerSwineford1939

model1 <- "
visual =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed =~ x7 + x8 + x9
"

モデル式の右辺はデータセットの列名である。モデルは''もしくは""で囲む。

あとはcfa()に作ったモデルとデータを読ませるだけである。

# モデルへのあてはめ
fit1 <- cfa(model = model1,
            data = dat
            )

modelという引数には書いたモデルを、dataにはデータセットを指定する。

summary(fit1)を用いると結果の要約をみることができるのでさっそく見てみよう。

> summary(fit1)
lavaan (0.6-1) converged normally after  35 iterations

  Number of observations                           301

  Estimator                                         ML
  Model Fit Test Statistic                      85.306
  Degrees of freedom                                24
  P-value (Chi-square)                           0.000

一番最初には

lavaanのバージョン
解が収束したか
何回の繰り返しが必要だったかが
分析で使用された観察数
使われた推定量
検定統計量、自由度、p値

などが表示されている。

次のブロックにはパラメータについての記述。得られたパラメータが期待値なのか観測値なのか、標準誤差が標準のものかブートストラップなのかなどが表示されるらしい。（詳しくはよく分からん）

Parameter Estimates:

  Information                                 Expected
  Information saturated (h1) model          Structured
  Standard Errors                             Standard

そこから先が各変数のパラメータの推定値とその標準誤差、Waldの統計量とそのp値が出力される。Waldの統計量とは推定値をその標準誤差で割ったもののようだ。

cfa()は各潜在変数の最初の因子負荷量を1に固定しているので、当然のことながら推定をおこなっていないそれらの標準誤差は存在しない。

Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)
  visual =~                                           
    x1                1.000                           
    x2                0.554    0.100    5.554    0.000
    x3                0.729    0.109    6.685    0.000
  textual =~                                          
    x4                1.000                           
    x5                1.113    0.065   17.014    0.000
    x6                0.926    0.055   16.703    0.000
  speed =~                                            
    x7                1.000                           
    x8                1.180    0.165    7.152    0.000
    x9                1.082    0.151    7.155    0.000

Covariances:
                   Estimate  Std.Err  z-value  P(>|z|)
  visual ~~                                           
    textual           0.408    0.074    5.552    0.000
    speed             0.262    0.056    4.660    0.000
  textual ~~                                          
    speed             0.173    0.049    3.518    0.000

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)
   .x1                0.549    0.114    4.833    0.000
   .x2                1.134    0.102   11.146    0.000
   .x3                0.844    0.091    9.317    0.000
   .x4                0.371    0.048    7.779    0.000
   .x5                0.446    0.058    7.642    0.000
   .x6                0.356    0.043    8.277    0.000
   .x7                0.799    0.081    9.823    0.000
   .x8                0.488    0.074    6.573    0.000
   .x9                0.566    0.071    8.003    0.000
    visual            0.809    0.145    5.564    0.000
    textual           0.979    0.112    8.737    0.000
    speed             0.384    0.086    4.451    0.000

ところで、モデルの当てはまりの良さを表す指標を得たい場合には次のように引数を指定する。

summary(fit1, fit.measures=TRUE)

そうすると、出力が表示される次のようなブロックが表示される。

Model test baseline model:

  Minimum Function Test Statistic              918.852
  Degrees of freedom                                36
  P-value                                        0.000

User model versus baseline model:

  Comparative Fit Index (CFI)                    0.931
  Tucker-Lewis Index (TLI)                       0.896

Loglikelihood and Information Criteria:

  Loglikelihood user model (H0)              -3737.745
  Loglikelihood unrestricted model (H1)      -3695.092

  Number of free parameters                         21
  Akaike (AIC)                                7517.490
  Bayesian (BIC)                              7595.339
  Sample-size adjusted Bayesian (BIC)         7528.739

Root Mean Square Error of Approximation:

  RMSEA                                          0.092
  90 Percent Confidence Interval          0.071  0.114
  P-value RMSEA <= 0.05                          0.001

Standardized Root Mean Square Residual:

  SRMR                                           0.065

一番最初にbeseline modelというものについての情報が表示される。これは独立モデル（independence model）とも呼ばれるもので、変数と変数の間の関連を一切想定しないモデルのことを指す。当然そんなモデルは当てはまりがとっても悪いので、そことの比較においてモデルの当てはまりを評価しようというのが基本の考え方だ。

次のブロックには、指定したモデルとベースラインモデルを比較した指標が２つ出力される。どちらも、指定されたモデルとベースラインモデルのχ二乗値を用いて計算する。CFIもTLIも1に近いほど、当てはまりがよいとされる。

その次のブロックは対数尤度と情報量基準に基づく指標。対数尤度尤度の方はちょっとよく知らないけれど、情報量基準の方は2つ以上のモデルを比較するときに使われるようである。小さいほど当てはまりが良いとされ絶対値に意味はない。

さらに次のブロックにはRMSEAとSRMRは0.1以上だと当てはまりがよくないとされるとか。

χ二乗値が計算のもとになるのだけれど、これが自由度とかサンプルサイズとかの影響を受けるから、それをうけないように指標を改良してきているようである。

モデルの比較

せっかくなので、モデルを作って比較してみる。以前の記事で探索的因子分析の際に、そういえば観測変数x9が因子2と因子3のどちらにも影響されていてどっちつかずだったので、x9の観測変数を削除したモデルと比較してみることにしよう。モデルを次のように書いて、同様に結果を出力する。

model2 <- "
visual =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed =~ x7 + x8
"

fit2 <- cfa(model = model2,
            data = dat
            )
summary(fit2, fit.measure=TRUE)

これを、最初に行ったモデル1の適合度の指標を比べてみる。

指標	モデル1	モデル2
chi²	85.306	44.409
Comparative Fit Index (CFI)	0.931	0.964
Tucker-Lewis Index(TLI)	0.896	0.941
AIC	7517.490	6742.935
Bayesian (BIC)	7595.339	6813.370
Sample-size adjusted Bayesian (BIC)	7528.739	6753.113
RMSEA	0.092	0.073
SRMR	0.065	0.052

モデルのデータへの当てはまりのみを見れば項目9を削除した方が良いことがわかる。もっとも、モデルの選択であったり項目の削除であったりは、適合度のみをもとに判断するのでなく、理論的な仮説に基づいてなされるべきであるので、モデルが想定する変数間の関係を変えて適合度を上げつつも、理論的な仮説と整合性のあるモデルを組まなければならないようだ。（だから今回のように内容を吟味することもなく数値だけで判断するのは過ちということである）