1-1-1 Mediation

讀進資料

資料 https://stats.oarc.ucla.edu/r/faq/how-can-i-perform-mediation-with-multilevel-data-method-2/

dta <- read.csv("mmedi1-1-1sim.csv",header=T)
head(dta)

##   id grp           x          m          y
## 1  1   1  1.54512054 0.10680165  0.5677760
## 2  2   1  2.27527200 2.11043937  1.2061451
## 3  3   1  0.78669920 0.03888269 -0.2612771
## 4  4   1 -0.06185872 0.47940645 -0.7589904
## 5  5   1  0.11725984 0.59082413  0.5190769
## 6  6   1  1.48092378 0.89094498 -0.6311193

載進套件

library(pacman)

## Warning: 套件 'pacman' 是用 R 版本 4.2.3 來建造的

p_load(devtools,lme4,reshape2,lmerTest)
install_github("marklhc/bootmlm")

## Skipping install of 'bootmlm' from a github remote, the SHA1 (bb5d477b) has not changed since last install.
##   Use `force = TRUE` to force installation

p_load(bootmlm)

Krull & MacKinnon (1999), 類似 Baro-Kenny method

rst1 <- lmer(y ~ x + (1 + x | grp), data = dta)
summary(rst1)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: y ~ x + (1 + x | grp)
##    Data: dta
## 
## REML criterion at convergence: 2435.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3156 -0.5552  0.0347  0.5920  2.9446 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  grp      (Intercept) 0.7470   0.8643        
##           x           0.2172   0.4661   -0.06
##  Residual             0.8226   0.9070        
## Number of obs: 800, groups:  grp, 100
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept) -0.02714    0.09598 99.57482  -0.283    0.778    
## x            0.68973    0.05865 99.87069  11.761   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##   (Intr)
## x -0.018

rst2 <- lmer(m ~ x + (1 + x | grp), data = dta)
summary(rst2)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: m ~ x + (1 + x | grp)
##    Data: dta
## 
## REML criterion at convergence: 2234.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6095 -0.6175  0.0217  0.6158  2.9143 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  grp      (Intercept) 0.7093   0.8422       
##           x           0.1167   0.3416   0.02
##  Residual             0.6450   0.8031       
## Number of obs: 800, groups:  grp, 100
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   0.09522    0.09162  99.60722   1.039    0.301    
## x             0.61145    0.04643 102.30830  13.169   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##   (Intr)
## x 0.049

rst3 <- lmer(y ~ m + x + (1 + m + x | grp), data = dta)
summary(rst3)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: y ~ m + x + (1 + m + x | grp)
##    Data: dta
## 
## REML criterion at convergence: 2044.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2655 -0.5846  0.0171  0.6039  2.3362 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr       
##  grp      (Intercept) 0.26147  0.5113              
##           m           0.12011  0.3466   -0.03      
##           x           0.03938  0.1984   -0.31 -0.04
##  Residual             0.50767  0.7125              
## Number of obs: 800, groups:  grp, 100
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept) -0.09374    0.06207 98.82689  -1.510    0.134    
## m            0.62507    0.04686 92.36026  13.339  < 2e-16 ***
## x            0.24971    0.03875 76.23636   6.444 9.49e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##   (Intr) m     
## m -0.045       
## x -0.084 -0.324

算出間接效果，並利用 Monte Carlo 法算出信賴區間

要先重排資料

newdta <- melt(dta, id.vars = c("id", "grp", "x", "m"),
   measure.vars = c("y", "m"), value.name = "z")
newdta <- within(newdta, {
  sy <- as.integer(variable == "y")
  sm <- as.integer(variable == "m")
})

dta[dta$id==1,]

##   id grp        x         m        y
## 1  1   1 1.545121 0.1068017 0.567776

newdta[newdta$id==1,]

##     id grp        x         m variable         z sm sy
## 1    1   1 1.545121 0.1068017        y 0.5677760  0  1
## 801  1   1 1.545121 0.1068017        m 0.1068017  1  0

有兩個 select variable, 用來切換依變項

把兩模型合一的模型

rst4 <- lmer(z ~ 0 + sm + sm:x + sy + sy:m + sy:x +
               (0 + sm + sm:x + sy + sy:m + sy:x | grp)+(0 + sm |id),
                data = newdta)
summary(rst4)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: z ~ 0 + sm + sm:x + sy + sy:m + sy:x + (0 + sm + sm:x + sy +  
##     sy:m + sy:x | grp) + (0 + sm | id)
##    Data: newdta
## 
## REML criterion at convergence: 4252.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3157 -0.5720  0.0282  0.5823  2.6167 
## 
## Random effects:
##  Groups   Name Variance Std.Dev. Corr                   
##  id       sm   0.13777  0.3712                          
##  grp      sm   0.67940  0.8243                          
##           sy   0.27030  0.5199    0.13                  
##           sm:x 0.12048  0.3471    0.06  0.07            
##           sy:m 0.11187  0.3345    0.03 -0.02  0.85      
##           x:sy 0.03244  0.1801   -0.04 -0.20 -0.34  0.09
##  Residual      0.50896  0.7134                          
## Number of obs: 1600, groups:  id, 800; grp, 100
## 
## Fixed effects:
##       Estimate Std. Error        df t value Pr(>|t|)    
## sm     0.09322    0.08943  99.02460   1.042    0.300    
## sy    -0.09685    0.06196  98.07399  -1.563    0.121    
## sm:x   0.61186    0.04649 101.26590  13.160   <2e-16 ***
## sy:m   0.61056    0.04554  92.08455  13.408   <2e-16 ***
## x:sy   0.22081    0.03725  71.15694   5.928    1e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      sm     sy     sm:x   sy:m  
## sy    0.104                     
## sm:x  0.078  0.044              
## sy:m  0.023 -0.040  0.465       
## x:sy -0.017 -0.026 -0.114 -0.285

算出間接效果，並利用 Monte Carlo 法算出信賴區間，來自http://www.quantpsy.org/medmc/medmc111.htm

擷取資訊部分來自 https://quantdev.ssri.psu.edu/sites/qdev/files/ILD_Ch07_2017_Within-PersonMedationWithMLM.html

FEmatrix <- fixef(rst4)
REmatrix <- as.data.frame(VarCorr(rst4))
ACM_FE <- vcov(rst4)
ACM_RE <- vcov_vc(rst4, sd_cor = FALSE, print_names = TRUE)


a <- FEmatrix['sm:x']
b <- FEmatrix['sy:m']
flag <- (REmatrix$var1=='sm:x' & REmatrix$var2=='sy:m')
flag[is.na(flag)] <- FALSE
covajbj <- REmatrix[flag,4]
Ind <- a*b+covajbj 
vara <- ACM_FE['sm:x','sm:x']
varb <- ACM_FE['sy:m','sy:m']
covab <- ACM_FE['sm:x','sy:m']
varcovajbj=ACM_RE['cov_sy:m.sm:x|grp','cov_sy:m.sm:x|grp']

rep=20000
conf=95
dvec=rnorm(rep)
avec=dvec*sqrt(vara)+a
bvec=dvec*covab/sqrt(vara)+sqrt(varb)*rnorm(rep,sd=sqrt(1-(covab^2)/(vara*varb)))+b
cvec=rnorm(rep)*sqrt(varcovajbj)+covajbj
ab=avec*bvec+cvec
low=(1-conf/100)/2
upp=((1-conf/100)/2)+(conf/100)
LL=quantile(ab,low)
UL=quantile(ab,upp)
LL4=format(LL,digits=4)
UL4=format(UL,digits=4)
Ind=format(Ind,digits=4)
c(IndirectEffeft = Ind, LowerLimit = LL4,UpperLimit=UL4)

## IndirectEffeft.sm:x     LowerLimit.2.5%    UpperLimit.97.5% 
##            "0.4725"            "0.3818"            "0.5739"

1-1-1 Mediation

CCP

2023-04-22

讀進資料

資料 https://stats.oarc.ucla.edu/r/faq/how-can-i-perform-mediation-with-multilevel-data-method-2/

載進套件

Krull & MacKinnon (1999), 類似 Baro-Kenny method

算出間接效果，並利用 Monte Carlo 法算出信賴區間

要先重排資料

把兩模型合一的模型

算出間接效果，並利用 Monte Carlo 法算出信賴區間，來自http://www.quantpsy.org/medmc/medmc111.htm

擷取資訊部分來自 https://quantdev.ssri.psu.edu/sites/qdev/files/ILD_Ch07_2017_Within-PersonMedationWithMLM.html