Bad (non-robust) standard-errors using "bootstrap" => CI's and p-values do not match

[LGraz]

Dear Team,
I appreciate this great product, THANKS!

I encountered a case where the confidence interval of or_fraction_mediated does not include 0, but the p-value is 0.9

I noticed that the SE estimate from the bootstrap is corrupted by outlier from the bootstraps (Hence the p-value might be too big). If thats so, I suggest using some more robust version (below I compare the SE with the MAD).

Here my minimal reproducible example using the data_WV_scaled.csv :

library(lavaan)

D <- read.csv("./data_WV_scaled.csv")

sem_mod_def <- {
"# Mediator
    F ~ df*D + of*O
    A ~ oa*O + da*D
    R ~ fr*F + ar*A
# direct effect
    R ~ dr*D + or*O
# Mediated effects
    dr_mediated := df*fr + da*ar
    or_mediated := of*fr + oa*ar
# total effect (direct + mediated)
    dr_total := dr + dr_mediated
    or_total := or + or_mediated
# mediated_effect / total_effect:
    dr_fraction_mediated := dr_mediated / dr_total
    or_fraction_mediated := or_mediated / or_total
"
}


# fit SEM and return parameters
set.seed(123)
object <- sem(sem_mod_def, D, se = "boot", # else "robust.sem"?
bootstrap = 2000)

BOOT <- lavaan:::lav_object_inspect_boot(object)
BOOT.def <- apply(BOOT, 1, object@[email protected])
or_fraction_mediated <- BOOT.def[6,]
mad(or_fraction_mediated)
#> [1] 0.38009
sd(or_fraction_mediated)
#> [1] 11.44
# ----> THATS A BIG DIFFERENCE !!!

parameterEstimates(object)
#>                     lhs op                  rhs                label    est     se      z pvalue ci.lower ci.upper
#> 1                     F  ~                    D                   df  0.411  0.023 17.665  0.000    0.362    0.453
#> 2                     F  ~                    O                   of -0.101  0.025 -4.021  0.000   -0.150   -0.051
#> 3                     A  ~                    O                   oa  0.479  0.025 19.148  0.000    0.431    0.528
#> 4                     A  ~                    D                   da -0.044  0.024 -1.839  0.066   -0.091    0.006
#> 5                     R  ~                    F                   fr  0.573  0.029 19.600  0.000    0.514    0.632
#> 6                     R  ~                    A                   ar -0.140  0.025 -5.534  0.000   -0.188   -0.090
#> 7                     R  ~                    D                   dr -0.024  0.024 -1.015  0.310   -0.069    0.023
#> 8                     R  ~                    O                   or  0.037  0.025  1.492  0.136   -0.011    0.085
#> 9                     F ~~                    F                       0.803  0.033 24.295  0.000    0.737    0.867
#> 10                    A ~~                    A                       0.759  0.025 30.706  0.000    0.708    0.806
#> 11                    R ~~                    R                       0.624  0.026 23.794  0.000    0.572    0.674
#> 12                    D ~~                    D                       0.999  0.000     NA     NA    0.999    0.999
#> 13                    D ~~                    O                      -0.206  0.000     NA     NA   -0.206   -0.206
#> 14                    O ~~                    O                       0.999  0.000     NA     NA    0.999    0.999
#> 15          dr_mediated :=          df*fr+da*ar          dr_mediated  0.242  0.019 12.966  0.000    0.204    0.278
#> 16          or_mediated :=          of*fr+oa*ar          or_mediated -0.125  0.020 -6.231  0.000   -0.163   -0.085
#> 17             dr_total :=       dr+dr_mediated             dr_total  0.218  0.025  8.742  0.000    0.168    0.267
#> 18             or_total :=       or+or_mediated             or_total -0.088  0.028 -3.161  0.002   -0.141   -0.033
#> 19 dr_fraction_mediated := dr_mediated/dr_total dr_fraction_mediated  1.110  0.122  9.118  0.000    0.907    1.380
#> 20 or_fraction_mediated := or_mediated/or_total or_fraction_mediated  1.422 11.440  0.124  0.901    0.909    3.317

Unrelated to this robustness issue: Would you recommend to me using se = "robust.sem" instead?

[yrosseel]

I would not recommend se = "robust.sem" here, as you define new variables that are clearly nonlinear functions of the original parameters.
I also hesitate to just replace 'sd' by 'mad' while calculating the bootstrap-based standard errors. Surely, sd() is very sensitive to outliers, but if an outlier effectively occurs every now and then, then this should also be reflected in the standard errors. I would have to dive into the literature to better understand the theoretical implications if mad() is used instead of sd() while computing bootstrap based standard errors.
I would rather try to pinpoint when/why you often get outlying parameter estimates. Perhaps using bounded estimation may prevent this? (Try running this with the argument bounds = "standard").

[LGraz]

bounds = "standard" does not have an effect.

pinpoint why outlying parameter:
or_fraction_mediated := or_mediated/or_total so if or_total is very close to 0 in one sample, then or_fraction_mediated explodes and hence also sd breaks down.

Simple exchanging sd() with mad() is not advisable.

What do you think about a warning if the sd and the mad are far away from each other? like:

# BOOT <- lavaan:::lav_object_inspect_boot(object)
# BOOT.def <- apply(BOOT, 1, object@[email protected])
sd_mad_ratio <- apply(BOOT.def, 1, sd) / apply(BOOT.def, 1, mad)
if(any(5 > sd_mad_ratio)){
    params_w_outliers <- names(sd_mad_ratio)[5 < sd_mad_ratio]
    warning(paste("The following bootstrap-parameters have a high
      ratio of standard deviation to median absolute deviation:", 
      params_w_outliers, 
      "\nP-values and Confidence Intervals might not match."))
}