本資料は、lavaan 公式サイト(https://lavaan.ugent.be/ )の内容を主に参照しながら、 lavaanを使った分析の操作について簡潔に解説したものである。lavaanは現在も開発中であり、最新の情報については上記URLを参照されたい。また、拙著「縦断データ分析:構造方程式モデリングによるアプローチ」(宇佐美, 2026)の補足資料としても使用しており、本資料内における「本書」はこれを指す。 lavaanについての引用は、以下の通り:
lavaan は R パッケージとして実装されている。したがって、lavaan をインストールする前に、まず新しい版(4.0.0 以上)の R をインストールしておく必要がある。RStudio および R は、ここ からダウンロードできる。
R Markdown は、本資料のように、本文とR
コードを一つのファイルにまとめ、HTML / PDF / Word
に出力できる形式である。このR
Markdownファイル(lavaan_tutorial_ja.Rmd)を Knit する(=文章と
Rコードをまとめて実行・変換する)には、RStudio
に加えて、rmarkdown と knitr
を利用する必要がある。入っていない場合は、R
コンソールで一度だけ次を実行する。
options(repos = c(CRAN = "https://cloud.r-project.org"))
install.packages(c("rmarkdown", "knitr"))
このファイルを例えば lavaan_tutorial_ja.Rmd
として保存したあと、RStudio 上で Knit を押すと HTML /
PDF / Word に出力できる。
あるいは、R コンソールから次のように実行する。
rmarkdown::render("lavaan_tutorial_ja.Rmd")
HTML に限定すれば次のようにも書ける。
rmarkdown::render("lavaan_tutorial_ja.Rmd", output_format = "html_document")
lavaan
パッケージがRにインストールされていない場合は、以下からインストールする。
if (!requireNamespace("lavaan", quietly = TRUE)) {
install.packages("lavaan", dependencies = TRUE)
}
library(lavaan) #lavaanパッケージのダウンロード(Rを開いたら逐一必要)
lavaan には、以下のような学習用の組み込みデータがいくつか含まれている:
data("HolzingerSwineford1939")
data("PoliticalDemocracy")
head(HolzingerSwineford1939)
id sex ageyr agemo school grade x1 x2 x3 x4 x5 x6
1 1 1 13 1 Pasteur 7 3.333333 7.75 0.375 2.333333 5.75 1.2857143
2 2 2 13 7 Pasteur 7 5.333333 5.25 2.125 1.666667 3.00 1.2857143
3 3 2 13 1 Pasteur 7 4.500000 5.25 1.875 1.000000 1.75 0.4285714
4 4 1 13 2 Pasteur 7 5.333333 7.75 3.000 2.666667 4.50 2.4285714
5 5 2 12 2 Pasteur 7 4.833333 4.75 0.875 2.666667 4.00 2.5714286
6 6 2 14 1 Pasteur 7 5.333333 5.00 2.250 1.000000 3.00 0.8571429
x7 x8 x9
1 3.391304 5.75 6.361111
2 3.782609 6.25 7.916667
3 3.260870 3.90 4.416667
4 3.000000 5.30 4.861111
5 3.695652 6.30 5.916667
6 4.347826 6.65 7.500000
head(PoliticalDemocracy)
y1 y2 y3 y4 y5 y6 y7 y8 x1
1 2.50 0.000000 3.333333 0.000000 1.250000 0.000000 3.726360 3.333333 4.442651
2 1.25 0.000000 3.333333 0.000000 6.250000 1.100000 6.666666 0.736999 5.384495
3 7.50 8.800000 9.999998 9.199991 8.750000 8.094061 9.999998 8.211809 5.961005
4 8.90 8.800000 9.999998 9.199991 8.907948 8.127979 9.999998 4.615086 6.285998
5 10.00 3.333333 9.999998 6.666666 7.500000 3.333333 9.999998 6.666666 5.863631
6 7.50 3.333333 6.666666 6.666666 6.250000 1.100000 6.666666 0.368500 5.533389
x2 x3
1 3.637586 2.557615
2 5.062595 3.568079
3 6.255750 5.224433
4 7.567863 6.267495
5 6.818924 4.573679
6 5.135798 3.892270
lavaan は、潜在変数モデルのための free open-source の R パッケージである。パス解析、確認的因子分析(CFA)、構造方程式モデリング(SEM)、潜在成長モデル(LGM)、交差遅延モデル群、媒介モデルなど、幅広いモデルを設定・推定できる。
一方、公式の tutorial では現状、以下のような制約が示されている:
lavaanでは、モデルを文字列として記述する。
y ~ x1 + x2 + x3 + x4
これは、y を従属変数、x1〜x4
を独立変数とする回帰式である。
すなわち、記号チルダ(“~’)が回帰のオペレータであり、このオペレータの左には従属変数(y)を与え,右には独立変数(x)を“+”オペレータで区切りながら与える。以下のように、潜在変数(f1,f2,f3)を独立変数として含む場合も基本的に同様である。
y ~ f1 + f2 + x1 + x2
f1 ~ f2 + f3
f2 ~ f3 + x1 + x2
HolzingerSwineford1939
データ(2つの学校[PasteurとGrant-White]の7年生と8年生の生徒の知能検査のスコア[x1-x9])に対する3因子の
CFA を行う。
visual(視覚):x1, x2,
x3textual(言語):x4, x5,
x6speed(速度):x7, x8,
x9モデル構文を一重引用符(’)を用いて入力する。
HS.model <- '
visual =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed =~ x7 + x8 + x9
'
次のようにして,CFAモデル(HS.model)をデータ(HolzingerSwineford1939)に当てはめる。最初の引数はユーザーが設定したモデルであり、2番目の引数は分析に用いるデータセットである。モデルを当てはめた後、関数summary()を用いて、分析結果の要約情報(推定値、標準誤差、z 値、p 値、適合度指標;本書の2.4-2.6節)を表示できる。
fit.cfa <- cfa(HS.model, data = HolzingerSwineford1939)
summary(fit.cfa, fit.measures = TRUE)
lavaan 0.6-21 ended normally after 35 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 21
Number of observations 301
Model Test User Model:
Test statistic 85.306
Degrees of freedom 24
P-value (Chi-square) 0.000
Model Test Baseline Model:
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
Akaike (AIC) 7517.490
Bayesian (BIC) 7595.339
Sample-size adjusted Bayesian (SABIC) 7528.739
Root Mean Square Error of Approximation:
RMSEA 0.092
90 Percent confidence interval - lower 0.071
90 Percent confidence interval - upper 0.114
P-value H_0: RMSEA <= 0.050 0.001
P-value H_0: RMSEA >= 0.080 0.840
Standardized Root Mean Square Residual:
SRMR 0.065
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
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
lavaan の出力は、大きく 3つの部分に分けて理解できる。第1に、ヘッダー部分である。 ここには、lavaan のバージョン、最適化が正常終了したかどうか、反復回数、使用した推定法(例:ML)、最適化アルゴリズム(例:NLMINB)、分析モデル内の母数の数、分析に実際に用いられたサンプルサイズなどが表示される。さらに、Model Test User Model: の箇所では、分析モデルに対する検定統計量、自由度、p値 が示される。
第2に、適合度指標の部分である。これは fit.measures = TRUE を付けたときに表示され、ベースラインモデルの検定結果、CFI、TLI、対数尤度、AIC、BIC、SABIC、RMSEA(とその信頼区間)、SRMRなどが出力される。したがって、この部分を見ることで、分析モデルがデータにどの程度適合しているかを総合的に確認できる。
第3に、母数推定値の部分である。ここには Latent Variables:, Covariances:, Variances: などの見出しの下に、分析モデルに含まれる自由母数と固定母数が表形式で示される。通常は、まず因子負荷、次に共分散、最後に分散や残差分散が表示される。この表のうち、特に重要なのは次の4つの列である。
Estimate 各母数の推定値、または固定した値を表す列である。たとえば因子負荷、因子間共分散、分散などの大きさがここに示される。
Std.Err 各母数の推定値の標準誤差を表す列である。推定値の不確実性の大きさを示している。
z-value Wald 統計量を表す列であり、推定値を標準誤差で割った値である。母数が 0 からどれだけ離れているかを標準誤差単位で見た指標といえる。
P(>|z|) 母集団においてその母数が 0 であるという帰無仮説を検定したときの p 値を表す列である。
加えて、分散の出力では若干の注意が必要である。Variances: の欄で、観測変数名の前に .(ドット) が付いている場合、それはその変数が内生変数であり、表示されているのは 残差分散である。一方、外生的な潜在変数にはドットが付かず、その場合は潜在変数の全分散が表示される。
PoliticalDemocracy
データを用いた、共通因子間の回帰モデル(本書2.2.5節も参照)を示す。このデータは、Bollen(1989)以来広く利用されてきたデータセットでもある。
y1 1960年の報道の自由に関する専門家の評価y2 1960年の政治的反対の自由y3 1960年の選挙の公平性y4 1960年に選出された立法府の有効性y5 1965年の報道の自由に関する専門家による評価y6 1965年の政治的反対の自由y7 1965年の選挙の公平性y8 1965年に選出された立法府の有効性x1 1960年の一人当たりの国民総生産(GNP)x2 1960年の一人当たりの無生物のエネルギー消費量x3 1960年の産業における労働力の割合ここで示している変数の意味(訳出)については、こちらのサイトを参考にさせて頂いた。
model.sem <- '
ind60 =~ x1 + x2 + x3
dem60 =~ y1 + y2 + y3 + y4
dem65 =~ y5 + y6 + y7 + y8
dem60 ~ ind60
dem65 ~ ind60 + dem60
y1 ~~ y5 #独自因子間の相関
y2 ~~ y4 + y6 #独自因子間の相関
y3 ~~ y7 #独自因子間の相関
y4 ~~ y8 #独自因子間の相関
y6 ~~ y8 #独自因子間の相関
'
fit.sem <- sem(model.sem, data = PoliticalDemocracy)
summary(fit.sem, standardized = TRUE, fit.measures = TRUE)
lavaan 0.6-21 ended normally after 68 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 31
Number of observations 75
Model Test User Model:
Test statistic 38.125
Degrees of freedom 35
P-value (Chi-square) 0.329
Model Test Baseline Model:
Test statistic 730.654
Degrees of freedom 55
P-value 0.000
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.995
Tucker-Lewis Index (TLI) 0.993
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -1547.791
Loglikelihood unrestricted model (H1) -1528.728
Akaike (AIC) 3157.582
Bayesian (BIC) 3229.424
Sample-size adjusted Bayesian (SABIC) 3131.720
Root Mean Square Error of Approximation:
RMSEA 0.035
90 Percent confidence interval - lower 0.000
90 Percent confidence interval - upper 0.092
P-value H_0: RMSEA <= 0.050 0.611
P-value H_0: RMSEA >= 0.080 0.114
Standardized Root Mean Square Residual:
SRMR 0.044
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
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
この例では、
=~:因子負荷~:回帰係数~~:分散・共分散が全て含まれている。y1~~y5は、2つの観測変数の残差が相関することを意味する。これは,各変数が対応する潜在変数だけからでは説明できない何か別の共通原因があると考えられるときにしばしば仮定される。今の場合、2つの変数は同じ種類の変数であり異なる時期(1960年と1965年)に測定されていることを反映している。y2~~y4+y6は、2つの表現y2~~y4とy2~~y6を結合したものであり、このような略記が可能である。
そして、引数standardized=TRUEは、標準化解についての出力を追加する.そのため、対応する2つの列が結果に追加されている(5,6列目)。第5列(Std.lv)は潜在変数だけが標準化され、第6列(Std.all)では潜在変数と観測変数がともに標準化された場合の結果である。
母数の固定、初期値、ラベル、制約の付与について説明する。
NA*HS.model2 <- '
visual =~ x1 + x2 + x3 #識別の観点から、最初の変数(x1)の負荷が1に固定される。
textual =~ x4 + x5 + x6 #識別の観点から、最初の変数(x4)の負荷が1に固定される。
speed =~ NA*x7 + x8 + x9 #最初の変数(x7)の負荷が1に固定されず自由推定される。
visual ~~ 0*speed #無相関を仮定
textual ~~ 0*speed #無相関を仮定
speed ~~ 1*speed #(識別のため、)speedの共通因子の分散が1であることを仮定。
'
fit.cfa2 <- cfa(HS.model2, data = HolzingerSwineford1939)
summary(fit.cfa2)
orthogonal = TRUE (直交解)と
std.lv = TRUE(標準化解)fit.ortho <- cfa(HS.model, data = HolzingerSwineford1939, orthogonal = TRUE)
fit.stdlv <- cfa(HS.model, data = HolzingerSwineford1939, std.lv = TRUE)
start()を用いて、()内の数値を初期値として設定できる。
model.start <- '
visual =~ x1 + start(0.8)*x2 + start(1.2)*x3
textual =~ x4 + start(0.5)*x5 + start(1.0)*x6
speed =~ x7 + start(0.7)*x8 + start(1.8)*x9
'
母数について名称を付与することで出力時に参照しやすくなる。また、後述する等値制約の際にもラベリングが利用できる。以下は、x3に対応する因子負荷の名称をmyLabelとしている。このように、ラベルは記号a-z, A-Zで始める必要がある。
model.label <- '
ind60 =~ x1 + x2 + myLabel*x3
dem60 =~ y1 + y2 + y3 + y4
dem65 =~ y5 + y6 + y7 + y8
dem60 ~ ind60
dem65 ~ ind60 + dem60
y1 ~~ y5
y2 ~~ y4 + y6
y3 ~~ y7
y4 ~~ y8
y6 ~~ y8
'
潜在成長モデル(latent growth model; 本書の4章-6章)は縦断データの分析に用いられ、測定値の時間的変化(軌跡)やその個人差を一次や二次などの関数によって記述することを目的とする。成長因子(growth factor)と呼ばれる共通因子を介して軌跡の個人差を記述できる。下記の例では,Demo.growthという人工データを利用している。これは、読解力尺度に関する標準化された得点のデータで、4時点で測定されている。このデータに対して線形LGM(本書の4章;切片因子iと傾き因子sからなる2つの成長因子を含むLGM)を当てはめる:
growth() を用いた基本例model.growth <- '
i =~ 1*t1 + 1*t2 + 1*t3 + 1*t4
s =~ 0*t1 + 1*t2 + 2*t3 + 3*t4
'
fit.growth <- growth(model.growth, data = Demo.growth)
summary(fit.growth)
lavaan 0.6-21 ended normally after 29 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 9
Number of observations 400
Model Test User Model:
Test statistic 8.069
Degrees of freedom 5
P-value (Chi-square) 0.152
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Latent Variables:
Estimate Std.Err z-value P(>|z|)
i =~
t1 1.000
t2 1.000
t3 1.000
t4 1.000
s =~
t1 0.000
t2 1.000
t3 2.000
t4 3.000
Covariances:
Estimate Std.Err z-value P(>|z|)
i ~~
s 0.618 0.071 8.686 0.000
Intercepts:
Estimate Std.Err z-value P(>|z|)
i 0.615 0.077 8.007 0.000
s 1.006 0.042 24.076 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
.t1 0.595 0.086 6.944 0.000
.t2 0.676 0.061 11.061 0.000
.t3 0.635 0.072 8.761 0.000
.t4 0.508 0.124 4.090 0.000
i 1.932 0.173 11.194 0.000
s 0.587 0.052 11.336 0.000
cfa() を用いた明示的な指定cfa()を用いた設定も可能である。一見冗長にはなるが、本書の中で提供されているコードは、表現の統一およびモデル上の仮定の明示の観点から、基本的にLGMにおいてもcfa()を利用している。
model.growth2 <- '
i =~ 1*t1 + 1*t2 + 1*t3 + 1*t4
s =~ 0*t1 + 1*t2 + 2*t3 + 3*t4
t1 ~ 0*1
t2 ~ 0*1
t3 ~ 0*1
t4 ~ 0*1
i ~ mui*1
s ~ mus*1
i ~~ phii*i
s ~~ phis*s
i ~~ phiis*s
t1 ~~ psi1*t1
t2 ~~ psi2*t2
t3 ~~ psi3*t3
t4 ~~ psi4*t4
'
fit.growth2 <- cfa(model.growth2, data = Demo.growth)
summary(fit.growth2)
lavaan 0.6-21 ended normally after 29 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 9
Number of observations 400
Model Test User Model:
Test statistic 8.069
Degrees of freedom 5
P-value (Chi-square) 0.152
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Latent Variables:
Estimate Std.Err z-value P(>|z|)
i =~
t1 1.000
t2 1.000
t3 1.000
t4 1.000
s =~
t1 0.000
t2 1.000
t3 2.000
t4 3.000
Covariances:
Estimate Std.Err z-value P(>|z|)
i ~~
s (phis) 0.618 0.071 8.686 0.000
Intercepts:
Estimate Std.Err z-value P(>|z|)
.t1 0.000
.t2 0.000
.t3 0.000
.t4 0.000
i (mui) 0.615 0.077 8.007 0.000
s (mus) 1.006 0.042 24.076 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
i (phii) 1.932 0.173 11.194 0.000
s (phis) 0.587 0.052 11.336 0.000
.t1 (psi1) 0.595 0.086 6.944 0.000
.t2 (psi2) 0.676 0.061 11.061 0.000
.t3 (psi3) 0.635 0.072 8.761 0.000
.t4 (psi4) 0.508 0.124 4.090 0.000
時不変的な説明変数(x1,x2)の他、時変的な説明変数(c1,c2,c3,c4)を含む場合の条件付LGM。
model.cg <- '
i =~ 1*t1 + 1*t2 + 1*t3 + 1*t4
s =~ 0*t1 + 1*t2 + 2*t3 + 3*t4
i ~ x1 + x2
s ~ x1 + x2
t1 ~ c1
t2 ~ c2
t3 ~ c3
t4 ~ c4
t1 ~ 0*1
t2 ~ 0*1
t3 ~ 0*1
t4 ~ 0*1
i ~ inti*1
s ~ ints*1
i ~~ phii*i
s ~~ phis*s
i ~~ phiis*s
t1 ~~ psi1*t1
t2 ~~ psi2*t2
t3 ~~ psi3*t3
t4 ~~ psi4*t4
'
fit.conditionalgrowth <- cfa(model.cg, data = Demo.growth)
summary(fit.conditionalgrowth)
lavaan 0.6-21 ended normally after 31 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 17
Number of observations 400
Model Test User Model:
Test statistic 26.059
Degrees of freedom 21
P-value (Chi-square) 0.204
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Latent Variables:
Estimate Std.Err z-value P(>|z|)
i =~
t1 1.000
t2 1.000
t3 1.000
t4 1.000
s =~
t1 0.000
t2 1.000
t3 2.000
t4 3.000
Regressions:
Estimate Std.Err z-value P(>|z|)
i ~
x1 0.608 0.060 10.134 0.000
x2 0.604 0.064 9.412 0.000
s ~
x1 0.262 0.029 9.198 0.000
x2 0.522 0.031 17.083 0.000
t1 ~
c1 0.143 0.050 2.883 0.004
t2 ~
c2 0.289 0.046 6.295 0.000
t3 ~
c3 0.328 0.044 7.361 0.000
t4 ~
c4 0.330 0.058 5.655 0.000
Covariances:
Estimate Std.Err z-value P(>|z|)
.i ~~
.s (phis) 0.075 0.040 1.855 0.064
Intercepts:
Estimate Std.Err z-value P(>|z|)
.t1 0.000
.t2 0.000
.t3 0.000
.t4 0.000
.i (inti) 0.580 0.062 9.368 0.000
.s (ints) 0.958 0.029 32.552 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
.i (phii) 1.079 0.112 9.609 0.000
.s (phis) 0.224 0.027 8.429 0.000
.t1 (psi1) 0.580 0.080 7.230 0.000
.t2 (psi2) 0.596 0.054 10.969 0.000
.t3 (psi3) 0.481 0.055 8.745 0.000
.t4 (psi4) 0.535 0.098 5.466 0.000
別添の(欠測値NAを含む、人工データである)
semexample.csvデータを利用した、交差遅延パネルモデルによる分析例を示す。利用者は、まずこの
semexample.csvデータを
ダウンロードし、lavaan_tutorial_ja.Rmd
と同じフォルダに保存しておくことを前提とする。
SLEEP <- read.csv("C:/Users/Satoshi Usami/Desktop/semexample.csv", header = TRUE) #分析者のPCの環境に合わせて変更する必要あり。データの読み込みの方法は、本資料の最後(21節)を参照のこと。
head(SLEEP)
ID Sex GHQ1 GHQ2 GHQ3 GHQ4 GHQ5 GHQ6 sleep1 sleep2 sleep3 sleep4 sleep5
1 1 NA 0 NA NA NA NA NA NA NA NA NA NA
2 2 1 3 NA NA NA NA NA 9.0 NA NA NA NA
3 3 1 0 NA NA NA NA NA 8.5 NA NA NA NA
4 4 NA 0 NA NA NA NA NA 9.0 NA NA NA NA
5 5 NA 0 NA NA NA NA NA 8.0 NA NA NA NA
6 6 NA 0 NA NA NA NA NA 8.5 NA NA NA NA
sleep6 bed1 bed2 bed3 bed4 bed5 bed6
1 NA NA NA NA NA NA NA
2 NA 21.0 NA NA NA NA NA
3 NA 21.0 NA NA NA NA NA
4 NA 21.0 NA NA NA NA NA
5 NA 21.5 NA NA NA NA NA
6 NA 21.5 NA NA NA NA NA
CLPMT4 <- '
sleep1 ~ muy1*1
sleep2 ~ inty2*1
sleep3 ~ inty3*1
sleep4 ~ inty4*1
GHQ1 ~ mux1*1
GHQ2 ~ intx2*1
GHQ3 ~ intx3*1
GHQ4 ~ intx4*1
sleep2 ~ betay2*sleep1 + gammay2*GHQ1
GHQ2 ~ betax2*GHQ1 + gammax2*sleep1
sleep3 ~ betay3*sleep2 + gammay3*GHQ2
GHQ3 ~ betax3*GHQ2 + gammax3*sleep2
sleep4 ~ betay4*sleep3 + gammay4*GHQ3
GHQ4 ~ betax4*GHQ3 + gammax4*sleep3
GHQ2 ~~ Omegax2*GHQ2
sleep2 ~~ Omegay2*sleep2
GHQ2 ~~ Omegaxy2*sleep2
GHQ3 ~~ Omegax3*GHQ3
sleep3 ~~ Omegay3*sleep3
GHQ3 ~~ Omegaxy3*sleep3
GHQ4 ~~ Omegax4*GHQ4
sleep4 ~~ Omegay4*sleep4
GHQ4 ~~ Omegaxy4*sleep4
'
fit.clpm <- sem(CLPMT4, data = SLEEP, missing = "fiml") #fimlは完全情報最尤推定法による欠測処理
summary(fit.clpm, fit.measures = TRUE)
lavaan 0.6-21 ended normally after 76 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 31
Used Total
Number of observations 826 1067
Number of missing patterns 30
Model Test User Model:
Test statistic 111.985
Degrees of freedom 13
P-value (Chi-square) 0.000
Model Test Baseline Model:
Test statistic 850.195
Degrees of freedom 28
P-value 0.000
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.880
Tucker-Lewis Index (TLI) 0.741
Robust Comparative Fit Index (CFI) 0.829
Robust Tucker-Lewis Index (TLI) 0.631
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -7388.811
Loglikelihood unrestricted model (H1) -7332.819
Akaike (AIC) 14839.622
Bayesian (BIC) 14985.837
Sample-size adjusted Bayesian (SABIC) 14887.392
Root Mean Square Error of Approximation:
RMSEA 0.096
90 Percent confidence interval - lower 0.080
90 Percent confidence interval - upper 0.113
P-value H_0: RMSEA <= 0.050 0.000
P-value H_0: RMSEA >= 0.080 0.951
Robust RMSEA 0.174
90 Percent confidence interval - lower 0.143
90 Percent confidence interval - upper 0.206
P-value H_0: Robust RMSEA <= 0.050 0.000
P-value H_0: Robust RMSEA >= 0.080 1.000
Standardized Root Mean Square Residual:
SRMR 0.095
Parameter Estimates:
Standard errors Standard
Information Observed
Observed information based on Hessian
Regressions:
Estimate Std.Err z-value P(>|z|)
sleep2 ~
sleep1 (bty2) 0.573 0.045 12.742 0.000
GHQ1 (gmmy2) -0.036 0.020 -1.811 0.070
GHQ2 ~
GHQ1 (btx2) 0.569 0.045 12.576 0.000
sleep1 (gmmx2) -0.032 0.103 -0.310 0.757
sleep3 ~
sleep2 (bty3) 0.440 0.041 10.637 0.000
GHQ2 (gmmy3) -0.015 0.019 -0.787 0.431
GHQ3 ~
GHQ2 (btx3) 0.724 0.046 15.738 0.000
sleep2 (gmmx3) -0.172 0.117 -1.469 0.142
sleep4 ~
sleep3 (bty4) 0.393 0.043 9.145 0.000
GHQ3 (gmmy4) -0.025 0.015 -1.655 0.098
GHQ4 ~
GHQ3 (btx4) 0.467 0.041 11.399 0.000
sleep3 (gmmx4) -0.208 0.119 -1.752 0.080
Covariances:
Estimate Std.Err z-value P(>|z|)
.sleep2 ~~
.GHQ2 (Omg2) -0.150 0.102 -1.463 0.144
.sleep3 ~~
.GHQ3 (Omg3) -0.340 0.114 -2.981 0.003
.sleep4 ~~
.GHQ4 (Omg4) -0.098 0.112 -0.872 0.383
Intercepts:
Estimate Std.Err z-value P(>|z|)
sleep1 (muy1) 7.476 0.043 174.757 0.000
.sleep2 (inty2) 2.911 0.348 8.364 0.000
.sleep3 (inty3) 3.907 0.310 12.608 0.000
.sleep4 (inty4) 4.039 0.315 12.815 0.000
GHQ1 (mux1) 2.028 0.099 20.557 0.000
.GHQ2 (intx2) 1.596 0.796 2.005 0.045
.GHQ3 (intx3) 2.355 0.883 2.666 0.008
.GHQ4 (intx4) 2.846 0.871 3.267 0.001
Variances:
Estimate Std.Err z-value P(>|z|)
.GHQ2 (Omgx2) 4.768 0.326 14.618 0.000
.sleep2 (Omgy2) 0.959 0.062 15.344 0.000
.GHQ3 (Omgx3) 5.262 0.411 12.790 0.000
.sleep3 (Omgy3) 0.890 0.059 15.189 0.000
.GHQ4 (Omgx4) 6.250 0.436 14.349 0.000
.sleep4 (Omgy4) 0.822 0.057 14.405 0.000
sleep1 1.077 0.063 17.200 0.000
GHQ1 5.796 0.344 16.870 0.000
交差遅延パネルモデルからの様々な派生モデルの種類やその概説、lavaanコードについては以下も参照のこと.
SEMでは、例えば因子分析モデルを扱う場合のように、共分散(相関)構造のみに関心がもたれる場合も多いが、例えば潜在成長モデル(LGM; 本書の4章)のように平均構造に関心がもたれることも多い。より一般に、平均共分散構造分析(MACS)としてのSEMでは、変数の切片や平均も推定の対象となる。先述した3因子のCFAモデルでは,観測変数の切片を次のようにして加えることができる。
HS.mean <- '
visual =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed =~ x7 + x8 + x9
#以下、切片の設定
x1 ~ 1
x2 ~ 1
x3 ~ 1
x4 ~ 1
x5 ~ 1
x6 ~ 1
x7 ~ 1
x8 ~ 1
x9 ~ 1
'
meanstructure = TRUEモデル構文で切片を明示しないかわりにcfa()においてmeanstructure=TRUEを加える方法もある:
fit.mean <- cfa(HS.model,
data = HolzingerSwineford1939,
meanstructure = TRUE)
summary(fit.mean)
lavaan 0.6-21 ended normally after 35 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 30
Number of observations 301
Model Test User Model:
Test statistic 85.306
Degrees of freedom 24
P-value (Chi-square) 0.000
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
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
Intercepts:
Estimate Std.Err z-value P(>|z|)
.x1 4.936 0.067 73.473 0.000
.x2 6.088 0.068 89.855 0.000
.x3 2.250 0.065 34.579 0.000
.x4 3.061 0.067 45.694 0.000
.x5 4.341 0.074 58.452 0.000
.x6 2.186 0.063 34.667 0.000
.x7 4.186 0.063 66.766 0.000
.x8 5.527 0.058 94.854 0.000
.x9 5.374 0.058 92.546 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
データセット内の群変数(集団=群を表す変数)の名前を引数groupを通して与えることで、多群モデルの分析(本書の6章)ができる。デフォルトでは、各群に対して同じ(つまり、群間の制約がない)モデルが当てはめられる。次の例では、2つの学校(PasteurとGrant-White)のデータセットに対して先述の3因子のCFAモデルを当てはめている。
fit.group <- cfa(HS.model,
data = HolzingerSwineford1939,
group = "school")
summary(fit.group)
lavaan 0.6-21 ended normally after 57 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 60
Number of observations per group:
Pasteur 156
Grant-White 145
Model Test User Model:
Test statistic 115.851
Degrees of freedom 48
P-value (Chi-square) 0.000
Test statistic for each group:
Pasteur 64.309
Grant-White 51.542
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Group 1 [Pasteur]:
Latent Variables:
Estimate Std.Err z-value P(>|z|)
visual =~
x1 1.000
x2 0.394 0.122 3.220 0.001
x3 0.570 0.140 4.076 0.000
textual =~
x4 1.000
x5 1.183 0.102 11.613 0.000
x6 0.875 0.077 11.421 0.000
speed =~
x7 1.000
x8 1.125 0.277 4.057 0.000
x9 0.922 0.225 4.104 0.000
Covariances:
Estimate Std.Err z-value P(>|z|)
visual ~~
textual 0.479 0.106 4.531 0.000
speed 0.185 0.077 2.397 0.017
textual ~~
speed 0.182 0.069 2.628 0.009
Intercepts:
Estimate Std.Err z-value P(>|z|)
.x1 4.941 0.095 52.249 0.000
.x2 5.984 0.098 60.949 0.000
.x3 2.487 0.093 26.778 0.000
.x4 2.823 0.092 30.689 0.000
.x5 3.995 0.105 38.183 0.000
.x6 1.922 0.079 24.321 0.000
.x7 4.432 0.087 51.181 0.000
.x8 5.563 0.078 71.214 0.000
.x9 5.418 0.079 68.440 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
.x1 0.298 0.232 1.286 0.198
.x2 1.334 0.158 8.464 0.000
.x3 0.989 0.136 7.271 0.000
.x4 0.425 0.069 6.138 0.000
.x5 0.456 0.086 5.292 0.000
.x6 0.290 0.050 5.780 0.000
.x7 0.820 0.125 6.580 0.000
.x8 0.510 0.116 4.406 0.000
.x9 0.680 0.104 6.516 0.000
visual 1.097 0.276 3.967 0.000
textual 0.894 0.150 5.963 0.000
speed 0.350 0.126 2.778 0.005
Group 2 [Grant-White]:
Latent Variables:
Estimate Std.Err z-value P(>|z|)
visual =~
x1 1.000
x2 0.736 0.155 4.760 0.000
x3 0.925 0.166 5.583 0.000
textual =~
x4 1.000
x5 0.990 0.087 11.418 0.000
x6 0.963 0.085 11.377 0.000
speed =~
x7 1.000
x8 1.226 0.187 6.569 0.000
x9 1.058 0.165 6.429 0.000
Covariances:
Estimate Std.Err z-value P(>|z|)
visual ~~
textual 0.408 0.098 4.153 0.000
speed 0.276 0.076 3.639 0.000
textual ~~
speed 0.222 0.073 3.022 0.003
Intercepts:
Estimate Std.Err z-value P(>|z|)
.x1 4.930 0.095 51.696 0.000
.x2 6.200 0.092 67.416 0.000
.x3 1.996 0.086 23.195 0.000
.x4 3.317 0.093 35.625 0.000
.x5 4.712 0.096 48.986 0.000
.x6 2.469 0.094 26.277 0.000
.x7 3.921 0.086 45.819 0.000
.x8 5.488 0.087 63.174 0.000
.x9 5.327 0.085 62.571 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
.x1 0.715 0.126 5.676 0.000
.x2 0.899 0.123 7.339 0.000
.x3 0.557 0.103 5.409 0.000
.x4 0.315 0.065 4.870 0.000
.x5 0.419 0.072 5.812 0.000
.x6 0.406 0.069 5.880 0.000
.x7 0.600 0.091 6.584 0.000
.x8 0.401 0.094 4.249 0.000
.x9 0.535 0.089 6.010 0.000
visual 0.604 0.160 3.762 0.000
textual 0.942 0.152 6.177 0.000
speed 0.461 0.118 3.910 0.000
lavaan では、外生的なカテゴリカル変数と内生的なカテゴリカル変数を区別して考える必要がある。
K-1
個のダミー変数に置き換える生データがなくても、標本共分散行列とサンプルサイズの情報が入手できればモデルの当てはめが可能である(例えば本書の2.4.4節)。
#共分散行列の読み込み
wheaton.cov <- lav_matrix_lower2full(
c(11.834,
6.947, 9.364,
6.819, 5.091, 12.532,
4.783, 5.028, 7.495, 9.986,
-3.839, -3.889, -3.841, -3.625, 9.610,
-21.899, -18.831, -21.748, -18.775, 35.522, 450.288)
)
#変数名の付与
colnames(wheaton.cov) <- rownames(wheaton.cov) <-
c("anomia67", "powerless67", "anomia71", "powerless71", "education", "sei")
#分析モデルの設定
wheaton.model <- '
ses =~ education + sei
alien67 =~ anomia67 + powerless67
alien71 =~ anomia71 + powerless71
alien71 ~ alien67 + ses
alien67 ~ ses
anomia67 ~~ anomia71
powerless67 ~~ powerless71
'
fit.cov <- sem(wheaton.model,
sample.cov = wheaton.cov,
sample.nobs = 932) #サンプルサイズを明示する。
summary(fit.cov, standardized = TRUE)
lavaan 0.6-21 ended normally after 84 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 17
Number of observations 932
Model Test User Model:
Test statistic 4.735
Degrees of freedom 4
P-value (Chi-square) 0.316
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
ses =~
education 1.000 2.607 0.842
sei 5.219 0.422 12.364 0.000 13.609 0.642
alien67 =~
anomia67 1.000 2.663 0.774
powerless67 0.979 0.062 15.895 0.000 2.606 0.852
alien71 =~
anomia71 1.000 2.850 0.805
powerless71 0.922 0.059 15.498 0.000 2.628 0.832
Regressions:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
alien71 ~
alien67 0.607 0.051 11.898 0.000 0.567 0.567
ses -0.227 0.052 -4.334 0.000 -0.207 -0.207
alien67 ~
ses -0.575 0.056 -10.195 0.000 -0.563 -0.563
Covariances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.anomia67 ~~
.anomia71 1.623 0.314 5.176 0.000 1.623 0.356
.powerless67 ~~
.powerless71 0.339 0.261 1.298 0.194 0.339 0.121
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.education 2.801 0.507 5.525 0.000 2.801 0.292
.sei 264.597 18.126 14.597 0.000 264.597 0.588
.anomia67 4.731 0.453 10.441 0.000 4.731 0.400
.powerless67 2.563 0.403 6.359 0.000 2.563 0.274
.anomia71 4.399 0.515 8.542 0.000 4.399 0.351
.powerless71 3.070 0.434 7.070 0.000 3.070 0.308
ses 6.798 0.649 10.475 0.000 1.000 1.000
.alien67 4.841 0.467 10.359 0.000 0.683 0.683
.alien71 4.083 0.404 10.104 0.000 0.503 0.503
連続変数のみを扱う場合、デフォルトは ML
である。その他、次の推定法が利用可能である(例えば、本書の2.4節や9.2節)。
# estimator = "ML"
# estimator = "GLS"
# estimator = "WLS"
# estimator = "DWLS"
# estimator = "ULS"
# estimator = "DLS"
# estimator = "PML"
ロバスト版としては、例えば次がある(標準誤差、カイ二乗検定統計量、およびそれに基づく適合度評価がロバストに補正され、通常、点推定値自体はMLと同じ)。
# estimator = "MLM"
# estimator = "MLMVS"
# estimator = "MLMV"
# estimator = "MLR"
lavaan
で最尤推定(estimator = "ML")またはそのロバスト版を用いるとき、デフォルトでは
“normal” アプローチが採用される。これは、標本共分散行列を
\(N\)
で割った形(いわゆる biased sample covariance
matrix)に基づいて内部的に計算し、カイ二乗統計量(尤度比検定統計量:-2logλ)も最小化された目的関数値に
\(N\)
を掛けて求めるという方法である。
これに対して、慣習的な理由で、不偏共分散行列、すなわち各要素を
\(N-1\)
で割った共分散行列に基づいて計算し、さらにカイ二乗統計量の計算でも
\(N-1\)を掛けることも多い(本書の2.4.1節;(N−1)fML)。その場合には、likelihood = "wishart"
を指定する。たとえば、先述のCFAモデルでは次のように書く。
fit <- cfa(
HS.model,
data = HolzingerSwineford1939,
likelihood = "wishart"
)
fit
lavaan 0.6-21 ended normally after 35 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 21
Number of observations 301
Model Test User Model:
Test statistic 85.022
Degrees of freedom 24
P-value (Chi-square) 0.000
likelihood = "wishart"
を指定すると、主として検定統計量のスケーリングが変わる。そのため、同じモデル・同じデータでも、lavaan
のデフォルト設定(normal)で得られるカイ二乗値と、Wishart
指定で得られるカイ二乗値はわずかに異なりうる(サンプルサイズが大きくなればこの違いは実用上無視してよい)。本書では便宜的に、3章以降、lavaan
のデフォルト設定(normal)の下でのカイ二乗値を示している。
FIML(例えば、本書の11.4節、11.6節)を使った欠測データの処理。
fit.missing <- sem(model.sem,
data = PoliticalDemocracy,
missing = "fiml")
summary(fit.missing)
blavaan (Merkle & Rosseel, 2018)および、組み込みデータセットPoliticalDemocracyと先述のmodel.sem(共通因子間の回帰モデル)を用いたベイズ推定の例を示す。blavaanは、lavaan に近い構文を用いながらSEMのベイズ推定が実装できるR パッケージであり、sem() に対応する bsem()や、他にもbcfa() などを使ってモデルを推定できる。推定ではJAGSやStanを利用したMCMCが用いられる。blavaanの詳細については上記サイトを参照されたい。
install.packages("blavaan")
パッケージ 'blavaan' は無事に展開され、MD5 サムもチェックされました
ダウンロードされたパッケージは、以下にあります
C:\Users\Satoshi Usami\AppData\Local\Temp\RtmpIxTDFF\downloaded_packages
library(blavaan)
#事前分布の設定。dpriors()では母数のクラスごとに同じ形式の分布を設定する
mydp <- dpriors(
nu = "normal(0,100)",
alpha = "normal(0,100)",
lambda = "normal(0,10)",
beta = "normal(0,10)",
theta = "gamma(1,0.01)",
psi = "gamma(1,0.01)"
)
fit.bayes <- bcfa(model.sem, #共通因子間の回帰モデル(6節)
data = PoliticalDemocracy,
dp = mydp, #mydpで事前分布を設定する。
sample = 2000)
SAMPLING FOR MODEL 'stanmarg' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.000957 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 9.57 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2500 [ 0%] (Warmup)
Chain 1: Iteration: 250 / 2500 [ 10%] (Warmup)
Chain 1: Iteration: 500 / 2500 [ 20%] (Warmup)
Chain 1: Iteration: 501 / 2500 [ 20%] (Sampling)
Chain 1: Iteration: 750 / 2500 [ 30%] (Sampling)
Chain 1: Iteration: 1000 / 2500 [ 40%] (Sampling)
Chain 1: Iteration: 1250 / 2500 [ 50%] (Sampling)
Chain 1: Iteration: 1500 / 2500 [ 60%] (Sampling)
Chain 1: Iteration: 1750 / 2500 [ 70%] (Sampling)
Chain 1: Iteration: 2000 / 2500 [ 80%] (Sampling)
Chain 1: Iteration: 2250 / 2500 [ 90%] (Sampling)
Chain 1: Iteration: 2500 / 2500 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 3.578 seconds (Warm-up)
Chain 1: 11.07 seconds (Sampling)
Chain 1: 14.648 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'stanmarg' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0.00032 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 3.2 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2500 [ 0%] (Warmup)
Chain 2: Iteration: 250 / 2500 [ 10%] (Warmup)
Chain 2: Iteration: 500 / 2500 [ 20%] (Warmup)
Chain 2: Iteration: 501 / 2500 [ 20%] (Sampling)
Chain 2: Iteration: 750 / 2500 [ 30%] (Sampling)
Chain 2: Iteration: 1000 / 2500 [ 40%] (Sampling)
Chain 2: Iteration: 1250 / 2500 [ 50%] (Sampling)
Chain 2: Iteration: 1500 / 2500 [ 60%] (Sampling)
Chain 2: Iteration: 1750 / 2500 [ 70%] (Sampling)
Chain 2: Iteration: 2000 / 2500 [ 80%] (Sampling)
Chain 2: Iteration: 2250 / 2500 [ 90%] (Sampling)
Chain 2: Iteration: 2500 / 2500 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 3.104 seconds (Warm-up)
Chain 2: 11.072 seconds (Sampling)
Chain 2: 14.176 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'stanmarg' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.000441 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 4.41 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2500 [ 0%] (Warmup)
Chain 3: Iteration: 250 / 2500 [ 10%] (Warmup)
Chain 3: Iteration: 500 / 2500 [ 20%] (Warmup)
Chain 3: Iteration: 501 / 2500 [ 20%] (Sampling)
Chain 3: Iteration: 750 / 2500 [ 30%] (Sampling)
Chain 3: Iteration: 1000 / 2500 [ 40%] (Sampling)
Chain 3: Iteration: 1250 / 2500 [ 50%] (Sampling)
Chain 3: Iteration: 1500 / 2500 [ 60%] (Sampling)
Chain 3: Iteration: 1750 / 2500 [ 70%] (Sampling)
Chain 3: Iteration: 2000 / 2500 [ 80%] (Sampling)
Chain 3: Iteration: 2250 / 2500 [ 90%] (Sampling)
Chain 3: Iteration: 2500 / 2500 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 3.596 seconds (Warm-up)
Chain 3: 11.83 seconds (Sampling)
Chain 3: 15.426 seconds (Total)
Chain 3:
Computing post-estimation metrics (including lvs if requested)...
summary(fit.bayes, fit.measures = TRUE)
blavaan 0.5.10 ended normally after 2000 iterations
Estimator BAYES
Optimization method MCMC
Number of model parameters 31
Number of observations 75
Statistic MargLogLik PPP
Value NA 0.505
Parameter Estimates:
Latent Variables:
Estimate Post.SD pi.lower pi.upper Rhat Prior
ind60 =~
x1 1.000
x2 2.303 0.154 2.020 2.625 1.001 normal(0,10)
x3 1.869 0.169 1.557 2.224 1.000 normal(0,10)
dem60 =~
y1 1.000
y2 1.354 0.215 0.965 1.805 1.003 normal(0,10)
y3 1.114 0.172 0.804 1.483 1.000 normal(0,10)
y4 1.364 0.187 1.043 1.781 1.001 normal(0,10)
dem65 =~
y5 1.000
y6 1.248 0.195 0.913 1.678 1.001 normal(0,10)
y7 1.334 0.185 1.016 1.750 1.001 normal(0,10)
y8 1.338 0.190 1.007 1.762 1.001 normal(0,10)
Regressions:
Estimate Post.SD pi.lower pi.upper Rhat Prior
dem60 ~
ind60 1.434 0.395 0.680 2.222 1.000 normal(0,10)
dem65 ~
ind60 0.511 0.229 0.052 0.964 1.000 normal(0,10)
dem60 0.882 0.112 0.682 1.119 1.000 normal(0,10)
Covariances:
Estimate Post.SD pi.lower pi.upper Rhat Prior
.y1 ~~
.y5 0.649 0.385 -0.034 1.489 1.001 lkj_corr(1)
.y2 ~~
.y4 1.236 0.683 -0.017 2.684 1.000 beta(1,1)
.y6 2.076 0.735 0.777 3.626 1.001 beta(1,1)
.y3 ~~
.y7 0.763 0.626 -0.379 2.129 1.000 lkj_corr(1)
.y4 ~~
.y8 0.240 0.467 -0.626 1.222 1.000 beta(1,1)
.y6 ~~
.y8 1.341 0.555 0.340 2.542 1.000 beta(1,1)
Variances:
Estimate Post.SD pi.lower pi.upper Rhat Prior
.x1 0.096 0.021 0.058 0.139 1.000 gamma(1,0.01)
.x2 0.057 0.059 0.003 0.213 1.001 gamma(1,0.01)
.x3 0.510 0.096 0.349 0.724 1.000 gamma(1,0.01)
.y1 2.103 0.509 1.240 3.230 1.001 gamma(1,0.01)
.y2 7.475 1.358 5.136 10.516 1.000 gamma(1,0.01)
.y3 5.355 1.024 3.676 7.662 1.000 gamma(1,0.01)
.y4 3.158 0.808 1.693 4.892 1.001 gamma(1,0.01)
.y5 2.508 0.521 1.622 3.700 1.000 gamma(1,0.01)
.y6 5.108 0.934 3.535 7.237 1.000 gamma(1,0.01)
.y7 3.647 0.766 2.359 5.393 1.000 gamma(1,0.01)
.y8 3.331 0.739 2.081 4.948 1.000 gamma(1,0.01)
ind60 0.420 0.086 0.279 0.613 1.000 gamma(1,0.01)
.dem60 3.604 0.938 2.010 5.686 1.001 gamma(1,0.01)
.dem65 0.055 0.094 0.003 0.344 1.000 gamma(1,0.01)
連続データに対する、(random interceptsを導入した、切片[平均部分]に関する集団レベルの変動[集団差]を考慮した)2レベルのSEMが実装できる(本書の8.3節)。以下は、Demo.twolevelのデータを用いた例。
Demo.twolevelデータの構造の確認。変数clusterによって集団が区別されている。
head(Demo.twolevel)
y1 y2 y3 y4 y5 y6 x1
1 0.2293216 1.3555232 -0.6911702 0.8028079 -0.3011085 -1.7260671 1.1739003
2 0.3085801 -1.8624397 -2.4179783 0.7659289 1.6750617 1.1764210 -1.0039958
3 0.2004934 -1.3400514 0.4376087 1.1974194 1.1951594 1.4988962 -0.4402545
4 1.0447982 -0.9624490 -0.4464898 -0.2027252 -0.4590574 1.1734061 -0.6253657
5 0.6881792 -0.4565633 -0.6422296 0.9900408 1.7662535 0.7944601 -0.8450025
6 -2.0687644 -0.5997856 0.3148418 0.6764432 -0.6519928 1.8405605 -0.7831784
x2 x3 w1 w2 cluster
1 -0.62315173 0.6470414 -0.2479975 -0.4989800 1
2 -0.56689380 0.0201264 -0.2479975 -0.4989800 1
3 -2.13432572 -0.4591246 -0.2479975 -0.4989800 1
4 -0.33688869 1.2852093 -0.2479975 -0.4989800 1
5 -0.04229954 1.5598970 -0.2479975 -0.4989800 1
6 -0.22441996 -0.3814231 -2.3219338 -0.6910567 2
model.ml <- '
#集団内(個人)レベルの設定
level: 1
fw =~ y1 + y2 + y3
fw ~ x1 + x2 + x3
#集団間(集団)レベルの設定
level: 2
fb =~ y1 + y2 + y3
fb ~ w1 + w2
'
fit.ml <- sem(model = model.ml,
data = Demo.twolevel,
cluster = "cluster")
summary(fit.ml)
lavaan 0.6-21 ended normally after 36 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 20
Number of observations 2500
Number of clusters [cluster] 200
Model Test User Model:
Test statistic 8.092
Degrees of freedom 10
P-value (Chi-square) 0.620
Parameter Estimates:
Standard errors Standard
Information Observed
Observed information based on Hessian
Level 1 [within]:
Latent Variables:
Estimate Std.Err z-value P(>|z|)
fw =~
y1 1.000
y2 0.774 0.034 22.671 0.000
y3 0.734 0.033 22.355 0.000
Regressions:
Estimate Std.Err z-value P(>|z|)
fw ~
x1 0.510 0.023 22.037 0.000
x2 0.407 0.022 18.273 0.000
x3 0.205 0.021 9.740 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
.y1 0.986 0.046 21.591 0.000
.y2 1.066 0.039 27.271 0.000
.y3 1.011 0.037 27.662 0.000
.fw 0.546 0.040 13.539 0.000
Level 2 [cluster]:
Latent Variables:
Estimate Std.Err z-value P(>|z|)
fb =~
y1 1.000
y2 0.717 0.052 13.824 0.000
y3 0.587 0.048 12.329 0.000
Regressions:
Estimate Std.Err z-value P(>|z|)
fb ~
w1 0.165 0.079 2.093 0.036
w2 0.131 0.076 1.715 0.086
Intercepts:
Estimate Std.Err z-value P(>|z|)
.y1 0.024 0.075 0.327 0.743
.y2 -0.016 0.060 -0.269 0.788
.y3 -0.042 0.054 -0.777 0.437
Variances:
Estimate Std.Err z-value P(>|z|)
.y1 0.058 0.047 1.213 0.225
.y2 0.120 0.031 3.825 0.000
.y3 0.149 0.028 5.319 0.000
.fb 0.899 0.118 7.592 0.000
従属変数Y、独立変数X、媒介変数Mの3変数がある状況を考える。数値例として、これら3変数を含むデータセットを作り、XがYに直接影響を与え(=直接効果)、XがMを媒介してYに影響を与えること(=間接効果)を仮定した媒介モデルを当てはめる:
set.seed(1234)
X <- rnorm(100)
M <- 0.5 * X + rnorm(100)
Y <- 0.7 * M + rnorm(100)
Data.med <- data.frame(X = X, Y = Y, M = M)
model.med <- '
Y ~ c*X
M ~ a*X
Y ~ b*M
ab := a*b
total := c + (a*b)
'
fit.med <- sem(model.med, data = Data.med)
summary(fit.med)
lavaan 0.6-21 ended normally after 1 iteration
Estimator ML
Optimization method NLMINB
Number of model parameters 5
Number of observations 100
Model Test User Model:
Test statistic 0.000
Degrees of freedom 0
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Regressions:
Estimate Std.Err z-value P(>|z|)
Y ~
X (c) 0.036 0.104 0.348 0.728
M ~
X (a) 0.474 0.103 4.613 0.000
Y ~
M (b) 0.788 0.092 8.539 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
.Y 0.898 0.127 7.071 0.000
.M 1.054 0.149 7.071 0.000
Defined Parameters:
Estimate Std.Err z-value P(>|z|)
ab 0.374 0.092 4.059 0.000
total 0.410 0.125 3.287 0.001
ブートストラップ標準誤差を使うなら次のようにする。
fit.med.boot <- sem(model.med,
data = Data.med,
se = "bootstrap",
bootstrap = 1000)
summary(fit.med.boot)
分析モデルにおいて(0 などに)固定された全ての母数について計算される、自由度1 のカイ二乗検定に基づく指標。もし固定せずに推定した場合に検定統計量がどの程度低下するか を近似的に求めている。
modindices(fit.cfa, sort = TRUE, maximum.number = 5)
lhs op rhs mi epc sepc.lv sepc.all sepc.nox
30 visual =~ x9 36.411 0.577 0.519 0.515 0.515
76 x7 ~~ x8 34.145 0.536 0.536 0.859 0.859
28 visual =~ x7 18.631 -0.422 -0.380 -0.349 -0.349
78 x8 ~~ x9 14.946 -0.423 -0.423 -0.805 -0.805
33 textual =~ x3 9.151 -0.272 -0.269 -0.238 -0.238
因子負荷(“=~”に相当する母数)だけに絞って確認するには次のようにする。
mi <- modindices(fit.cfa)
mi[mi$op == "=~", ]
lhs op rhs mi epc sepc.lv sepc.all sepc.nox
25 visual =~ x4 1.211 0.077 0.069 0.059 0.059
26 visual =~ x5 7.441 -0.210 -0.189 -0.147 -0.147
27 visual =~ x6 2.843 0.111 0.100 0.092 0.092
28 visual =~ x7 18.631 -0.422 -0.380 -0.349 -0.349
29 visual =~ x8 4.295 -0.210 -0.189 -0.187 -0.187
30 visual =~ x9 36.411 0.577 0.519 0.515 0.515
31 textual =~ x1 8.903 0.350 0.347 0.297 0.297
32 textual =~ x2 0.017 -0.011 -0.011 -0.010 -0.010
33 textual =~ x3 9.151 -0.272 -0.269 -0.238 -0.238
34 textual =~ x7 0.098 -0.021 -0.021 -0.019 -0.019
35 textual =~ x8 3.359 -0.121 -0.120 -0.118 -0.118
36 textual =~ x9 4.796 0.138 0.137 0.136 0.136
37 speed =~ x1 0.014 0.024 0.015 0.013 0.013
38 speed =~ x2 1.580 -0.198 -0.123 -0.105 -0.105
39 speed =~ x3 0.716 0.136 0.084 0.075 0.075
40 speed =~ x4 0.003 -0.005 -0.003 -0.003 -0.003
41 speed =~ x5 0.201 -0.044 -0.027 -0.021 -0.021
42 speed =~ x6 0.273 0.044 0.027 0.025 0.025
等値制約を外した場合の検討には lavTestScore()
が有用である。
coef()母数の点推定値の抽出。
coef(fit.sem)
ind60=~x2 ind60=~x3 dem60=~y2 dem60=~y3 dem60=~y4 dem65=~y6
2.180 1.819 1.257 1.058 1.265 1.186
dem65=~y7 dem65=~y8 dem60~ind60 dem65~ind60 dem65~dem60 y1~~y5
1.280 1.266 1.483 0.572 0.837 0.624
y2~~y4 y2~~y6 y3~~y7 y4~~y8 y6~~y8 x1~~x1
1.313 2.153 0.795 0.348 1.356 0.082
x2~~x2 x3~~x3 y1~~y1 y2~~y2 y3~~y3 y4~~y4
0.120 0.467 1.891 7.373 5.067 3.148
y5~~y5 y6~~y6 y7~~y7 y8~~y8 ind60~~ind60 dem60~~dem60
2.351 4.954 3.431 3.254 0.448 3.956
dem65~~dem65
0.172
parameterEstimates()母数の点推定値やSE, 検定統計量等の抽出。
parameterEstimates(fit.sem)
lhs op rhs est se z pvalue ci.lower ci.upper
1 ind60 =~ x1 1.000 0.000 NA NA 1.000 1.000
2 ind60 =~ x2 2.180 0.139 15.742 0.000 1.909 2.452
3 ind60 =~ x3 1.819 0.152 11.967 0.000 1.521 2.116
4 dem60 =~ y1 1.000 0.000 NA NA 1.000 1.000
5 dem60 =~ y2 1.257 0.182 6.889 0.000 0.899 1.614
6 dem60 =~ y3 1.058 0.151 6.987 0.000 0.761 1.354
7 dem60 =~ y4 1.265 0.145 8.722 0.000 0.981 1.549
8 dem65 =~ y5 1.000 0.000 NA NA 1.000 1.000
9 dem65 =~ y6 1.186 0.169 7.024 0.000 0.855 1.517
10 dem65 =~ y7 1.280 0.160 8.002 0.000 0.966 1.593
11 dem65 =~ y8 1.266 0.158 8.007 0.000 0.956 1.576
12 dem60 ~ ind60 1.483 0.399 3.715 0.000 0.701 2.265
13 dem65 ~ ind60 0.572 0.221 2.586 0.010 0.139 1.006
14 dem65 ~ dem60 0.837 0.098 8.514 0.000 0.645 1.030
15 y1 ~~ y5 0.624 0.358 1.741 0.082 -0.079 1.326
16 y2 ~~ y4 1.313 0.702 1.871 0.061 -0.063 2.689
17 y2 ~~ y6 2.153 0.734 2.934 0.003 0.715 3.591
18 y3 ~~ y7 0.795 0.608 1.308 0.191 -0.396 1.986
19 y4 ~~ y8 0.348 0.442 0.787 0.431 -0.519 1.215
20 y6 ~~ y8 1.356 0.568 2.386 0.017 0.242 2.470
21 x1 ~~ x1 0.082 0.019 4.184 0.000 0.043 0.120
22 x2 ~~ x2 0.120 0.070 1.718 0.086 -0.017 0.256
23 x3 ~~ x3 0.467 0.090 5.177 0.000 0.290 0.643
24 y1 ~~ y1 1.891 0.444 4.256 0.000 1.020 2.762
25 y2 ~~ y2 7.373 1.374 5.366 0.000 4.680 10.066
26 y3 ~~ y3 5.067 0.952 5.324 0.000 3.202 6.933
27 y4 ~~ y4 3.148 0.739 4.261 0.000 1.700 4.596
28 y5 ~~ y5 2.351 0.480 4.895 0.000 1.410 3.292
29 y6 ~~ y6 4.954 0.914 5.419 0.000 3.162 6.746
30 y7 ~~ y7 3.431 0.713 4.814 0.000 2.034 4.829
31 y8 ~~ y8 3.254 0.695 4.685 0.000 1.893 4.615
32 ind60 ~~ ind60 0.448 0.087 5.173 0.000 0.279 0.618
33 dem60 ~~ dem60 3.956 0.921 4.295 0.000 2.151 5.762
34 dem65 ~~ dem65 0.172 0.215 0.803 0.422 -0.249 0.593
standardizedSolution()標準化解の抽出。
standardizedSolution(fit.sem)
lhs op rhs est.std se z pvalue ci.lower ci.upper
1 ind60 =~ x1 0.920 0.023 40.005 0.000 0.875 0.965
2 ind60 =~ x2 0.973 0.016 59.110 0.000 0.941 1.005
3 ind60 =~ x3 0.872 0.031 28.081 0.000 0.811 0.933
4 dem60 =~ y1 0.850 0.042 20.200 0.000 0.768 0.933
5 dem60 =~ y2 0.717 0.064 11.157 0.000 0.591 0.843
6 dem60 =~ y3 0.722 0.063 11.439 0.000 0.599 0.846
7 dem60 =~ y4 0.846 0.043 19.462 0.000 0.761 0.931
8 dem65 =~ y5 0.808 0.048 16.847 0.000 0.714 0.902
9 dem65 =~ y6 0.746 0.058 12.924 0.000 0.633 0.859
10 dem65 =~ y7 0.824 0.045 18.392 0.000 0.736 0.911
11 dem65 =~ y8 0.828 0.045 18.477 0.000 0.740 0.916
12 dem60 ~ ind60 0.447 0.103 4.323 0.000 0.244 0.649
13 dem65 ~ ind60 0.182 0.070 2.587 0.010 0.044 0.320
14 dem65 ~ dem60 0.885 0.051 17.405 0.000 0.786 0.985
15 y1 ~~ y5 0.296 0.140 2.117 0.034 0.022 0.570
16 y2 ~~ y4 0.273 0.120 2.265 0.024 0.037 0.508
17 y2 ~~ y6 0.356 0.097 3.690 0.000 0.167 0.545
18 y3 ~~ y7 0.191 0.133 1.432 0.152 -0.070 0.452
19 y4 ~~ y8 0.109 0.131 0.832 0.405 -0.147 0.365
20 y6 ~~ y8 0.338 0.110 3.079 0.002 0.123 0.553
21 x1 ~~ x1 0.154 0.042 3.637 0.000 0.071 0.237
22 x2 ~~ x2 0.053 0.032 1.661 0.097 -0.010 0.116
23 x3 ~~ x3 0.239 0.054 4.419 0.000 0.133 0.346
24 y1 ~~ y1 0.277 0.072 3.865 0.000 0.136 0.417
25 y2 ~~ y2 0.486 0.092 5.269 0.000 0.305 0.666
26 y3 ~~ y3 0.478 0.091 5.242 0.000 0.299 0.657
27 y4 ~~ y4 0.285 0.073 3.875 0.000 0.141 0.429
28 y5 ~~ y5 0.347 0.078 4.478 0.000 0.195 0.499
29 y6 ~~ y6 0.443 0.086 5.149 0.000 0.275 0.612
30 y7 ~~ y7 0.322 0.074 4.359 0.000 0.177 0.466
31 y8 ~~ y8 0.315 0.074 4.242 0.000 0.169 0.460
32 ind60 ~~ ind60 1.000 0.000 NA NA 1.000 1.000
33 dem60 ~~ dem60 0.800 0.092 8.669 0.000 0.619 0.981
34 dem65 ~~ dem65 0.039 0.048 0.808 0.419 -0.056 0.134
fitted()推定されたモデルの共分散構造(および平均構造)の出力。
fitted(fit.cfa)
$cov
x1 x2 x3 x4 x5 x6 x7 x8 x9
x1 1.358
x2 0.448 1.382
x3 0.590 0.327 1.275
x4 0.408 0.226 0.298 1.351
x5 0.454 0.252 0.331 1.090 1.660
x6 0.378 0.209 0.276 0.907 1.010 1.196
x7 0.262 0.145 0.191 0.173 0.193 0.161 1.183
x8 0.309 0.171 0.226 0.205 0.228 0.190 0.453 1.022
x9 0.284 0.157 0.207 0.188 0.209 0.174 0.415 0.490 1.015
resid() / residuals()残差平均ベクトル・残差共分散行列・残差相関行列(本書の2.5.2節)
resid(fit.cfa, type = "standardized")
$type
[1] "standardized"
$cov
x1 x2 x3 x4 x5 x6 x7 x8 x9
x1 0.000
x2 -1.996 0.000
x3 -0.997 2.689 0.000
x4 2.679 -0.284 -1.899 0.000
x5 -0.359 -0.591 -4.157 1.545 0.000
x6 2.155 0.681 -0.711 -2.588 0.942 0.000
x7 -3.773 -3.654 -1.858 0.865 -0.842 -0.326 0.000
x8 -1.380 -1.119 -0.300 -2.021 -1.099 -0.641 4.823 0.000
x9 4.077 1.606 3.518 1.225 1.701 1.423 -2.325 -4.132 0.000
vcov()推定値の漸近共分散行列(本書の2.4.3節)
vcov(fit.sem)
i60=~2 i60=~3 d60=~2 d60=~3 d60=~4 d65=~6 d65=~7 d65=~8 d60~60
ind60=~x2 0.019
ind60=~x3 0.010 0.023
dem60=~y2 0.000 0.000 0.033
dem60=~y3 0.000 0.000 0.008 0.023
dem60=~y4 0.000 0.000 0.015 0.009 0.021
dem65=~y6 0.000 0.000 0.010 0.003 0.004 0.028
dem65=~y7 0.000 0.000 0.004 0.006 0.004 0.013 0.026
dem65=~y8 0.000 0.000 0.004 0.003 0.005 0.017 0.013 0.025
dem60~ind60 0.009 0.007 -0.012 -0.010 -0.012 -0.004 -0.005 -0.005 0.159
dem65~ind60 0.003 0.003 -0.001 -0.002 -0.001 -0.005 -0.006 -0.005 -0.011
dem65~dem60 0.000 0.000 0.004 0.004 0.004 -0.006 -0.007 -0.007 -0.003
y1~~y5 0.000 0.000 0.009 0.007 0.010 0.009 0.010 0.010 -0.007
y2~~y4 0.000 0.000 -0.025 -0.001 -0.017 -0.002 -0.002 -0.002 0.010
y2~~y6 0.000 0.000 -0.009 0.000 0.002 -0.007 0.000 0.003 0.000
y3~~y7 0.000 0.000 0.000 -0.014 0.000 0.000 -0.011 0.000 0.003
y4~~y8 0.000 0.000 -0.001 -0.001 -0.010 -0.001 -0.001 -0.010 0.005
y6~~y8 0.000 0.000 0.000 0.000 0.000 -0.013 -0.001 -0.011 0.000
x1~~x1 0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
x2~~x2 -0.004 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
x3~~x3 0.001 -0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.000
y1~~y1 0.000 0.000 0.015 0.013 0.017 0.006 0.007 0.007 -0.012
y2~~y2 0.000 0.000 -0.041 -0.001 -0.010 -0.008 -0.002 0.001 0.009
y3~~y3 0.000 0.000 -0.001 -0.023 0.000 -0.001 -0.007 0.000 0.006
y4~~y4 0.000 0.000 -0.015 -0.002 -0.028 -0.003 -0.002 -0.008 0.014
y5~~y5 0.000 0.000 0.005 0.005 0.006 0.012 0.013 0.014 -0.005
y6~~y6 0.000 0.000 -0.004 0.000 0.001 -0.020 -0.001 -0.006 0.000
y7~~y7 0.000 0.000 -0.001 -0.007 0.000 -0.001 -0.019 0.000 0.002
y8~~y8 0.000 0.000 -0.001 0.000 -0.006 -0.009 -0.001 -0.021 0.003
ind60~~ind60 -0.006 -0.004 0.000 0.000 0.000 0.000 0.000 0.000 -0.003
dem60~~dem60 0.001 0.000 -0.063 -0.053 -0.064 -0.021 -0.023 -0.023 0.039
dem65~~dem65 0.000 0.000 0.000 -0.001 0.000 -0.002 -0.003 -0.002 0.001
dm65~n60 dm65~d60 y1~~y5 y2~~y4 y2~~y6 y3~~y7 y4~~y8 y6~~y8 x1~~x1
ind60=~x2
ind60=~x3
dem60=~y2
dem60=~y3
dem60=~y4
dem65=~y6
dem65=~y7
dem65=~y8
dem60~ind60
dem65~ind60 0.049
dem65~dem60 -0.006 0.010
y1~~y5 0.000 -0.003 0.128
y2~~y4 -0.013 0.005 -0.029 0.493
y2~~y6 0.000 0.000 0.000 0.020 0.538
y3~~y7 0.006 -0.003 -0.012 -0.028 0.000 0.369
y4~~y8 0.003 -0.002 -0.017 -0.010 -0.064 -0.016 0.196
y6~~y8 0.002 0.003 -0.018 -0.034 0.009 -0.020 0.000 0.323
x1~~x1 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
x2~~x2 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -0.001
x3~~x3 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
y1~~y1 -0.008 0.007 0.086 -0.025 0.000 -0.023 -0.025 -0.003 0.000
y2~~y2 -0.013 0.005 -0.027 0.486 0.504 -0.025 -0.060 -0.034 0.000
y3~~y3 -0.002 0.001 -0.018 -0.019 0.000 0.230 -0.018 -0.004 0.000
y4~~y4 -0.011 0.003 -0.044 0.273 -0.037 -0.041 0.113 -0.010 0.000
y5~~y5 -0.002 -0.005 0.092 -0.009 0.000 -0.018 -0.018 -0.008 0.000
y6~~y6 0.002 0.003 -0.015 -0.028 0.321 -0.017 -0.040 0.281 0.000
y7~~y7 0.005 0.002 -0.024 -0.001 0.000 0.166 -0.023 -0.005 0.000
y8~~y8 0.004 0.002 -0.030 -0.015 -0.035 -0.031 0.105 0.227 0.000
ind60~~ind60 -0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
dem60~~dem60 0.032 -0.030 -0.042 0.013 0.000 0.017 0.018 0.004 0.000
dem65~~dem65 0.008 -0.004 0.012 -0.028 0.000 0.017 0.008 -0.026 0.000
x2~~x2 x3~~x3 y1~~y1 y2~~y2 y3~~y3 y4~~y4 y5~~y5 y6~~y6 y7~~y7
ind60=~x2
ind60=~x3
dem60=~y2
dem60=~y3
dem60=~y4
dem65=~y6
dem65=~y7
dem65=~y8
dem60~ind60
dem65~ind60
dem65~dem60
y1~~y5
y2~~y4
y2~~y6
y3~~y7
y4~~y8
y6~~y8
x1~~x1
x2~~x2 0.005
x3~~x3 -0.002 0.008
y1~~y1 0.000 0.000 0.198
y2~~y2 0.000 0.000 -0.021 1.888
y3~~y3 0.000 0.000 -0.021 -0.016 0.906
y4~~y4 0.000 0.000 -0.046 0.131 -0.034 0.546
y5~~y5 0.000 0.000 0.044 -0.007 -0.012 -0.024 0.231
y6~~y6 0.000 0.000 -0.001 0.140 -0.002 -0.025 -0.005 0.836
y7~~y7 0.000 0.000 -0.015 0.003 0.067 -0.019 -0.016 -0.003 0.508
y8~~y8 0.000 0.000 -0.019 -0.034 -0.015 0.038 -0.019 0.102 -0.020
ind60~~ind60 0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
dem60~~dem60 -0.001 0.000 -0.086 0.010 0.013 0.027 -0.027 0.002 0.012
dem65~~dem65 0.000 0.000 -0.010 -0.029 -0.004 -0.022 -0.010 -0.025 -0.019
y8~~y8 i60~~6 d60~~6 d65~~6
ind60=~x2
ind60=~x3
dem60=~y2
dem60=~y3
dem60=~y4
dem65=~y6
dem65=~y7
dem65=~y8
dem60~ind60
dem65~ind60
dem65~dem60
y1~~y5
y2~~y4
y2~~y6
y3~~y7
y4~~y8
y6~~y8
x1~~x1
x2~~x2
x3~~x3
y1~~y1
y2~~y2
y3~~y3
y4~~y4
y5~~y5
y6~~y6
y7~~y7
y8~~y8 0.482
ind60~~ind60 0.000 0.008
dem60~~dem60 0.015 0.000 0.849
dem65~~dem65 -0.023 0.000 0.013 0.046
fitMeasures()適合度指標の出力(本書の2.5.1節)
fitMeasures(fit.cfa, c("cfi", "rmsea", "srmr"))
cfi rmsea srmr
0.931 0.092 0.065
lavInspect() / inspect()当てはめたモデルオブジェクトの中身を調べたり、内部情報を取り出したりするための関数。例えば、設定した分析モデルの係数行列や分散共分散行列、母数の設定の具体を確認したいときに利用できる。
lavInspect(fit.cfa)
$lambda
visual textul speed
x1 0 0 0
x2 1 0 0
x3 2 0 0
x4 0 0 0
x5 0 3 0
x6 0 4 0
x7 0 0 0
x8 0 0 5
x9 0 0 6
$theta
x1 x2 x3 x4 x5 x6 x7 x8 x9
x1 7
x2 0 8
x3 0 0 9
x4 0 0 0 10
x5 0 0 0 0 11
x6 0 0 0 0 0 12
x7 0 0 0 0 0 0 13
x8 0 0 0 0 0 0 0 14
x9 0 0 0 0 0 0 0 0 15
$psi
visual textul speed
visual 16
textual 19 17
speed 20 21 18
inspect(fit.cfa)
$lambda
visual textul speed
x1 0 0 0
x2 1 0 0
x3 2 0 0
x4 0 0 0
x5 0 3 0
x6 0 4 0
x7 0 0 0
x8 0 0 5
x9 0 0 6
$theta
x1 x2 x3 x4 x5 x6 x7 x8 x9
x1 7
x2 0 8
x3 0 0 9
x4 0 0 0 10
x5 0 0 0 0 11
x6 0 0 0 0 0 12
x7 0 0 0 0 0 0 13
x8 0 0 0 0 0 0 0 14
x9 0 0 0 0 0 0 0 0 15
$psi
visual textul speed
visual 16
textual 19 17
speed 20 21 18
初期値を確認したい場合は次のようにする。
lavInspect(fit.cfa, what = "start")
$lambda
visual textul speed
x1 1.000 0.000 0.000
x2 0.778 0.000 0.000
x3 1.107 0.000 0.000
x4 0.000 1.000 0.000
x5 0.000 1.133 0.000
x6 0.000 0.924 0.000
x7 0.000 0.000 1.000
x8 0.000 0.000 1.225
x9 0.000 0.000 0.854
$theta
x1 x2 x3 x4 x5 x6 x7 x8 x9
x1 0.679
x2 0.000 0.691
x3 0.000 0.000 0.637
x4 0.000 0.000 0.000 0.675
x5 0.000 0.000 0.000 0.000 0.830
x6 0.000 0.000 0.000 0.000 0.000 0.598
x7 0.000 0.000 0.000 0.000 0.000 0.000 0.592
x8 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.511
x9 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.508
$psi
visual textul speed
visual 0.05
textual 0.00 0.05
speed 0.00 0.00 0.05
公式 tutorial では、モデル内に含まれる多重指標(e.g.,心理尺度)に関する因子構造が事前に明確でない場合に、因子の回転を用いて構造を探索しながら(モデル全体を)推定する枠組みとして ESEM / EFA (探索的構造方程式モデリング、探索的因子分析)が紹介されている。具体例は公式 tutorial の ESEM / EFA セクションを参照されたい。
read.csv() によるデータの読み込み本資料のCLPMの例では、別添の semexample.csv
(人工データ)を利用する。まず利用者はこのファイルをダウンロードし、lavaan_tutorial_ja.Rmd
と同じフォルダに保存する。そうすれば、read.csv()
を用いて次のように読み込める。
mydata <- read.csv("semexample.csv", header = TRUE)
head(mydata)
ここで、"semexample.csv"
は読み込むCSVファイル名を表しており、header = TRUE は 1
行目を変数名として扱うことを意味する。head(mydata)
は、読み込んだデータの冒頭数行を確認するための関数である。
上のような書き方ができるのは、semexample.csv を
lavaan_tutorial_ja.Rmd
と同じフォルダに置いているためである。別の場所に保存した場合には、その保存場所に応じてパスを書き換える必要がある。
Windows の場合、たとえば semexample.csv
をデスクトップに保存したなら、次のように書ける。
mydata <- read.csv("C:/Users/YourName/Desktop/semexample.csv", header = TRUE)
head(mydata)
宇佐美のPCの場合、以下の通り:
mydata <- read.csv("C:/Users/Satoshi Usami/Desktop/semexample.csv", header = TRUE)
head(mydata)
ID Sex GHQ1 GHQ2 GHQ3 GHQ4 GHQ5 GHQ6 sleep1 sleep2 sleep3 sleep4 sleep5
1 1 NA 0 NA NA NA NA NA NA NA NA NA NA
2 2 1 3 NA NA NA NA NA 9.0 NA NA NA NA
3 3 1 0 NA NA NA NA NA 8.5 NA NA NA NA
4 4 NA 0 NA NA NA NA NA 9.0 NA NA NA NA
5 5 NA 0 NA NA NA NA NA 8.0 NA NA NA NA
6 6 NA 0 NA NA NA NA NA 8.5 NA NA NA NA
sleep6 bed1 bed2 bed3 bed4 bed5 bed6
1 NA NA NA NA NA NA NA
2 NA 21.0 NA NA NA NA NA
3 NA 21.0 NA NA NA NA NA
4 NA 21.0 NA NA NA NA NA
5 NA 21.5 NA NA NA NA NA
6 NA 21.5 NA NA NA NA NA
Mac の場合は、たとえばホームディレクトリの Documents フォルダ内に保存したとき、次のような形になる。
mydata <- read.csv("/Users/YourName/Documents/semexample.csv",
header = TRUE)
head(mydata)
semexample.csv
をそのまま利用する限り通常は特別な指定は不要だが、環境によって文字コードの指定が必要になる場合には、必要に応じて
fileEncoding を指定する。
mydata <- read.csv("semexample.csv",
header = TRUE,
fileEncoding = "UTF-8")
head(mydata)
Windows の環境では、CSV ファイルが CP932 や
Shift-JIS
系で保存されていることもあり、その場合には次のような指定が必要になることがある。
mydata <- read.csv("semexample.csv",
header = TRUE,
fileEncoding = "CP932")
head(mydata)
RStudio では、Project を作成して作業する方法が便利である。Project は作業ディレクトリと結び付けられており、分析ごとに 1 つの Project を作成して、その中にデータ・スクリプト・出力ファイルをまとめて管理する方法が勧められている。Project の基本的な説明や作成方法については、Posit 公式ガイドの Projects in RStudio IDE および Get Started を参照のこと。
たとえば、次のようなフォルダ構成を考える。
project_folder/
lavaan_tutorial_ja.Rmd
semexample.csv
このとき、lavaan_tutorial_ja.Rmd と
semexample.csv が同じ Project
内の同じフォルダにあれば、次のように書ける。
mydata <- read.csv("semexample.csv", header = TRUE)
head(mydata)
この方法の利点は、他の PC にファイル一式を移しても、フォルダ構成が同じであればそのまま動きやすいことである。授業資料や共同研究でファイルを共有する場合にも有用である。
Bollen, K.A. (1989). Structural Equations with Latent Variables. New York: Wiley.
Merkle, E.C., & Rosseel, Y. (2018). blavaan: Bayesian structural equation models via parameter expansion. Journal of Statistical Software, 85, 1–30. https://www.jstatsoft.org/article/view/v085i04.
R Core Team (2024). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/
Rosseel, Y. (2012). lavaan: An R Package for Structural Equation Modeling. Journal of Statistical Software, 48(2), 1–36.
宇佐美慧 (2026). 縦断データ分析:構造方程式モデリングによるアプローチ 東京大学出版会