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年の産業界における労働力の割合
分析モデル
| 記号 | 意味 |
|---|---|
| =~ | 測定方程式 |
| ~ | 構造方程式 |
| ~~ | 共変関係 |
まずはモデルを組み立てる。
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
