中介效果
一個獨變項，一個中介變項
#Baron & Kenny（1986） 的四步驟，檢驗被霸凌是體指數到憂鬱的中介
m1 <- lm(Y ~ X, data=dta)
m2 <- lm(M1 ~ X, data=dta)
m3 <- lm(M2 ~ X, data=dta)
m4 <- lm(Y ~ M1, data=dta)
m5 <- lm(Y ~ M2, data=dta)
m6 <- lm(Y ~ M1+M2, data=dta)
m7 <- lm(Y ~ X+M1, data=dta)
m8 <- lm(Y ~ X+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) |
16.99 *** |
7.67 *** |
7.55 *** |
13.46 *** |
12.84 *** |
|
[15.91,18.06] |
[6.93,8.40] |
[6.88,8.22] |
[12.04,14.89] |
[11.37,14.31] |
| X |
0.69 *** |
0.64 *** |
0.68 *** |
|
|
|
[0.59,0.78] |
[0.58,0.71] |
[0.62,0.74] |
|
|
| M1 |
|
|
|
0.73 *** |
|
|
|
|
|
[0.64,0.83] |
|
| M2 |
|
|
|
|
0.76 *** |
|
|
|
|
|
[0.67,0.86] |
| N |
500 |
500 |
500 |
500 |
500 |
| R2 |
0.28 |
0.42 |
0.50 |
0.31 |
0.32 |
|
*** 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) |
9.21 *** |
13.21 *** |
13.04 *** |
9.26 *** |
|
[7.69,10.74] |
[11.83,14.59] |
[11.60,14.48] |
[7.66,10.87] |
| M1 |
0.50 *** |
0.49 *** |
|
0.49 *** |
|
[0.40,0.60] |
[0.37,0.61] |
|
[0.38,0.61] |
| M2 |
0.53 *** |
|
0.52 *** |
0.52 *** |
|
[0.43,0.63] |
|
[0.39,0.66] |
[0.40,0.65] |
| X |
|
0.37 *** |
0.33 *** |
0.01 |
|
|
[0.25,0.49] |
[0.20,0.46] |
[-0.13,0.16] |
| N |
500 |
500 |
500 |
500 |
| R2 |
0.44 |
0.36 |
0.35 |
0.44 |
| * |
** p < 0.001; |
** p < 0.01; |
* p < 0.05. |
|
## 以 SEM 中的徑路分析方式，分析中介
model1 <- '
Y ~ c*X+b1*M1+b2*M2
M1 ~ a1*X
M2 ~ a2*X
indirect1 := a1*b1
indirect2 := a2*b2
total := c + (a1*b1) + (a2*b2)
proportion := (indirect1+indirect2)/total
'
#徑路分析報表
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 8
##
## Number of observations 500
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 1
## P-value (Chi-square) 0.993
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## Y ~
## X (c) 0.014 0.072 0.189 0.850
## M1 (b1) 0.492 0.058 8.518 0.000
## M2 (b2) 0.523 0.064 8.186 0.000
## M1 ~
## X (a1) 0.644 0.034 19.037 0.000
## M2 ~
## X (a2) 0.681 0.031 22.271 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .Y 25.153 1.591 15.811 0.000
## .M1 15.063 0.953 15.811 0.000
## .M2 12.325 0.780 15.811 0.000
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|)
## indirect1 0.317 0.041 7.775 0.000
## indirect2 0.356 0.046 7.683 0.000
## total 0.687 0.049 13.896 0.000
## proportion 0.980 0.104 9.414 0.000
#以拔靴法看徑路係數與中介效果信賴區間
set.seed(1234)
fit <- lavaan::sem(model1, data=dta, test="bootstrap", bootstrap=501)
parameterEstimates(fit,ci=TRUE,boot.ci.type="bca.simple")
## lhs op rhs label est se z
## 1 Y ~ X c 0.014 0.072 0.189
## 2 Y ~ M1 b1 0.492 0.058 8.518
## 3 Y ~ M2 b2 0.523 0.064 8.186
## 4 M1 ~ X a1 0.644 0.034 19.037
## 5 M2 ~ X a2 0.681 0.031 22.271
## 6 Y ~~ Y 25.153 1.591 15.811
## 7 M1 ~~ M1 15.063 0.953 15.811
## 8 M2 ~~ M2 12.325 0.780 15.811
## 9 X ~~ X 26.358 0.000 NA
## 10 indirect1 := a1*b1 indirect1 0.317 0.041 7.775
## 11 indirect2 := a2*b2 indirect2 0.356 0.046 7.683
## 12 total := c+(a1*b1)+(a2*b2) total 0.687 0.049 13.896
## 13 proportion := (indirect1+indirect2)/total proportion 0.980 0.104 9.414
## pvalue ci.lower ci.upper
## 1 0.85 -0.128 0.155
## 2 0.00 0.379 0.606
## 3 0.00 0.398 0.648
## 4 0.00 0.577 0.710
## 5 0.00 0.621 0.741
## 6 0.00 22.035 28.271
## 7 0.00 13.196 16.930
## 8 0.00 10.797 13.853
## 9 NA 26.358 26.358
## 10 0.00 0.237 0.397
## 11 0.00 0.265 0.447
## 12 0.00 0.590 0.783
## 13 0.00 0.776 1.184
#畫圖看模型與估計值
lavaanPlot::lavaanPlot(model = fit,
edge_options = list(color = "grey"),
coefs = TRUE,
stand = TRUE)