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':
##
## ggsaveeffectSizes <- 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.976115964powerData <- 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 766You 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()
# })