Bias corrected CI for mediation with multiply imputed data (SEM.MI vs. FIML)?

[zrc2...@gmail.com]

Hello,

I am using the "sem.mi" function to run a series of mediation models for imputed data. Importantly, I used mice to impute the entire dataset first and ran models based on the imputed data. I hope to get the bias corrected ci's and was not sure how to get the complete set of results you'd get from "parameterEstimates". Thank you in advance for any help/insight you can offer.

(1) MI dataset

HSMiss <- HolzingerSwineford1939[ , c(paste("x", 1:9, sep = ""),"ageyr","agemo","school")]

impHS=mice(HolzingerSwineford1939, m=20, seed=123, print=TRUE)

imps=impHS$imputations

View(HolzingerSwineford1939

       + )

HS.model <- ' # direct effect

+ x9 ~ c*x3

+ # mediator

+ x6 ~ a*x3

+ x9 ~ b*x6

+ # indirect effect (a*b)

+ ab := a*b

+ # total effect

+ total := c + (a*b)

+ '

out2=sem.mi(HS.model, data=impHS)

summary(out2)

THIS IS ALL THAT WAS PRODUCED FOR MI DATASET

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

    lhs op     rhs label std.lv std.all std.nox

1    x9  ~      x3     c  0.094   0.106   0.094

2    x6  ~      x3     a  0.060   0.062   0.055

3    x9  ~      x6     b  0.033   0.036   0.036

4    x9 ~~      x9        0.472   0.465   0.465

5    x6 ~~      x6        0.606   0.506   0.506

6    x3 ~~      x3        0.669   0.524   0.669

7    ab :=     a*b    ab  0.002   0.002   0.002

8 total := c+(a*b) total  0.096   0.108   0.096

(2) FIML

targetFit=sem(HS.model, HSMiss, meanstructure=TRUE, std.lv=TRUE, missing="fiml",se="bootstrap", bootstrap=5000)

summary(targetFit, rsquare=TRUE, standardized=TRUE)

HERE, I GET CI's and other output along with it when I use the same command

> parameterEstimates(targetFit, boot.ci.type = "bca.simple", standardized = TRUE)

     lhs op     rhs label   est    se      z pvalue ci.lower ci.upper std.lv std.all std.nox

1     x9  ~      x3     c 0.266 0.047  5.651  0.000    0.176    0.361  0.266   0.298   0.264

2     x6  ~      x3     a 0.191 0.060  3.201  0.001    0.077    0.313  0.191   0.198   0.175

3     x9  ~      x6     b 0.143 0.055  2.603  0.009    0.037    0.253  0.143   0.155   0.155

4     x9 ~~      x9       0.882 0.071 12.499  0.000    0.762    1.047  0.882   0.869   0.869

5     x6 ~~      x6       1.150 0.105 10.904  0.000    0.962    1.378  1.150   0.961   0.961

6     x3 ~~      x3       1.275 0.000     NA     NA    1.275    1.275  1.275   1.000   1.275

7     x9 ~1               4.463 0.161 27.672  0.000    4.142    4.771  4.463   4.430   4.430

8     x6 ~1               1.755 0.128 13.686  0.000    1.496    2.004  1.755   1.604   1.604

9     x3 ~1               2.250 0.000     NA     NA    2.250    2.250  2.250   1.993   2.250

10    ab :=     a*b    ab 0.027 0.015  1.812  0.070    0.006    0.068  0.027   0.031   0.027

11 total := c+(a*b) total 0.293 0.049  6.043  0.000    0.201    0.391  0.293   0.329   0.291

[Mauricio Garnier-Villarreal]

I am not sure the lavaan.mi object can report the pooled bootstrap CI. Also, this would take a long time, as you would need to run the 2000 bootstrap (for example) for each imputed data set.

I would recommend to use the MonteCarlo sampling method, here is an example on how to do this with MI

# impose missing data for example

HSMiss <- HolzingerSwineford1939[ , c(paste("x", 1:9, sep = ""),
                                      "ageyr","agemo","school")]

set.seed(12345)
HSMiss$x5 <- ifelse(HSMiss$x5 <= quantile(HSMiss$x5, .3), NA, HSMiss$x5)
age <- HSMiss$ageyr + HSMiss$agemo/12
HSMiss$x9 <- ifelse(age <= quantile(age, .3), NA, HSMiss$x9)

## specify CFA model from lavaan's ?cfa help page
HS.model <- '
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed   =~ x7 + x8 + x9

textual ~ a*visual
speed ~ cp*visual + b*textual
ind := a*b
'


## ... or impute missing data first
library(Amelia)
set.seed(12345)
HS.amelia <- amelia(HSMiss, m = 20, noms = "school", p2s = FALSE)
imps <- HS.amelia$imputations

out2 <- cfa.mi(HS.model, data = imps, std.lv=T, verbose=T)
sm <- data.frame(summary(out2, ci=T))

abs <- sm[sm$label %in%c("a","b") ,"est"] ## save the parameter estimates a and b

ACM <- vcov(out2, type="pooled")[c("a","b"),c("a","b")] # get the pooled coavriance matrix of the parameters

monteCarloMed('a*b', coef1=abs[1], coef2=abs[2], ACM=ACM) ## estimate the montecarlo sampling method