中介效果

一個獨變項，一個中介變項

#Baron & Kenny（1986） 的四步驟，檢驗被霸凌是體指數到憂鬱的中介
#程式報表5.3
m1 <- lm(憂鬱 ~ 體指數, data=dta)
m2 <- lm(被霸凌 ~ 體指數, data=dta) 
m3 <- lm(憂鬱 ~ 被霸凌, data=dta)
m4 <- lm(憂鬱 ~ 被霸凌 + 體指數, data=dta)

options(huxtable.knitr_output_format="md")
jtools::export_summs(m1,m2,m3,m4,
                     model.names = c("憂鬱", "被霸凌","憂鬱","憂鬱"),
                     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)	12.31 ***	1.41 ***	11.70 ***	9.53 ***
	[9.85,14.77]	[0.91,1.90]	[11.13,12.26]	[7.26,11.81]
體指數	0.17 **	0.03 *		0.11
	[0.05,0.28]	[0.01,0.05]		[-0.00,0.21]
被霸凌			1.99 ***	1.98 ***
			[1.79,2.19]	[1.78,2.18]
N	2000	2000	2000	2000
R2	0.00	0.00	0.16	0.16
**	* p < 0.001;	** p < 0.01;	* p < 0.05.

## 以 SEM 中的徑路分析方式，分析中介
model1 <- '
  憂鬱 ~ c*體指數+b*被霸凌
  被霸凌 ~ a*體指數 
  indirect := a*b
  total := c + (a*b)
  proportion := indirect/total
'

#徑路分析報表
#程式報表5.5
fit <- lavaan::sem(model1, data=dta)
summary(fit)

## lavaan 0.6.15 ended normally after 8 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                         5
## 
##   Number of observations                          2000
## 
## 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|)
##   憂鬱 ~                                           
##     體指數    c    0.106    0.055    1.928    0.054
##     被霸凌    b    1.978    0.102   19.444    0.000
##   被霸凌 ~                                         
##     體指數    a    0.030    0.012    2.504    0.012
## 
## Variances:
##                  Estimate  Std.Err  z-value  P(>|z|)
##    .憂鬱           83.081    2.627   31.623    0.000
##    .被霸凌          4.013    0.127   31.623    0.000
## 
## Defined Parameters:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     indirect          0.060    0.024    2.484    0.013
##     total             0.166    0.060    2.769    0.006
##     proportion        0.361    0.152    2.373    0.018

#Sobel test
a <- 0.030    
ase <- 0.012
b <- 1.978
bse <- 0.102

#計算中介效果與標準誤，並進行檢驗
ab <- a*b
abse <- sqrt(a^2 * bse^2 + b^2 * ase^2)
c(ab=ab, z = ab/abse, p = 2 * (1 - pnorm(abs(ab/abse))))

##         ab          z          p 
## 0.05934000 2.47948058 0.01315739

#以拔靴法看徑路係數與中介效果信賴區間
#程式報表5.6
set.seed(1234)
fit <- lavaan::sem(model1, data=dta, test="bootstrap", bootstrap=501)

## Warning in lav_model_test(lavmodel = lavmodel, lavpartable = lavpartable, :
## lavaan WARNING: 5 bootstrap runs failed or did not converge.

parameterEstimates(fit,ci=TRUE,boot.ci.type="bca.simple")

##          lhs op            rhs      label    est    se      z pvalue ci.lower
## 1       憂鬱  ~         體指數          c  0.106 0.055  1.928  0.054   -0.002
## 2       憂鬱  ~         被霸凌          b  1.978 0.102 19.444  0.000    1.779
## 3     被霸凌  ~         體指數          a  0.030 0.012  2.504  0.012    0.007
## 4       憂鬱 ~~           憂鬱            83.081 2.627 31.623  0.000   77.932
## 5     被霸凌 ~~         被霸凌             4.013 0.127 31.623  0.000    3.764
## 6     體指數 ~~         體指數            13.724 0.000     NA     NA   13.724
## 7   indirect :=            a*b   indirect  0.060 0.024  2.484  0.013    0.013
## 8      total :=        c+(a*b)      total  0.166 0.060  2.769  0.006    0.049
## 9 proportion := indirect/total proportion  0.361 0.152  2.373  0.018    0.063
##   ci.upper
## 1    0.214
## 2    2.178
## 3    0.054
## 4   88.230
## 5    4.261
## 6   13.724
## 7    0.107
## 8    0.284
## 9    0.658

#畫圖看模型與估計值
#圖5.5
lavaanPlot::lavaanPlot(model = fit,
                       edge_options = list(color = "grey"), 
                       coefs = TRUE,
                       stand = TRUE)

中介分析

資料與管理

中介效果