lavaanパッケージで共分散構造分析を行う[R]

lavaanパッケージ習得のためのお勉強ノート。

http://lavaan.ugent.be/tutorial/sem.htmllavaan.ugent.be

library(lavaan)

データの詳細

lavaanに同梱されている選挙に関するデータを使用する。Bollenが1989年に書いた本で使われていたデータ。

11変数の75個の観測からなるデータフレーム。

  • y1 1960年の報道の自由に関する専門家の評価
  • y2 1960年の政治的反対運動の自由
  • y3 1960年の選挙の公平性
  • y4 1960年の選出立法府の有効性
  • y5 1965年の報道の自由に関する専門家の評価
  • y6 1965年の政治的反対運動の自由
  • y7 1965年の選挙の公平性について
  • Y8 1965年の選出立法府の有効性
  • x1 1960年の一人当たり国民総生産(GNP)
  • x2 1960年の一人当たりの無生物エネルギー消費量
  • x3 1960年の産業界における労働力の割合

分析モデル

f:id:iDES:20200522224331p:plain

記号 意味
=~ 測定方程式
~ 構造方程式
~~ 共変関係

まずはモデルを組み立てる。

model <- '
  # measurement model
    ind60 =~ x1 + x2 + x3
    dem60 =~ y1 + y2 + y3 + y4
    dem65 =~ y5 + y6 + y7 + y8
    
  # regressions
    dem60 ~ ind60
    dem65 ~ ind60 + dem60
    
  # residual correlations
    y1 ~~ y5
    y2 ~~ y4 + y6
    y3 ~~ y7
    y4 ~~ y8
    y6 ~~ y8
'

組んだモデルをsemコマンドで分析し、summaryで結果表示させる。standardized=TRUEとオプションを指定しておくと標準化した値を出力してくれる。

fit <- sem(model=model, data=PoliticalDemocracy, estimator="ML")
summary(fit, standardized=TRUE) # 標準化

結果。

lavaan 0.6-5 ended normally after 68 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of free parameters                         31
                                                      
  Number of observations                            75
                                                      
Model Test User Model:
                                                      
  Test statistic                                38.125
  Degrees of freedom                                35
  P-value (Chi-square)                           0.329

Parameter Estimates:

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

Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  ind60 =~                                                              
    x1                1.000                               0.670    0.920
    x2                2.180    0.139   15.742    0.000    1.460    0.973
    x3                1.819    0.152   11.967    0.000    1.218    0.872
  dem60 =~                                                              
    y1                1.000                               2.223    0.850
    y2                1.257    0.182    6.889    0.000    2.794    0.717
    y3                1.058    0.151    6.987    0.000    2.351    0.722
    y4                1.265    0.145    8.722    0.000    2.812    0.846
  dem65 =~                                                              
    y5                1.000                               2.103    0.808
    y6                1.186    0.169    7.024    0.000    2.493    0.746
    y7                1.280    0.160    8.002    0.000    2.691    0.824
    y8                1.266    0.158    8.007    0.000    2.662    0.828

Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  dem60 ~                                                               
    ind60             1.483    0.399    3.715    0.000    0.447    0.447
  dem65 ~                                                               
    ind60             0.572    0.221    2.586    0.010    0.182    0.182
    dem60             0.837    0.098    8.514    0.000    0.885    0.885

Covariances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
 .y1 ~~                                                                 
   .y5                0.624    0.358    1.741    0.082    0.624    0.296
 .y2 ~~                                                                 
   .y4                1.313    0.702    1.871    0.061    1.313    0.273
   .y6                2.153    0.734    2.934    0.003    2.153    0.356
 .y3 ~~                                                                 
   .y7                0.795    0.608    1.308    0.191    0.795    0.191
 .y4 ~~                                                                 
   .y8                0.348    0.442    0.787    0.431    0.348    0.109
 .y6 ~~                                                                 
   .y8                1.356    0.568    2.386    0.017    1.356    0.338

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .x1                0.082    0.019    4.184    0.000    0.082    0.154
   .x2                0.120    0.070    1.718    0.086    0.120    0.053
   .x3                0.467    0.090    5.177    0.000    0.467    0.239
   .y1                1.891    0.444    4.256    0.000    1.891    0.277
   .y2                7.373    1.374    5.366    0.000    7.373    0.486
   .y3                5.067    0.952    5.324    0.000    5.067    0.478
   .y4                3.148    0.739    4.261    0.000    3.148    0.285
   .y5                2.351    0.480    4.895    0.000    2.351    0.347
   .y6                4.954    0.914    5.419    0.000    4.954    0.443
   .y7                3.431    0.713    4.814    0.000    3.431    0.322
   .y8                3.254    0.695    4.685    0.000    3.254    0.315
    ind60             0.448    0.087    5.173    0.000    1.000    1.000
   .dem60             3.956    0.921    4.295    0.000    0.800    0.800
   .dem65             0.172    0.215    0.803    0.422    0.039    0.039

統計量を出力

summary(fit, fit.measure=TRUE, standardized=TRUE)

fit.measure=TRUEをつけると統計量も出力できる。統計量だけ知りたい場合はfitMeasuresコマンドを使う。

fitMeasures(fit)

結果。

               npar                fmin               chisq                  df 
             31.000               0.254              38.125              35.000 
             pvalue      baseline.chisq         baseline.df     baseline.pvalue 
              0.329             730.654              55.000               0.000 
                cfi                 tli                nnfi                 rfi 
              0.995               0.993               0.993               0.918 
                nfi                pnfi                 ifi                 rni 
              0.948               0.603               0.996               0.995 
               logl   unrestricted.logl                 aic                 bic 
          -1547.791           -1528.728            3157.582            3229.424 
             ntotal                bic2               rmsea      rmsea.ci.lower 
             75.000            3131.720               0.035               0.000 
     rmsea.ci.upper        rmsea.pvalue                 rmr          rmr_nomean 
              0.092               0.611               0.276               0.276 
               srmr        srmr_bentler srmr_bentler_nomean                crmr 
              0.044               0.044               0.044               0.049 
        crmr_nomean          srmr_mplus   srmr_mplus_nomean               cn_05 
              0.049               0.045               0.045              98.970 
              cn_01                 gfi                agfi                pgfi 
            113.803               0.923               0.854               0.489 
                mfi                ecvi 
              0.979               1.335 

cfiのみ出力する場合は下記。

fitMeasures(fit, "cfi")

必要な統計量だけ出力する場合にはC()で囲う。

fitMeasures(fit, c("chisq", "df", "pvalue", "cfi", "rmsea"))

matrixで出力する場合にはoutput = "matrix"とオプションをつける。デフォルトはvector

fitMeasures(fit, c("chisq", "df", "pvalue", "cfi", "rmsea"),output = "matrix")

結果。

chisq  38.125
df     35.000
pvalue  0.329
cfi     0.995
rmsea   0.035

vectorよりもmatrixの方があとあと加工がしやすい気がする。

χ2検定は観測対象数の影響を受けやすく,観測対象数が大きくなると棄却されやすくなってしまいます。χ2値だけを利用してモデル適合の検討を行わないようにしましょう。
CFI(comparative fit index:比較適合度指標)は,その値が0から1までの範囲に収まる指標です。1に近いほど適合が良く,一般には, 0.95以上あればよいモデルであると判断します。TLI(Tucker-Lewis Index)も1に近いほど適合が良いと判断する指標ですが,CFIと異なり, 1を超える場合もあります。
RMSEAは頻繁に利用される指標で, 0.05以下であれば当てはまりが良く, 0.1以上であれば当てはまりが悪いと判断します。また,SRMR(standardized root mean square residual :標準化残差平方平均)は,標本の分散・共分散とモデルにより再現される分散・共分散の差である残差によってモデル適合を検討するための指標です。下限値は0で, 0に近いほど適合していると判断します。(豊田秀樹『共分散構造分析[R編]』p.34)

パラメーターを出力

共分散のみ単純に表示させる場合はcoefコマンド。

coef(object=fit)

もう少し詳しいものを出したいときはparameterEstimatesコマンドを使う。
95%区間を出力する場合にはci=T、標準化した値を出す場合にはstandardized=Tをつける。

parameterEstimates(object=fit, ci=T, standardized=T)

オプションでoutput="data.frame"を指定すると出力結果をデータフレームで返せる。こちらも加工するときに便利そうだ。出力が多いので出力値を選ぶ方法はないか調べてみたもののなかった。

とはいえ、データフレームとして出力できるのであれば、一旦格納してしまえばどうにでもない。

df1 <- parameterEstimates(object=fit, ci=T, standardized=T, output="data.frame")
df1$est

例えばestであれば、その値だけ取り出すことができる。

 [1] 1.00000000 2.18036773 1.81851130 1.00000000 1.25674650 1.05771677 1.26478654 1.00000000 1.18569630 1.27951218 1.26594698 1.48300054 0.57233619 0.83734481
[15] 0.62367107 1.31311258 2.15286127 0.79496028 0.34822604 1.35616712 0.08154935 0.11980648 0.46670258 1.89139552 7.37286854 5.06746210 3.14790480 2.35097047
[29] 4.95396775 3.43137392 3.25408501 0.44843715 3.95603311 0.17248133