データの準備
data("mtcars")
write.csv(mtcars, "mtcars.csv")
head(mtcars) str(mtcars)
今回は2つのグループを作るため、軽い車と重い車の2値データを作成し、軽い車をGroup1、重い車をGroup2として分析を行う。
# データの読み込み data(mtcars) # wtの中央値を計算 wt_median <- median(mtcars$wt) # wt_group変数を作成 mtcars$wt_group <- ifelse(mtcars$wt > wt_median, 1, 0) # 結果の確認 head(mtcars)
> head(mtcars)
mpg cyl disp hp drat wt qsec vs am gear carb wt_group
Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4 0
Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4 0
Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1 0
Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1 0
Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2 1
Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1 1
wt_groupがグループを示す変数である。
多母集団同時分析
library(lavaan)
model.mtcars <- '
mpg ~ cyl + hp + wt + am
'
fit.mtcars <- sem(model.mtcars,
data = mtcars,
group = "wt_group",
missing = "ML",
estimator = "MLR")
summary(fit.mtcars, standardized = TRUE)
Estimator ML
Optimization method NLMINB
Number of model parameters 12
Number of observations per group:
0 16
1 16
Number of missing patterns per group:
0 1
1 1
Model Test User Model:
Standard Scaled
Test Statistic 0.000 0.000
Degrees of freedom 0 0
Test statistic for each group:
0 0.000 0.000
1 0.000 0.000
Parameter Estimates:
Standard errors Sandwich
Information bread Observed
Observed information based on Hessian
Group 1 [0]:
Regressions:
Estimate Std.Err z-value P(>|z|)
mpg ~
cyl 0.119 0.708 0.169 0.866
hp -0.044 0.018 -2.397 0.017
wt -5.475 1.523 -3.594 0.000
am -0.359 2.122 -0.169 0.866
Intercepts:
Estimate Std.Err z-value P(>|z|)
.mpg 42.430 4.352 9.750 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
.mpg 5.429 1.678 3.235 0.001
Group 2 [1]:
Regressions:
Estimate Std.Err z-value P(>|z|)
mpg ~
cyl -0.207 0.700 -0.296 0.767
hp -0.020 0.014 -1.460 0.144
wt -2.045 0.793 -2.580 0.010
am 1.702 1.896 0.898 0.369
Intercepts:
Estimate Std.Err z-value P(>|z|)
.mpg 29.019 3.446 8.421 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
.mpg 2.690 0.732 3.674 0.000
グループ1の hp の Std.all の値は -0.446、グループ2の hp の Std.all の値は -0.424で、これらの値の比をとると、-0.446 / -0.424 = 1.052 となる。 したがって、グループ2に比べて、グループ1の hp の mpg への影響は約 1.05倍 であると言える。 つまり、軽い車（グループ1）では、重い車（グループ2）に比べて、馬力が増加すると燃費が悪化する。
グループ間における特定の変数の比較
Group1(軽い車)ではhs(馬力)のP値が有意で、Group2(重い車)では有意ではなかった。 この結果でも軽い車と重い車でhsの影響が異なると言ってもよいのだが、直接比較するためにWald検定をする方法がある。
# 等値制約モデルを定義 (hp以外)
model.equal.hs <- '
# グループ間でパス係数を等値制約 (hp以外)
mpg ~ c(a1, a1)*cyl + c(b1, b2)*hp + c(c1, c1)*wt + c(d1, d1)*am
'
# 等値制約モデルを適用
fit.equal.hs <- sem(model.equal.hs,
data = mtcars,
group = "wt_group",
missing = "ML",
estimator = "MLR")
# スコア検定 (hpのパス係数を等値制約)
lavTestWald(fit.equal.hs, constraints = "b1 == b2")
$stat [1] 8.903259
$df [1] 1
$p.value [1] 0.002846622
$se [1] "robust.huber.white"
hs(馬力)は2つのグループ間で hp のパス係数に有意な差があることが示された。
cyl(シリンダー数)は2つのグループ間における回帰分析で有意な差がないことが示されている。 2つのグループ間で比較した時はどうだろうか。
# 等値制約モデルを定義 (cyl以外)
model.equal.cyl <- '
# グループ間でパス係数を等値制約 (cyl以外)
mpg ~ c(a1, a2)*cyl + c(b1, b1)*hp + c(c1, c1)*wt + c(d1, d1)*am
'
# 等値制約モデルを適用
fit.equal.cyl <- sem(model.equal.cyl,
data = mtcars,
group = "wt_group",
missing = "ML",
estimator = "MLR")
# スコア検定 (hpのパス係数を等値制約)
lavTestWald(fit.equal.cyl, constraints = "a1 == a2")
$stat [1] 8.052777
$df [1] 1
$p.value [1] 0.004543395
$se [1] "robust.huber.white"
回帰式では有意ではなかったが、2つのグループ間でシリンダー数のパス係数に有意な差があることが示された。Group1(軽い車)ではシリンダー数が燃費によい影響があり、Group2(重い車)では燃費に悪い影響があり、グループ間比較をすると有意に差があることがわかった。
