R/medserial.R In pwr2ppl: Power Analyses for Common Designs (Power to the People)

Documented in medserial

```#'Compute Power for Serial Mediation Effects
#'Requires correlations between all variables as sample size.
#'This approach calculates power for the serial mediation using
#'joint significance (recommended)
#'@param rxy Correlation between DV (y) and predictor (x)
#'@param rxm1 Correlation between predictor (x) and first mediator (m1)
#'@param rxm2 Correlation between predictor (x) and second mediator (m2)
#'@param rym1 Correlation between DV (y) and first mediator (m1)
#'@param rym2 Correlation between DV (y) and second mediator (m2)
#'@param rm1m2 Correlation first mediator (m1) and second mediator (m2)
#'@param alpha Type I error (default is .05)
#'@param rep number of repetitions (1000 is default)
#'@param n sample size
#'@examples
#'\donttest{medserial(rxm1=.3, rxm2=.3, rxy=-.35,
#'rym1=-.5,rym2=-.5, rm1m2=.7,n=150)}
#'@return Power for Serial Mediated (Indirect) Effects
#'@export
#'
#'

medserial<-function(rxm1,rxm2,rxy,rm1m2,rym1,rym2,n,alpha=.05, rep=1000)
{
V1<-NA;V2<-NA;V3<-NA;V4<-NA
pop <- MASS::mvrnorm(100000, mu = c(0,0,0,0),
Sigma = matrix(c(1.0,rxm1,rxm2, rxy,
rxm1,1.0,rm1m2, rym1,
rxm2,rm1m2,1.0,rym2,
rxy, rym1, rym2,1.0), ncol = 4),
empirical = TRUE)
pop<-as.data.frame(pop)
pop<-dplyr::rename(pop, x = V1, m1 = V2, m2 = V3, y = V4)
set.seed(1234)
nruns = rep
a = numeric(nruns)
b = numeric(nruns)
pa1<-NA
pa2<-NA
pb1<-NA
pb2<-NA
pd<-NA
for (i in 1:nruns)
{
samp <- pop[ sample(nrow(pop), n), ]
pathsa1<-summary(lm(m1~x,data=samp))
pa1[i]<-pathsa1\$coefficients[2,4]
pathsa2d<-summary(lm(m2~x+m1,data=samp))
pa2[i]<-pathsa2d\$coefficients[2,4]
pd[i]<-pathsa2d\$coefficients[3,4]
pathsb<-summary(lm(y~x+m1+m2,data=samp))
pb1[i]<-pathsb\$coefficients[3,4]
pb2[i]<-pathsb\$coefficients[4,4]

power = data.frame(pa1=pa1,pa2=pa2, pb1=pb1,pb2=pb2,pd=pd)
power\$jsa1["JointSiga1"]<-NA
power\$jsa1[pa1 < alpha] <- 1
power\$jsa1[pa1 >= alpha] <- 0
power\$jsa2["JointSiga2"]<-NA
power\$jsa2[pa2 < alpha] <- 1
power\$jsa2[pa2 >= alpha] <- 0
power\$jsb1["JointSigb1"]<-NA
power\$jsb1[pb1 < alpha] <- 1
power\$jsb1[pb1 >= alpha] <- 0
power\$jsb2["JointSigb2"]<-NA
power\$jsb2[pb2 < alpha] <- 1
power\$jsb2[pb2 >= alpha] <- 0
power\$jsd["JointSigd"]<-NA
power\$jsd[pd < alpha] <- 1
power\$jsd[pd >= alpha] <- 0
power\$jointm1<-power\$jsa1*power\$jsb1
power\$jointm2<-power\$jsa2*power\$jsb2
power\$jointm12<-power\$jsa1*power\$jsd*power\$jsb2
JSm1<-mean(power\$jointm1)
JSm2<-mean(power\$jointm2)
JSm12<-mean(power\$jointm12)}
message("Power for n = ", n,", mediator 1", " is ", JSm1)
message("Power for n = ", n,", mediator 2", " is ", JSm2)
message("Power for n = ", n,", serial mediation", " is ", JSm12)
result <- data.frame(matrix(ncol = 4))
colnames(result) <- c( "n","Power Mediator 1 (x-ma-mb)", "Power Mediator 2 (ma-mb-y)","Power Double (x-ma-mb-y)")
result[, 1]<-n
result[, 2]<-JSm1
result[, 3]<-JSm2
result[, 4]<-JSm12

output<-na.omit(result)
rownames(output)<- c()
invisible(output)

}
```

Try the pwr2ppl package in your browser

Any scripts or data that you put into this service are public.

pwr2ppl documentation built on Sept. 6, 2022, 5:06 p.m.