今回も因子分析ですが，前回と違って 確証的因子分析(Confirmatory Factor Analysis,CFA) ,というものをやってみましょう。
前回の因子分析は，探索的因子分析(Exploratory Factor Analysis,EFA) と言われます。これは，事前に因子の数や因子構造がわかっていなかったため，手探りで＝探索的に，因子分析をやってみた，ということです。EFAでは，
1.因子数を探して自分で決めて， 2.どの項目がまとまるのかを調べて（その時に分かりにくかったら回転をかけるなどして見やすくして）， 3.まとまりを見つけては「これだ」と名前をつける，
という手順が必要でした。
今回は，事前に因子の数や構造に仮定があります。この仮定は，理論的に導かれるものであったり，先行研究から考えられるものであったりします。
いずれにせよ，この項目はこれを測定している，あの項目はこれを測定している，という事前の モデル がある場合です。
例えばBIG 5の場合はしっかりした因子構造の仮定に基づき，尺度化されています。
library(psych)
data(bfi)
help(bfi)今回はこのモデルに沿った因子分析をしてみよう，というのです。
実践にあたってはlavaanパッケージを使用します。持っていない人はインストールしておいてください。
lavaanとはLAtent VAriable ANalysis,つまり「潜在変数分析」という意味です。
library(lavaan)## This is lavaan 0.5-16
## lavaan is BETA software! Please report any bugs.このパッケージの使い方は，複数行からなるモデルとデータを与えて，分析させるという方法をとります。 やってみましょう。
model <- '
A =~ A1 + A2 + A3 + A4 + A5
E =~ E1 + E2 + E3 + E4 + E5
O =~ O1 + O2 + O3 + O4 + O5
N =~ N1 + N2 + N3 + N4 + N5
C =~ C1 + C2 + C3 + C4 + C5
'
result.fit <- cfa(model,data=bfi)
summary(result.fit,standardized=T)## lavaan (0.5-16) converged normally after 55 iterations
##
## Used Total
## Number of observations 2436 2800
##
## Estimator ML
## Minimum Function Test Statistic 4165.467
## Degrees of freedom 265
## P-value (Chi-square) 0.000
##
## Parameter estimates:
##
## Information Expected
## Standard Errors Standard
##
## Estimate Std.err Z-value P(>|z|) Std.lv Std.all
## Latent variables:
## A =~
## A1 1.000 0.484 0.344
## A2 -1.579 0.108 -14.650 0.000 -0.764 -0.648
## A3 -2.030 0.134 -15.093 0.000 -0.983 -0.749
## A4 -1.564 0.115 -13.616 0.000 -0.757 -0.510
## A5 -1.804 0.121 -14.852 0.000 -0.873 -0.687
## E =~
## E1 1.000 0.920 0.564
## E2 1.226 0.051 23.899 0.000 1.128 0.699
## E3 -0.921 0.041 -22.431 0.000 -0.847 -0.627
## E4 -1.121 0.047 -23.977 0.000 -1.031 -0.703
## E5 -0.808 0.039 -20.648 0.000 -0.743 -0.553
## O =~
## O1 1.000 0.635 0.564
## O2 -1.020 0.068 -14.962 0.000 -0.648 -0.418
## O3 1.373 0.072 18.942 0.000 0.872 0.724
## O4 0.437 0.048 9.160 0.000 0.277 0.233
## O5 -0.960 0.060 -16.056 0.000 -0.610 -0.461
## N =~
## N1 1.000 1.300 0.825
## N2 0.947 0.024 39.899 0.000 1.230 0.803
## N3 0.884 0.025 35.919 0.000 1.149 0.721
## N4 0.692 0.025 27.753 0.000 0.899 0.573
## N5 0.628 0.026 24.027 0.000 0.816 0.503
## C =~
## C1 1.000 0.680 0.551
## C2 1.148 0.057 20.152 0.000 0.781 0.592
## C3 1.036 0.054 19.172 0.000 0.705 0.546
## C4 -1.421 0.065 -21.924 0.000 -0.967 -0.702
## C5 -1.489 0.072 -20.694 0.000 -1.012 -0.620
##
## Covariances:
## A ~~
## E 0.304 0.025 12.293 0.000 0.683 0.683
## O -0.093 0.011 -8.446 0.000 -0.303 -0.303
## N 0.141 0.018 7.712 0.000 0.223 0.223
## C -0.110 0.012 -9.254 0.000 -0.334 -0.334
## E ~~
## O -0.265 0.021 -12.347 0.000 -0.453 -0.453
## N 0.292 0.032 9.131 0.000 0.244 0.244
## C -0.224 0.020 -11.121 0.000 -0.357 -0.357
## O ~~
## N -0.093 0.022 -4.138 0.000 -0.112 -0.112
## C 0.130 0.014 9.190 0.000 0.301 0.301
## N ~~
## C -0.250 0.025 -10.117 0.000 -0.283 -0.283
##
## Variances:
## A1 1.745 0.052 1.745 0.882
## A2 0.807 0.028 0.807 0.580
## A3 0.754 0.032 0.754 0.438
## A4 1.632 0.051 1.632 0.740
## A5 0.852 0.032 0.852 0.528
## E1 1.814 0.058 1.814 0.682
## E2 1.332 0.049 1.332 0.512
## E3 1.108 0.038 1.108 0.607
## E4 1.088 0.041 1.088 0.506
## E5 1.251 0.040 1.251 0.694
## O1 0.865 0.032 0.865 0.682
## O2 1.990 0.063 1.990 0.826
## O3 0.691 0.039 0.691 0.476
## O4 1.346 0.040 1.346 0.946
## O5 1.380 0.045 1.380 0.788
## N1 0.793 0.037 0.793 0.320
## N2 0.836 0.036 0.836 0.356
## N3 1.222 0.043 1.222 0.481
## N4 1.654 0.052 1.654 0.672
## N5 1.969 0.060 1.969 0.747
## C1 1.063 0.035 1.063 0.697
## C2 1.130 0.039 1.130 0.650
## C3 1.170 0.039 1.170 0.702
## C4 0.960 0.040 0.960 0.507
## C5 1.640 0.059 1.640 0.615
## A 0.234 0.030 1.000 1.000
## E 0.846 0.062 1.000 1.000
## O 0.404 0.033 1.000 1.000
## N 1.689 0.073 1.000 1.000
## C 0.463 0.036 1.000 1.000多くの結果が出てきます。読み方に注意しつつ，探索的因子分析の結果との共通点・相違点をよく確認しておきましょう。
事前に仮説＝モデルをもってデータに当てはめるわけですから，仮説が間違っていたらうまくいきません。 仮説とデータの相性が良かったのかどうか，これを確認するのが 適合度 です。 適合度は数字として出てきます。複数の適合度指標があり，総合的に判断します。
fitMeasures(result.fit)## fmin chisq df
## 0.855 4165.467 265.000
## pvalue baseline.chisq baseline.df
## 0.000 18222.116 300.000
## baseline.pvalue cfi tli
## 0.000 0.782 0.754
## nnfi rfi nfi
## 0.754 0.741 0.771
## pnfi ifi rni
## 0.681 0.783 0.782
## logl unrestricted.logl npar
## -99840.238 -97757.504 60.000
## aic bic ntotal
## 199800.476 200148.363 2436.000
## bic2 rmsea rmsea.ci.lower
## 199957.729 0.078 0.076
## rmsea.ci.upper rmsea.pvalue rmr
## 0.080 0.000 0.157
## rmr_nomean srmr srmr_bentler
## 0.157 0.075 0.075
## srmr_bentler_nomean srmr_bollen srmr_bollen_nomean
## 0.075 0.075 0.075
## srmr_mplus srmr_mplus_nomean cn_05
## 0.075 0.075 178.764
## cn_01 gfi agfi
## 189.003 0.862 0.830
## pgfi mfi ecvi
## 0.703 0.449 1.759有名なのは，GFI.AGFI,RMSEA,CFI,AICなどです。 それぞれを個別に呼び出すこともできます。
fitMeasures(result.fit,"AGFI")## agfi
## 0.83fitMeasures(result.fit,c("GFI","RMSEA"))## gfi rmsea
## 0.862 0.078あるいは，結果と合わせて出力することもできます。こちらの方が一般的で見やすいですね。
summary(result.fit,standardized=T,fit.measures=T)## lavaan (0.5-16) converged normally after 55 iterations
##
## Used Total
## Number of observations 2436 2800
##
## Estimator ML
## Minimum Function Test Statistic 4165.467
## Degrees of freedom 265
## P-value (Chi-square) 0.000
##
## Model test baseline model:
##
## Minimum Function Test Statistic 18222.116
## Degrees of freedom 300
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.782
## Tucker-Lewis Index (TLI) 0.754
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -99840.238
## Loglikelihood unrestricted model (H1) -97757.504
##
## Number of free parameters 60
## Akaike (AIC) 199800.476
## Bayesian (BIC) 200148.363
## Sample-size adjusted Bayesian (BIC) 199957.729
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.078
## 90 Percent Confidence Interval 0.076 0.080
## P-value RMSEA <= 0.05 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.075
##
## Parameter estimates:
##
## Information Expected
## Standard Errors Standard
##
## Estimate Std.err Z-value P(>|z|) Std.lv Std.all
## Latent variables:
## A =~
## A1 1.000 0.484 0.344
## A2 -1.579 0.108 -14.650 0.000 -0.764 -0.648
## A3 -2.030 0.134 -15.093 0.000 -0.983 -0.749
## A4 -1.564 0.115 -13.616 0.000 -0.757 -0.510
## A5 -1.804 0.121 -14.852 0.000 -0.873 -0.687
## E =~
## E1 1.000 0.920 0.564
## E2 1.226 0.051 23.899 0.000 1.128 0.699
## E3 -0.921 0.041 -22.431 0.000 -0.847 -0.627
## E4 -1.121 0.047 -23.977 0.000 -1.031 -0.703
## E5 -0.808 0.039 -20.648 0.000 -0.743 -0.553
## O =~
## O1 1.000 0.635 0.564
## O2 -1.020 0.068 -14.962 0.000 -0.648 -0.418
## O3 1.373 0.072 18.942 0.000 0.872 0.724
## O4 0.437 0.048 9.160 0.000 0.277 0.233
## O5 -0.960 0.060 -16.056 0.000 -0.610 -0.461
## N =~
## N1 1.000 1.300 0.825
## N2 0.947 0.024 39.899 0.000 1.230 0.803
## N3 0.884 0.025 35.919 0.000 1.149 0.721
## N4 0.692 0.025 27.753 0.000 0.899 0.573
## N5 0.628 0.026 24.027 0.000 0.816 0.503
## C =~
## C1 1.000 0.680 0.551
## C2 1.148 0.057 20.152 0.000 0.781 0.592
## C3 1.036 0.054 19.172 0.000 0.705 0.546
## C4 -1.421 0.065 -21.924 0.000 -0.967 -0.702
## C5 -1.489 0.072 -20.694 0.000 -1.012 -0.620
##
## Covariances:
## A ~~
## E 0.304 0.025 12.293 0.000 0.683 0.683
## O -0.093 0.011 -8.446 0.000 -0.303 -0.303
## N 0.141 0.018 7.712 0.000 0.223 0.223
## C -0.110 0.012 -9.254 0.000 -0.334 -0.334
## E ~~
## O -0.265 0.021 -12.347 0.000 -0.453 -0.453
## N 0.292 0.032 9.131 0.000 0.244 0.244
## C -0.224 0.020 -11.121 0.000 -0.357 -0.357
## O ~~
## N -0.093 0.022 -4.138 0.000 -0.112 -0.112
## C 0.130 0.014 9.190 0.000 0.301 0.301
## N ~~
## C -0.250 0.025 -10.117 0.000 -0.283 -0.283
##
## Variances:
## A1 1.745 0.052 1.745 0.882
## A2 0.807 0.028 0.807 0.580
## A3 0.754 0.032 0.754 0.438
## A4 1.632 0.051 1.632 0.740
## A5 0.852 0.032 0.852 0.528
## E1 1.814 0.058 1.814 0.682
## E2 1.332 0.049 1.332 0.512
## E3 1.108 0.038 1.108 0.607
## E4 1.088 0.041 1.088 0.506
## E5 1.251 0.040 1.251 0.694
## O1 0.865 0.032 0.865 0.682
## O2 1.990 0.063 1.990 0.826
## O3 0.691 0.039 0.691 0.476
## O4 1.346 0.040 1.346 0.946
## O5 1.380 0.045 1.380 0.788
## N1 0.793 0.037 0.793 0.320
## N2 0.836 0.036 0.836 0.356
## N3 1.222 0.043 1.222 0.481
## N4 1.654 0.052 1.654 0.672
## N5 1.969 0.060 1.969 0.747
## C1 1.063 0.035 1.063 0.697
## C2 1.130 0.039 1.130 0.650
## C3 1.170 0.039 1.170 0.702
## C4 0.960 0.040 0.960 0.507
## C5 1.640 0.059 1.640 0.615
## A 0.234 0.030 1.000 1.000
## E 0.846 0.062 1.000 1.000
## O 0.404 0.033 1.000 1.000
## N 1.689 0.073 1.000 1.000
## C 0.463 0.036 1.000 1.000適合度をよくするために，モデルを修正することも少なくありません。 一般的に，GFIやCFIは0.9以上が望ましいとされていますから，これを目指してモデルを改善していくことにします。
例えば，係数の小さかったA1やO4を外し，相関が低いところをゼロにしてみます。
model2 <- '
A =~ A2 + A3 + A4 + A5
E =~ E1 + E2 + E3 + E4 + E5
O =~ O1 + O2 + O3 + O5
N =~ N1 + N2 + N3 + N4 + N5
C =~ C1 + C2 + C3 + C4 + C5
A ~~ 0*N
E ~~ 0*N
'
result.fit2 <- cfa(model2,data=bfi)
fitMeasures(result.fit,c("GFI","AGFI","TLI","CFI","AIC","BIC"))## gfi agfi tli cfi aic bic
## 8.620e-01 8.300e-01 7.540e-01 7.820e-01 1.998e+05 2.001e+05fitMeasures(result.fit2,c("GFI","AGFI","TLI","CFI","AIC","BIC"))## gfi agfi tli cfi aic bic
## 8.750e-01 8.450e-01 7.730e-01 8.010e-01 1.853e+05 1.856e+05幾分改良されたようですが，まだまだ改良できるでしょうかね。
どこを改善すればよいかを算出してくれる関数もあります。modindices関数ですが，大量になるので一部だけ表示させます。
md <- modindices(result.fit2)
head(md[order(md$mi,decreasing=T),])## lhs op rhs mi epc sepc.lv sepc.all sepc.nox
## 1 N1 ~~ N2 388.0 0.849 0.849 0.351 0.351
## 2 E =~ N4 220.6 0.489 0.448 0.285 0.285
## 3 N3 ~~ N4 156.5 0.437 0.437 0.175 0.175
## 4 O =~ E3 128.1 0.638 0.401 0.297 0.297
## 5 O =~ E4 127.9 -0.680 -0.427 -0.291 -0.291
## 6 C =~ E5 126.1 0.504 0.339 0.252 0.252これに従って，妥当なところは モデルを変更してみましょう。
model3 <- '
A =~ A2 + A3 + A4 + A5
E =~ E1 + E2 + E3 + E4 + E5
O =~ O1 + O2 + O3 + O5
N =~ N1 + N2 + N3 + N4 + N5
C =~ C1 + C2 + C3 + C4 + C5
A ~~ 0*N
E ~~ 0*N
N1 ~~ N2
N3 ~~ N4
'
result.fit3 <- cfa(model3,data=bfi)
fitMeasures(result.fit2,c("GFI","AGFI","TLI","CFI","AIC","BIC"))## gfi agfi tli cfi aic bic
## 8.750e-01 8.450e-01 7.730e-01 8.010e-01 1.853e+05 1.856e+05fitMeasures(result.fit3,c("GFI","AGFI","TLI","CFI","AIC","BIC"))## gfi agfi tli cfi aic bic
## 8.880e-01 8.590e-01 7.920e-01 8.200e-01 1.850e+05 1.853e+05こうして，より良いモデルに近づけていくのです。