Compute the difference between any two effects and form its percentile bootstrap confidence interval

[sfcheung]

Some users would like to compute and test the difference between any two effects.

For now, this cannot be done because the outputs of indirect_effect() are supposed to have one x-variable and one y-variable each. The + and - operator, when used on the outputs, will also be a similar object with one x-variable and one y-variable. If we allow the operators to work on two paths with different x and/or y variable(s), we need to make drastic changes to the indirect-class object.

Nevertheless, there are cases in which this is what users need.

Before we have decided on a long term solution, the following ad hoc function, any_diff_boot, can be used this way. It compute the difference between two effects, that of object1 - that of object2, and form the bootstrap CIs, if they have been requested when computing the two effects. Note: Users need to ensure that they are computed on the same set of bootstrap samples, which should be the case if the same object is used in the boot_out argument.

any_diff_boot <- function(object1,
                          object2,
                          level = .95) {
    boot_est1 <- object1$boot_indirect
    boot_est2 <- object2$boot_indirect
    if (is.null(boot_est1) || is.null(boot_est2)) {
        stop("Bootstrap estimates not available one or both objects.")
      }
    nboot <- length(boot_est1)
    if (nboot != length(boot_est2)) {
        stop("The numbers of bootstrap samples are not identical.")
      }
    boot_diff <- boot_est1 - boot_est2
    diff <- object1$indirect - object2$indirect
    tmp <- list(t = matrix(boot_diff, nrow = nboot, ncol = 1),
                t0 = diff,
                R = nboot)
    boot_ci0 <- boot::boot.ci(tmp, conf = level, type = "perc")
    boot_ci1 <- boot_ci0$percent[4:5]
    names(boot_ci1) <- paste0(formatC(c(100 * (1 - level) / 2,
                                  100 * (1 - (1 - level) / 2)), 2,
                                  format = "f"), "%")
    boot_se <- stats::sd(boot_diff, na.rm = TRUE)
    out <- c(Difference = diff, boot_ci1)
    out
  }

The following illustrates how it is used:

# Source:
# https://sfcheung.github.io/manymome/articles/manymome.html#mediation-only
library(manymome)
library(lavaan)
#> This is lavaan 0.6-16
#> lavaan is FREE software! Please report any bugs.
library(semhelpinghands)

dat <- data_serial

mod_med <- "
m1 ~ x + c1 + c2
m2 ~ m1 + x + c1 + c2
y ~ m2 + m1 + bx*x + bc1*c1 + c2
diff := bx - bc1
"
fit_med <- sem(model = mod_med,
               data = dat,
               fixed.x = FALSE)
fit_med_boot <- sem(model = mod_med,
                    data = dat,
                    fixed.x = FALSE,
                    se = "boot",
                    bootstrap = 1000,
                    iseed = 12345)
boot_out <- do_boot(fit_med,
                    R = 1000,
                    seed = 12345)
#> 5 processes started to run bootstrapping.
#> The expected CPU time is about 0 second(s).
tmp <- indirect_effect(x = "x",
                       y = "y",
                       fit = fit_med,
                       standardized_x = TRUE,
                       standardized_y = TRUE)
std_x_y <- indirect_effect(x = "x",
                           y = "y",
                           fit = fit_med,
                           standardized_x = TRUE,
                           standardized_y = TRUE,
                           boot_ci = TRUE,
                           boot_out = boot_out)
std_c1_y <- indirect_effect(x = "c1",
                            y = "y",
                            fit = fit_med,
                            standardized_x = TRUE,
                            standardized_y = TRUE,
                            boot_ci = TRUE,
                            boot_out = boot_out)
std_x_y
#> 
#> ==  Effect ==
#>                                       
#>  Path:               x -> y           
#>  Effect              0.234            
#>  95.0% Bootstrap CI: [-0.017 to 0.468]
#> 
#> Computation Formula:
#>   (b.y~x)*sd_x/sd_y
#> Computation:
#>   (0.49285)*(0.95010)/(1.99960)
#> 
#> Percentile confidence interval formed by nonparametric bootstrapping
#> with 1000 bootstrap samples.
std_c1_y
#> 
#> ==  Effect ==
#>                                       
#>  Path:               c1 -> y          
#>  Effect              0.049            
#>  95.0% Bootstrap CI: [-0.157 to 0.245]
#> 
#> Computation Formula:
#>   (b.y~c1)*sd_c1/sd_y
#> Computation:
#>   (0.09884)*(0.99841)/(1.99960)
#> 
#> Percentile confidence interval formed by nonparametric bootstrapping
#> with 1000 bootstrap samples.

# Show that the estimates are the same
print(standardizedSolution(fit_med)[10:11, 1:5], nd = 7)
#>    lhs op rhs label   est.std
#> 10   y  ~   x    bx 0.2341759
#> 11   y  ~  c1   bc1 0.0493495
coef(std_x_y)
#>       y~x 
#> 0.2341759
coef(std_c1_y)
#>       y~c1 
#> 0.04934946

# Show that the differences computed are the same
print(standardizedSolution(fit_med)[22, 1:5], nd = 8)
#>     lhs op    rhs label   est.std
#> 22 diff := bx-bc1  diff 0.1848264
coef(std_x_y) - coef(std_c1_y)
#>       y~x 
#> 0.1848264
any_diff_boot(std_x_y, std_c1_y)
#> Difference      2.50%     97.50% 
#>  0.1848264 -0.1042021  0.4646119

# Show that the bootstrap CIs are the same
std_boot_ci <- standardizedSolution_boot_ci(fit_med_boot)
print(std_boot_ci[22, c(1:5, 11, 12)], nd = 8)
#>     lhs op    rhs label   est.std boot.ci.lower boot.ci.upper
#> 22 diff := bx-bc1  diff 0.1848264    -0.1042021     0.4646119
any_diff_boot(std_x_y, std_c1_y)
#> Difference      2.50%     97.50% 
#>  0.1848264 -0.1042021  0.4646119

^{Created on 2023-11-18 with reprex v2.0.2}

Certainly, if the goal is to test the difference between two standardized indirect effects, lavaan can already do this: label the two paths, define the difference, and form the bootstrap CI for the standardized solution (this needs to be done using functions such as semhelpinghands::standardizedSolution_boot_ci() because the CIs in lavaan::standardizedSolution() are not bootstrap CIs even if bootstrapping is requested).

Nevertheless, if do_boot() is used, then indirect_effect() and do_boot() can be used, without the need to do bootstrapping when fitting the model by lavaan::sem().