Power to Detect “Mediation”






The title says “mediation” because it is incredibly difficult to find valid evidence of mediation. A ‘significant’ statistical test of mediation doesn’t mean a lot. Keep in mind that mediation is always a strict statement about causal effects unfolding over time.

Bullock, J. G., Green, D. P., & Ha, S. E. (2010). Yes, but what’s the mechanism? (don’t expect an easy answer). Journal of Personality and Social Psychology, 98(4), 550–558. http://doi.org/10.1037/a0018933

You’ll definitely want to check out docs for lavaan and simsem.

setwd("C:/Users/flournoy/Documents")

#install.packages('lavaan')
library(lavaan)
## This is lavaan 0.5-18
## lavaan is BETA software! Please report any bugs.
#install.packages('simsem')
library(simsem)
##  
## ###############################################################################################
## This is simsem 0.5-8
## simsem is BETA software! Please report any bugs.
## simsem was developed at the University of Kansas Center for Research Methods and Data Analysis.
## ###############################################################################################
#install.packages('dplyr')
library(dplyr)
## Warning: package 'dplyr' was built under R version 3.2.2
## 
## Attaching package: 'dplyr'
## 
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## 
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
#install.packages('ggplot2')
library(ggplot2)
# install.packages('snow')
# install.packages('cowplot')
library(cowplot)
## Warning: package 'cowplot' was built under R version 3.2.2
## 
## Attaching package: 'cowplot'
## 
## The following object is masked from 'package:ggplot2':
## 
##     ggsave
effectSizes <- c(.1,  .3, .5)

modelEffects <- expand.grid(effectSizes, effectSizes, effectSizes)
names(modelEffects) <- c('a', 'b', 'c')

models <- modelEffects %>%
    rowwise() %>%
    do({
        genModel <- paste0('# direct effect
Y ~ ', .$c, '*X
# mediator
M ~ ', .$a, '*X
Y ~ ', .$b, '*M
X ~~ 1*X
Y ~~ 1*Y
M ~~ 1*M
')
        testModel <-'# 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)
'
        data.frame(a=.$a, b=.$b, c=.$c, 
               gen=genModel, test=testModel, stringsAsFactors=F)
    })

head(models)
## Source: local data frame [6 x 5]
## 
##       a     b     c
##   (dbl) (dbl) (dbl)
## 1   0.1   0.1   0.1
## 2   0.3   0.1   0.1
## 3   0.5   0.1   0.1
## 4   0.1   0.3   0.1
## 5   0.3   0.3   0.1
## 6   0.5   0.3   0.1
## Variables not shown: gen (chr), test (chr)
allModelPowerSim <- models %>%
    rowwise() %>%
    do({
        manySims <- sim(NULL, model=.$test[1], n=seq(50, 1500, 50), generate=.$gen[1], 
            lavaanfun='sem', multicore=F) # Enable multicore on your personal computer
        data_frame(a=.$a, b=.$b, c=.$c, powersims=list(manySims))
    })
## Loading required package: parallel
## Progress: 1 / 30 
## Progress: 2 / 30 
## Progress: 3 / 30 
## Progress: 4 / 30 
## Progress: 5 / 30 
## Progress: 6 / 30 
...
.....
...... etc
# #saving the above
# saveRDS(allModelPowerSim, 'power_simulations.RData')

#loading the above
# allModelPowerSim <- readRDS('power_simulations.RData')

#?getPower
as.data.frame(aPowerTable.forab <- getPower(allModelPowerSim$powersims[[1]], 
                       nVal=seq(50, 1500, 25),
                       powerParam='ab'))
##    iv.N logical...j.
## 1    50  0.007163322
## 2    75  0.008304711
## 3   100  0.009626203
## 4   125  0.011155612
## 5   150  0.012924843
## 6   175  0.014970419
## 7   200  0.017334055
## 8   225  0.020063279
## 9   250  0.023212065
## 10  275  0.026841493
## 11  300  0.031020395
## 12  325  0.035825952
## 13  350  0.041344201
## 14  375  0.047670399
## 15  400  0.054909144
## 16  425  0.063174174
## 17  450  0.072587721
## 18  475  0.083279281
## 19  500  0.095383654
## 20  525  0.109038097
## 21  550  0.124378459
## 22  575  0.141534184
## 23  600  0.160622154
## 24  625  0.181739433
## 25  650  0.204955129
## 26  675  0.230301757
## 27  700  0.257766693
## 28  725  0.287284487
## 29  750  0.318730945
## 30  775  0.351919917
## 31  800  0.386603657
## 32  825  0.422477334
## 33  850  0.459187871
## 34  875  0.496346777
## 35  900  0.533546084
## 36  925  0.570376050
## 37  950  0.606443030
## 38  975  0.641385878
## 39 1000  0.674889467
## 40 1025  0.706694361
## 41 1050  0.736602177
## 42 1075  0.764476719
## 43 1100  0.790241409
## 44 1125  0.813873825
## 45 1150  0.835398288
## 46 1175  0.854877415
## 47 1200  0.872403450
## 48 1225  0.888089991
## 49 1250  0.902064541
## 50 1275  0.914462123
## 51 1300  0.925420057
## 52 1325  0.935073878
## 53 1350  0.943554304
## 54 1375  0.950985109
## 55 1400  0.957481764
## 56 1425  0.963150690
## 57 1450  0.968088972
## 58 1475  0.972384436
## 59 1500  0.976115964
powerData <- allModelPowerSim %>% rowwise() %>%
    do({
        aSimPower <- as.data.frame(getPower(.$powersims, 
                                            nVal=seq(50, 1500, 25),
                                            powerParam='ab'))
        data_frame(a=.$a, b=.$b, c=.$c, 
                   alab=paste0('a=',.$a), 
                   blab=paste0('b=',.$b), 
                   clab=paste0('c=',.$c), 
               N=aSimPower[,1],
               ab=aSimPower[,2])
    })

print(powerData)
## Source: local data frame [1,593 x 8]
## Groups: <by row>
## 
##        a     b     c  alab  blab  clab     N          ab
##    (dbl) (dbl) (dbl) (chr) (chr) (chr) (dbl)       (dbl)
## 1    0.1   0.1   0.1 a=0.1 b=0.1 c=0.1    50 0.007163322
## 2    0.1   0.1   0.1 a=0.1 b=0.1 c=0.1    75 0.008304711
## 3    0.1   0.1   0.1 a=0.1 b=0.1 c=0.1   100 0.009626203
## 4    0.1   0.1   0.1 a=0.1 b=0.1 c=0.1   125 0.011155612
## 5    0.1   0.1   0.1 a=0.1 b=0.1 c=0.1   150 0.012924843
## 6    0.1   0.1   0.1 a=0.1 b=0.1 c=0.1   175 0.014970419
## 7    0.1   0.1   0.1 a=0.1 b=0.1 c=0.1   200 0.017334055
## 8    0.1   0.1   0.1 a=0.1 b=0.1 c=0.1   225 0.020063279
## 9    0.1   0.1   0.1 a=0.1 b=0.1 c=0.1   250 0.023212065
## 10   0.1   0.1   0.1 a=0.1 b=0.1 c=0.1   275 0.026841493
## ..   ...   ...   ...   ...   ...   ...   ...         ...
ggplot(powerData, aes(x=N, y=ab))+
    geom_point(aes(color=clab),alpha=.5)+
  facet_grid(alab ~ blab, as.table=F)+
  geom_hline(y=.8,color='red',alpha=.5)

aPowerTable <- getPower(allModelPowerSim$powersims[[1]])
#?findPower
findPower(aPowerTable, 
         iv="N",
         power=.8)
##     c     a     b  Y~~Y  M~~M    ab total 
##   766   716   730   Inf   Inf  1111   766

You could try to make a data table of the power for finding the true mediated effect, ab, at each effect size level using ‘do’, below (uncomment, first, of course).

# allModelPowerSim %>% 
#   group_by(a, b, c) %>%
#   do({
#     aPowerTable <- getPower()
#     theNfor80percPower <- findPower()
#     data.frame()
#   })



Post a comment

You may use the following HTML:
<a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>