Step 0_A priori power analysis

library(lavaan)

This is lavaan 0.6-21
lavaan is FREE software! Please report any bugs.

library(MASS)
library(simsem)

#################################################################

This is simsem 0.5-17

simsem is BETA software! Please report any bugs.

simsem was first developed at the University of Kansas Center for

Research Methods and Data Analysis, under NSF Grant 1053160.

#################################################################

Attaching package: 'simsem'

The following object is masked from 'package:lavaan':

    inspect

Alternative a) Monte Carlo simulation based on Lavaan

1.Enter effect values for each path —-

# --- TRUE POPULATION VALUES
true_params <- list(
  # Main effects of dummies on IM
  b_r1e1 = 0.25, #indirect effects
  b_r0e2 = 0.25,
  b_r1e2 = 0.25,
  b_r0e3 = 0.25,
  b_r1e3 = 0.25,
  
  # Main effects of TR, ADT, and PEB on IM
  b_TR  = 0.25,
  b_ADT = 0.25,
  b_PEB = 0.25,
  
  # EEF regressions
  c_r1e1 = 0.37,
  c_r0e2 = 0.37,
  c_r1e2 = 0.37,
  c_r0e3 = 0.37,
  c_r1e3 = 0.37,
  c_IM   = 0.57, #indirect effects
  c_TR   = 0.20,
  c_ADT  = 0.20,
  c_PEB  = 0.40,
  
  # EEC regressions
  d_r1e1 = 0.37,
  d_r0e2 = 0.37,
  d_r1e2 = 0.37,
  d_r0e3 = 0.37,
  d_r1e3 = 0.37,
  d_IM   = 0.57, #indirect effects
  d_TR   = 0.20,
  d_ADT  = 0.20,
  d_PEB  = 0.40,
  
  # Residual covariance/correlation between EEF and EEC
  d_EEF_EEC = 0.35
)

2. Simulate dataset —-

simulate_one_dataset <- function(N, true_params) {
  
  with(true_params, {
    
    # Exogenous continuous variables
    TR      <- rnorm(N, mean = 0, sd = 1)
    ADT     <- rnorm(N, mean = 0, sd = 1)
    PEB_yes <- rnorm(N, mean = 0, sd = 1)
    
    # Binary predictors
    r1e1 <- rbinom(N, size = 1, prob = 0.5)
    r0e2 <- rbinom(N, size = 1, prob = 0.5)
    r1e2 <- rbinom(N, size = 1, prob = 0.5)
    r0e3 <- rbinom(N, size = 1, prob = 0.5)
    r1e3 <- rbinom(N, size = 1, prob = 0.5)
    
    # IM
    lin_IM <-
      b_r1e1 * r1e1 +
      b_r0e2 * r0e2 +
      b_r1e2 * r1e2 +
      b_r0e3 * r0e3 +
      b_r1e3 * r1e3 +
      b_TR   * TR +
      b_ADT  * ADT +
      b_PEB  * PEB_yes
    
    IM <- lin_IM + rnorm(N, mean = 0, sd = 1)
    
    # EEF
    lin_EEF <-
      c_r1e1 * r1e1 +
      c_r0e2 * r0e2 +
      c_r1e2 * r1e2 +
      c_r0e3 * r0e3 +
      c_r1e3 * r1e3 +
      c_IM   * IM +
      c_TR   * TR +
      c_ADT  * ADT +
      c_PEB  * PEB_yes
    
    # EEC
    lin_EEC <-
      d_r1e1 * r1e1 +
      d_r0e2 * r0e2 +
      d_r1e2 * r1e2 +
      d_r0e3 * r0e3 +
      d_r1e3 * r1e3 +
      d_IM   * IM +
      d_TR   * TR +
      d_ADT  * ADT +
      d_PEB  * PEB_yes
    
    # Correlated residuals for EEF and EEC
    residuals <- MASS::mvrnorm(
      n = N,
      mu = c(0, 0),
      Sigma = matrix(
        c(
          1,           d_EEF_EEC,
          d_EEF_EEC,  1
        ),
        nrow = 2,
        byrow = TRUE
      )
    )
    
    # Outcomes
    EEF <- lin_EEF + residuals[, 1]
    EEC <- lin_EEC + residuals[, 2]
    
    # Return dataset
    dat <- data.frame(
      r1e1, r0e2, r1e2, r0e3, r1e3,
      IM, EEF, EEC, TR, ADT, PEB_yes
    )
    
    dat
  })
}

3. SEM model

MediationModel <- '
  #### structural regressions
  IM ~ b_r1e1*r1e1 +
       b_r0e2*r0e2 +
       b_r1e2*r1e2 +
       b_r0e3*r0e3 +
       b_r1e3*r1e3 +
       b_TR*TR +
       b_ADT*ADT +
       b_PEB*PEB_yes

  EEF ~ c_r1e1*r1e1 +
        c_r0e2*r0e2 +
        c_r1e2*r1e2 +
        c_r0e3*r0e3 +
        c_r1e3*r1e3 +
        c_IM*IM +
        c_TR*TR +
        c_ADT*ADT +
        c_PEB*PEB_yes

  EEC ~ d_r1e1*r1e1 +
        d_r0e2*r0e2 +
        d_r1e2*r1e2 +
        d_r0e3*r0e3 +
        d_r1e3*r1e3 +
        d_IM*IM +
        d_TR*TR +
        d_ADT*ADT +
        d_PEB*PEB_yes
      
        EEF ~~ EEC
        
  # Indirect effects through IM -> EEF
  ind_r1e1_EEF := b_r1e1*c_IM
  ind_r0e2_EEF := b_r0e2*c_IM
  ind_r1e2_EEF := b_r1e2*c_IM
  ind_r0e3_EEF := b_r0e3*c_IM
  ind_r1e3_EEF := b_r1e3*c_IM

  # Indirect effects through IM -> EEC
  ind_r1e1_EEC := b_r1e1*d_IM
  ind_r0e2_EEC := b_r0e2*d_IM
  ind_r1e2_EEC := b_r1e2*d_IM
  ind_r0e3_EEC := b_r0e3*d_IM
  ind_r1e3_EEC := b_r1e3*d_IM
'

4. Monte Carlo for mediation effects

estimate_power_for_mediation <- function(
    N, true_params, MediationModel,
    params_to_check = c("b_r1e1", "b_r0e2", "b_r1e2", "b_r0e3", "b_r1e3"),
    nRep = 200, alpha = 0.05
) {
  
  pvals_list <- lapply(params_to_check, function(x) numeric(nRep))
  converged  <- logical(nRep)
  
  for (i in 1:nRep) {
    dat <- simulate_one_dataset(N, true_params)
    
    fit <- try(
      sem(MediationModel, data = dat, estimator = "MLR"),
      silent = TRUE
    )
    
    if (inherits(fit, "try-error") || !lavInspect(fit, "converged")) {
      converged[i] <- FALSE
      for (k in seq_along(params_to_check)) pvals_list[[k]][i] <- NA
    } else {
      converged[i] <- TRUE
      pe <- parameterEstimates(fit)
      for (k in seq_along(params_to_check)) {
        lab <- params_to_check[k]
        
        row <- pe[
          pe$lhs == lab &
            pe$op == ":=",
        ]
        
        pvals_list[[k]][i] <-
          if (nrow(row) == 1) row$pvalue else NA_real_
      }
    }
    
    if (i %% 20 == 0) cat("N =", N, "- Rep", i, "of", nRep, "done\n")
  }
  
  power_vec <- numeric(length(params_to_check))
  names(power_vec) <- params_to_check
  
    for (k in seq_along(params_to_check)) {
    pvec <- pvals_list[[k]]
    ok <- !is.na(pvec)
    
    if (any(ok)) {
      power_vec[k] <- mean(pvec[ok] < alpha)
    } else {
      power_vec[k] <- NA_real_
    }
  }
  
  list(
    power      = power_vec,
    n_success  = sum(converged),
    n_total    = nRep
  )
}

5. Power analysis for different sample sizes —-

5a) Indirect effects

for (N in Ns) {
  res <- results[[as.character(N)]]
  cat("N =", N,
      "→ power (ind_r1e1_EEF) ≈", round(res$power["ind_r1e1_EEF"], 3),
      "; power (ind_r0e2_EEF) ≈", round(res$power["ind_r0e2_EEF"], 3),
      "; power (ind_r1e2_EEF) ≈", round(res$power["ind_r1e2_EEF"], 3),
      "; power (ind_r0e3_EEF) ≈", round(res$power["ind_r0e3_EEF"], 3),
      "; power (ind_r1e3_EEF) ≈", round(res$power["ind_r1e3_EEF"], 3),
      "; power (ind_r1e1_EEC) ≈", round(res$power["ind_r1e1_EEC"], 3),
      "; power (ind_r0e2_EEC) ≈", round(res$power["ind_r0e2_EEC"], 3),
      "; power (ind_r1e2_EEC) ≈", round(res$power["ind_r1e2_EEC"], 3),
      "; power (ind_r0e3_EEC) ≈", round(res$power["ind_r0e3_EEC"], 3),
      "; power (ind_r1e3_EEC) ≈", round(res$power["ind_r1e3_EEC"], 3),
      "\n")
}

N = 200 → power (ind_r1e1_EEF) ≈ 0.392 ; power (ind_r0e2_EEF) ≈ 0.409 ; power (ind_r1e2_EEF) ≈ 0.398 ; power (ind_r0e3_EEF) ≈ 0.426 ; power (ind_r1e3_EEF) ≈ 0.422 ; power (ind_r1e1_EEC) ≈ 0.39 ; power (ind_r0e2_EEC) ≈ 0.409 ; power (ind_r1e2_EEC) ≈ 0.4 ; power (ind_r0e3_EEC) ≈ 0.422 ; power (ind_r1e3_EEC) ≈ 0.418 
N = 250 → power (ind_r1e1_EEF) ≈ 0.495 ; power (ind_r0e2_EEF) ≈ 0.483 ; power (ind_r1e2_EEF) ≈ 0.486 ; power (ind_r0e3_EEF) ≈ 0.489 ; power (ind_r1e3_EEF) ≈ 0.473 ; power (ind_r1e1_EEC) ≈ 0.5 ; power (ind_r0e2_EEC) ≈ 0.484 ; power (ind_r1e2_EEC) ≈ 0.482 ; power (ind_r0e3_EEC) ≈ 0.49 ; power (ind_r1e3_EEC) ≈ 0.477 
N = 270 → power (ind_r1e1_EEF) ≈ 0.535 ; power (ind_r0e2_EEF) ≈ 0.501 ; power (ind_r1e2_EEF) ≈ 0.51 ; power (ind_r0e3_EEF) ≈ 0.522 ; power (ind_r1e3_EEF) ≈ 0.537 ; power (ind_r1e1_EEC) ≈ 0.538 ; power (ind_r0e2_EEC) ≈ 0.505 ; power (ind_r1e2_EEC) ≈ 0.518 ; power (ind_r0e3_EEC) ≈ 0.523 ; power (ind_r1e3_EEC) ≈ 0.536 

5b) Dependent variables

Alternative b) SIMSEM