Mediation power analysis

The following example obtains sample size estimates for a simple mediation model to obtain a power of \(1 - \beta = .80\). Specifically, four sample sizes (\(N\)) are solved by varying the coefficients associated with the equations

1. \(M = a*X + e\)
2. \(Y = b*M + c^\prime * X + e\)

For simplicity the independent variable \(X\) is constrained to be a dichotomous variable, though this is not a requirement and users are free to change this distribution to suit their applications. Residuals (\(e\)) are assumed to be independent and distributed \(e \sim N(0,1)\). The simple mediation model itself is fitted using the lavaan package’s maximum-likelihood estimation criteria.

SimSolve code

library(SimDesign)

Design <- createDesign(N = NA,
                       a = sqrt(.35),
                       b = sqrt(c(.02, .07, .15, .35)),
                       cprime=.39)
Design    # solve for NA's

# A tibble: 4 × 4
      N       a       b cprime
  <dbl>   <dbl>   <dbl>  <dbl>
1    NA 0.59161 0.14142   0.39
2    NA 0.59161 0.26458   0.39
3    NA 0.59161 0.38730   0.39
4    NA 0.59161 0.59161   0.39

#~~~~~~~~~~~~~~~~~~~~~~~~
#### Step 2 --- Define generate, analyse, and summarise functions

Generate <- function(condition, fixed_objects = NULL) {
    Attach(condition)
    X <- rep(0:1, each=N)
    M <- a*X + rnorm(N)
    Y <- b*M + cprime*X + rnorm(N)
    dat <- data.frame(X, Y, M)
    dat
}

Analyse <- function(condition, dat, fixed_objects = NULL) {
    model <- ' # direct effect
             Y ~ c*X
           # mediator
             M ~ a*X
             Y ~ b*M
           # indirect effect (a*b)
             ab := a*b
           # total effect
             total := c + (a*b)
         '
    fit <- lavaan::sem(model, data=dat)
    if(!lavInspect(fit, 'converged')) stop('Model did not converge')
    PE <- parameterEstimates(fit)
    ret <- PE$pvalue[PE$lhs == 'ab']   # joint test
    ret
}

Summarise <- function(condition, results, fixed_objects = NULL) {
    ret <- EDR(results, alpha = .05)
    ret
}

Solve each \(N\) given the rows in Design

# Initial search between N = [10,500] for each row using the default
# integer solver (integer = TRUE)

# In this example, b = target power of 80%
interval <- c(20, 1000) # needless wide for most, but shows the point

# Search is terminated when either 100 iterations reached or the
# prediction CI is within [.795, .805]
solved <- SimSolve(design=Design, b=.8, interval=interval,
                   generate=Generate, analyse=Analyse, summarise=Summarise,
                   packages='lavaan', parallel=TRUE, 
                   ncores = ceiling(parallel::detectCores()/2),
                   verbose=FALSE, check.interval=FALSE,
                   maxiter=100, predCI.tol=.01)

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

solved

# A tibble: 4 × 4
        N       a       b cprime
    <dbl>   <dbl>   <dbl>  <dbl>
1 262.73  0.59161 0.14142   0.39
2  87.086 0.59161 0.26458   0.39
3  49.470 0.59161 0.38730   0.39
4  33.637 0.59161 0.59161   0.39

Additional information about the solutions

summary(solved)

$condition_1
$root
[1] 262.7264

$predCI.root
  CI_2.5  CI_97.5 
258.5502 266.9366 

$b
[1] 0.8

$predCI.b
[1] 0.7951403 0.8047727

$terminated_early
[1] TRUE

$time
[1] 02m 3.49s

$iterations
[1] 94

$total.replications
[1] 29900

$tab
           y   x reps
11 0.7516129 236  310
13 0.7863636 250  660
15 0.7823529 253  340
17 0.7701493 255  670
18 0.8152174 256  460
19 0.7974576 257 1180
20 0.7942308 258  520
21 0.7852459 259 1220
22 0.7913333 260 3000
23 0.7996721 261 3050
24 0.7983294 262 4190
25 0.8077143 263 3500
26 0.8130769 264 1300
27 0.8014388 265 1390
28 0.8216216 266  370
29 0.7924138 267 1450
30 0.7964029 268 1390
32 0.8413043 274  460
33 0.8187500 275  480
34 0.8318182 277  440
35 0.7900000 278  400
36 0.8023810 280  420
41 0.8219512 291  410


$condition_2
$root
[1] 87.08581

$predCI.root
  CI_2.5  CI_97.5 
86.12608 87.94639 

$b
[1] 0.8

$predCI.b
[1] 0.7950895 0.8048218

$terminated_early
[1] FALSE

$time
[1] 02m 1.36s

$iterations
[1] 100

$total.replications
[1] 32900

$tab
           y  x reps
5  0.7833333 82  660
6  0.7722628 83 1370
7  0.7878549 84 6340
8  0.7920635 85 3780
9  0.7988506 86 4350
10 0.8023256 87 4730
11 0.7948571 88 3500
12 0.8113333 89 6000
14 0.8017857 91  560
15 0.8620000 93  500


$condition_3
$root
[1] 49.46952

$predCI.root
  CI_2.5  CI_97.5 
49.10490 49.82108 

$b
[1] 0.8

$predCI.b
[1] 0.7947319 0.8051660

$terminated_early
[1] FALSE

$time
[1] 02m 1.50s

$iterations
[1] 100

$total.replications
[1] 32900

$tab
          y  x  reps
2 0.7763158 46   380
3 0.7632653 47   490
4 0.7768638 48 11670
5 0.7922311 49  5020
6 0.8111801 50 12880
7 0.8136000 51  1250


$condition_4
$root
[1] 33.63727

$predCI.root
  CI_2.5  CI_97.5 
33.50299 33.77255 

$b
[1] 0.8

$predCI.b
[1] 0.7951342 0.8047786

$terminated_early
[1] TRUE

$time
[1] 02m 0.83s

$iterations
[1] 98

$total.replications
[1] 31900

$tab
          y  x  reps
1 0.7355072 32  6900
2 0.7712743 33  9260
3 0.8133072 34 10220
4 0.8430876 35  4340

plot(solved, 1)

plot(solved, 2)

plot(solved, 3)

# also can plot median history and estimate precision
plot(solved, 1, type = 'history')

plot(solved, 1, type = 'density')

Warning in density.default(x, weights = reps/sum(reps)): Selecting bandwidth
*not* using 'weights'