中介效果
一個獨變項，一個中介變項
#Baron & Kenny（1986） 的四步驟，檢驗被霸凌是體指數到憂鬱的中介
m1 <- lm(Y ~ X, data=dta)
m2 <- lm(M1 ~ X, data=dta)
m3 <- lm(M2 ~ M1, data=dta)
m4 <- lm(M2 ~ X+M1, data=dta)
m5 <- lm(Y ~ M2, data=dta)
m6 <- lm(Y ~ M1, data=dta)
m7 <- lm(Y ~ X+M1, data=dta)
m8 <- lm(Y ~ X+M2, data=dta)
m7 <- lm(Y ~ M1+M2, data=dta)
m9 <- lm(Y ~ X+M1+M2, data=dta)
options(huxtable.knitr_output_format="md")
jtools::export_summs(m1,m2,m3,m4,m5,
model.names = c("Y", "M1","M2","Y","Y"),
error_format = "[{conf.low},{conf.high}]")
## Registered S3 methods overwritten by 'broom':
## method from
## tidy.glht jtools
## tidy.summary.glht jtools
## Warning in to_md.huxtable(structure(list(names = c("", "(Intercept)", "", :
## Markdown cannot handle cells with colspan/rowspan > 1
## Warning in to_md.huxtable(structure(list(names = c("", "(Intercept)", "", :
## Can't vary column alignment in markdown; using first row
| (Intercept) |
19.17 *** |
8.19 *** |
7.51 *** |
7.52 *** |
10.47 *** |
|
[17.84,20.49] |
[7.32,9.06] |
[6.57,8.46] |
[6.55,8.49] |
[8.72,12.23] |
| X |
0.59 *** |
0.58 *** |
|
-0.00 |
|
|
[0.47,0.70] |
[0.51,0.65] |
|
[-0.08,0.08] |
|
| M1 |
|
|
0.67 *** |
0.68 *** |
|
|
|
|
[0.61,0.74] |
[0.60,0.75] |
|
| M2 |
|
|
|
|
0.87 *** |
|
|
|
|
|
[0.77,0.97] |
| N |
500 |
500 |
500 |
500 |
500 |
| R2 |
0.17 |
0.32 |
0.48 |
0.48 |
0.38 |
|
*** p < 0. |
001; ** p < |
0.01; * p < |
0.05. |
|
jtools::export_summs(m6,m7,m8,m9,
model.names = c("Y", "Y","Y","Y"),
error_format = "[{conf.low},{conf.high}]")
## Warning in to_md.huxtable(structure(list(names = c("", "(Intercept)", "", :
## Markdown cannot handle cells with colspan/rowspan > 1
## Warning in to_md.huxtable(structure(list(names = c("", "(Intercept)", "", :
## Can't vary column alignment in markdown; using first row
| (Intercept) |
12.40 *** |
9.13 *** |
9.32 *** |
8.88 *** |
|
[11.01,13.80] |
[7.51,10.76] |
[7.56,11.07] |
[7.23,10.53] |
| M1 |
0.91 *** |
0.61 *** |
|
0.56 *** |
|
[0.81,1.00] |
[0.49,0.73] |
|
[0.43,0.70] |
| M2 |
|
0.44 *** |
0.75 *** |
0.44 *** |
|
|
[0.31,0.56] |
[0.65,0.86] |
[0.31,0.56] |
| X |
|
|
0.29 *** |
0.09 |
|
|
|
[0.19,0.40] |
[-0.02,0.20] |
| N |
500 |
500 |
500 |
500 |
| R2 |
0.43 |
0.48 |
0.41 |
0.48 |
| ** |
* p < 0.001; * |
* p < 0.01; * |
p < 0.05. |
|
## 以 SEM 中的徑路分析方式，分析中介
model1 <- '
M1 ~ a*X
M2 ~ b2*X+b1*M1
Y ~ c1*X+c2*M1+c3*M2
indirectXM1M2 := a*b1*c2
indirectXM2 := b2*c3
totalX := c1 + (a*b1*c2) + (b2*c3)
proportion := (indirectXM1M2+indirectXM2)/totalX
'
#徑路分析報表
fit <- lavaan::sem(model1, data=dta)
summary(fit)
## lavaan 0.6.15 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 9
##
## Number of observations 500
##
## 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|)
## M1 ~
## X (a) 0.580 0.038 15.267 0.000
## M2 ~
## X (b2) -0.002 0.039 -0.060 0.952
## M1 (b1) 0.676 0.038 17.650 0.000
## Y ~
## X (c1) 0.092 0.055 1.658 0.097
## M1 (c2) 0.561 0.068 8.195 0.000
## M2 (c3) 0.436 0.063 6.940 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .M1 17.409 1.101 15.811 0.000
## .M2 12.765 0.807 15.811 0.000
## .Y 25.156 1.591 15.811 0.000
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|)
## indirectXM1M2 0.220 0.033 6.683 0.000
## indirectXM2 -0.001 0.017 -0.060 0.952
## totalX 0.311 0.054 5.789 0.000
## proportion 0.705 0.140 5.034 0.000
#以拔靴法看徑路係數與中介效果信賴區間
set.seed(1234)
fit <- lavaan::sem(model1, data=dta, test="bootstrap", bootstrap=501)
## Warning in lav_model_test(lavmodel = lavmodel, lavpartable = lavpartable, :
## lavaan WARNING: 2 bootstrap runs failed or did not converge.
parameterEstimates(fit,ci=TRUE,boot.ci.type="bca.simple")
## lhs op rhs label est
## 1 M1 ~ X a 0.580
## 2 M2 ~ X b2 -0.002
## 3 M2 ~ M1 b1 0.676
## 4 Y ~ X c1 0.092
## 5 Y ~ M1 c2 0.561
## 6 Y ~ M2 c3 0.436
## 7 M1 ~~ M1 17.409
## 8 M2 ~~ M2 12.765
## 9 Y ~~ Y 25.156
## 10 X ~~ X 24.128
## 11 indirectXM1M2 := a*b1*c2 indirectXM1M2 0.220
## 12 indirectXM2 := b2*c3 indirectXM2 -0.001
## 13 totalX := c1+(a*b1*c2)+(b2*c3) totalX 0.311
## 14 proportion := (indirectXM1M2+indirectXM2)/totalX proportion 0.705
## se z pvalue ci.lower ci.upper
## 1 0.038 15.267 0.000 0.505 0.654
## 2 0.039 -0.060 0.952 -0.080 0.075
## 3 0.038 17.650 0.000 0.601 0.751
## 4 0.055 1.658 0.097 -0.017 0.200
## 5 0.068 8.195 0.000 0.427 0.695
## 6 0.063 6.940 0.000 0.313 0.559
## 7 1.101 15.811 0.000 15.251 19.567
## 8 0.807 15.811 0.000 11.183 14.347
## 9 1.591 15.811 0.000 22.037 28.274
## 10 0.000 NA NA 24.128 24.128
## 11 0.033 6.683 0.000 0.155 0.285
## 12 0.017 -0.060 0.952 -0.035 0.033
## 13 0.054 5.789 0.000 0.205 0.416
## 14 0.140 5.034 0.000 0.430 0.979
#畫圖看模型與估計值
lavaanPlot::lavaanPlot(model = fit,
edge_options = list(color = "grey"),
coefs = TRUE,
stand = TRUE)