# R/powerMediation.R In powerMediation: Power/Sample Size Calculation for Mediation Analysis

```# Models (Sobel, MacKinnon et al.):
# mi=theta0+theta1*xi+ei, ei~N(0, sigma.e^2)
# yi=gamma+lambda*mi+lambda2*xi+epsiloni, epsiloni~N(0, sigma.epsilone^2)
# H0: theta_1=0 and lambda=0 versus theta_1 neq 0 or lambda neq 0
#
# Model (Vittinghoff and McCulloch, 2009)
# (1) y=a0+a1*x+epsilon
# (2) y=b0+b1*x+b2*m+e, e~N(0, sigma.e^2)
# H0: b2=0 vs Ha: b2 neq 0
#
#
# p-value and CI for mediation effect based on Sobel's test
testMediation.Sobel <- function(theta.1.hat, lambda.hat,
sigma.theta1, sigma.lambda, alpha=0.05)
{
part1<-(theta.1.hat*sigma.lambda)^2
part2<-(lambda.hat*sigma.theta1)^2

sigma.hat<-sqrt(part1+part2)
Z.obs<-theta.1.hat*lambda.hat/sigma.hat

pval<-2*(1-pnorm(abs(Z.obs)))

za2<-qnorm(1-alpha/2)

t1<-theta.1.hat*lambda.hat
t2<-za2*sigma.hat

CI.low<-t1-t2
CI.upp<-t1+t2

res<-list(pval=pval, CI.low=CI.low, CI.upp=CI.upp)
return(res)
}

# p-value and CI for mediation effect based on MacKinnon's test
testMediation.MacKinnon <- function(n, theta.1.hat, lambda.hat,
sigma.theta1.hat, sigma.lambda.hat, alpha=0.05)
{
numer<-theta.1.hat*lambda.hat
denom<-sigma.theta1.hat*sigma.lambda.hat
C.obs<-numer/denom

Z1<-theta.1.hat/sigma.theta1.hat
Z2<-lambda.hat/sigma.lambda.hat
tmp<-CDFprod(z=abs(C.obs), xmux=0, xmuy=0, rho=0, verbose=FALSE)\$prob
pval<-2*(1-tmp)

c1<-qprod(p=alpha/2, xmux=Z1, xmuy=Z2, rho=0, verbose=FALSE)\$z-C.obs
c2<-qprod(p=1-alpha/2, xmux=Z1, xmuy=Z2, rho=0, verbose=FALSE)\$z-C.obs

CI.low<-numer-c2*denom
CI.upp<-numer-c1*denom

res<-list(pval=pval, CI.low=CI.low, CI.upp=CI.upp)
return(res)
}

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

# based on MacKinnon et al.'s Test
# mi = theta0 + theta1*xi + ei, ei~N(0, sigma.e^2)
# yi = gamma + lambda * mi + lambda2*xi+epsiloni, epsiloni~N(0, sigma.epsilon^2)
# rho.mx=cor(x, m), var(m)=sigma.m^2, var(x)=sigma.x^2
#
# H0: theta1*lambda = 0 vs Ha: theta1*lambda=theta.1a*lambda.a neq 0
powerMediation.MacKinnon <- function(n, theta.1a, lambda.a, sigma.x,
sigma.m, sigma.e, sigma.epsilon, alpha=0.05, verbose=TRUE)
{
rho2.mx = (theta.1a*sigma.x/sigma.m)^2
if(rho2.mx >= 1 || rho2.mx < 0)
{
stop("rho2.mx should be in the range [0, 1)")
}

numer <- n * theta.1a * lambda.a * sigma.x * sigma.m * sqrt(1-rho2.mx)
denom <- sigma.e * sigma.epsilon
delta <- numer / denom

# expection of Z1 and Z2
EZ1=theta.1a*sigma.x*sqrt(n)/sigma.e
EZ2=lambda.a*sigma.m*sqrt(n*(1-rho2.mx))/sigma.epsilon

alpha2<-alpha/2
ca2 <- qprod(p = 1 - alpha2, xmux = 0, xmuy = 0, rho = 0, verbose = FALSE)
ca2 <- ca2\$z

p1 <- CDFprod(z=ca2-delta, xmux = EZ1, xmuy = EZ2, rho = 0, verbose = FALSE)
p1 <- p1\$prob
p2 <- CDFprod(z=-ca2-delta, xmux = EZ1, xmuy = EZ2, rho = 0,
verbose = FALSE)
p2 <- p2\$prob

power<-1-p1+p2

res<-list(power=power, delta=delta)
if(verbose)
{
print(res\$power)
}
invisible(res)
}

# based on Sobel First-Order Test
# mi = theta0 + theta1*xi + ei, ei~N(0, sigma.e^2)
# yi = gamma + lambda * mi + epsiloni, epsiloni~N(0, sigma.epsilon^2)
# H0: theta1*lambda = 0 vs Ha: theta1*lambda=theta.1a*lambda.a neq 0
powerMediation.Sobel <- function(n, theta.1a, lambda.a, sigma.x,
sigma.m, sigma.epsilon, alpha=0.05, verbose=TRUE)
{
rho2.mx = (theta.1a*sigma.x/sigma.m)^2
if(rho2.mx >= 1 || rho2.mx < 0)
{
stop("rho2.mx should be in the range [0, 1)")
}

numer<-sqrt(n)*theta.1a*lambda.a
# sigma.e^2=sigma.m^2*(1-rho.mx^2)
denom<-sqrt(theta.1a^2*sigma.epsilon^2/(sigma.m^2*(1-rho2.mx))+lambda.a^2*sigma.m^2*(1-rho2.mx)/sigma.x^2)
delta<-numer/denom

alpha2<-alpha/2
za2<-qnorm(1-alpha2)
power<-1-pnorm(za2-delta)+pnorm(-za2-delta)

res<-list(power=power, delta=delta)
if(verbose)
{
print(res\$power)
}
invisible(res)
}

# cumulative distribution of the product of two normal random variables
# that bivariate normal distribution
# wrapper for Alan Miller's fortran code fnprod
# http://jblevins.org/mirror/amiller/
# pr(X*Y < z)
CDFprod <- function(z, xmux, xmuy, rho, verbose = TRUE)
{
if(rho > 1 || rho < -1)
{
stop("rho should be in the range [-1, 1]")
}

res <- .Fortran("fnprod",
as.double(xmux),
as.double(xmuy),
as.double(rho),
as.double(z),
ier = as.integer(0),
abserr = as.double(0),
last = as.integer(0),
PACKAGE = "powerMediation")

# prob = Pr(X*Y < z), where (X, Y) has bivariate normal distr.
ier = res\$ier,
abserr = res\$abserr,
last = res\$last)

if(verbose)
{
print(resList\$prob)
}
return(resList)
}

tmpCDFprod <- function(z, p, xmux, xmuy, rho)
{
tmp <- CDFprod(z, xmux, xmuy, rho, verbose = FALSE)
res <- tmp\$prob - p
return(res)
}

# quantile of the product of two normal random variables
# that bivariate normal distribution.
# find z such that pr(X*Y < z)=p
qprod <- function(p,  xmux, xmuy, rho, z.lower=-1.0e+30, z.upper=1.0e+30,
verbose = TRUE)
{
if(p > 1 || p < 0)
{
stop("p should be in the range [0, 1]")
}
if(rho > 1 || rho < -1)
{
stop("rho should be in the range [-1, 1]")
}
if(z.lower >= z.upper)
{
stop("z.lower must be < z.upper")
}

res.uniroot <- uniroot(f = tmpCDFprod,
interval=c(z.lower, z.upper),
p = p, xmux = xmux, xmuy = xmuy, rho = rho)

z = res.uniroot\$root
if(verbose)
{
print(z)
}

res<-list(z = z, res.uniroot = res.uniroot)

invisible(res)
}

tmpSS.mediation.Sobel<-function(n, power, theta.1a, lambda.a,
sigma.x, sigma.m, sigma.epsilon,
alpha=0.05, verbose=FALSE)
{
tmppower<-powerMediation.Sobel(n=n, theta.1a=theta.1a,
lambda.a=lambda.a, sigma.x=sigma.x, sigma.m=sigma.m,
sigma.epsilon=sigma.epsilon,
alpha=alpha, verbose=verbose)
res<- tmppower\$power-power
return(res)
}

#####################################################
ssMediation.Sobel <- function(power, theta.1a, lambda.a,
sigma.x, sigma.m, sigma.epsilon,
n.lower=1, n.upper=1.0e+30,
alpha = 0.05, verbose=TRUE)
{
rho2.mx = (theta.1a*sigma.x/sigma.m)^2
if(n.lower< 1)
{
stop("n.lower must be >= 1")
}
if(n.lower >= n.upper)
{
stop("n.lower must be < n.upper")
}
if(power<0 || power > 1)
{
stop("power must be in the range [0, 1]")
}
if(rho2.mx >= 1 || rho2.mx < 0)
{
stop("rho2.mx should be in the range [0, 1)")
}

res.uniroot<-uniroot(f=tmpSS.mediation.Sobel,
interval=c(n.lower, n.upper),
power=power, theta.1a=theta.1a, lambda.a=lambda.a,
sigma.x=sigma.x, sigma.m=sigma.m,
sigma.epsilon=sigma.epsilon,
alpha=alpha, verbose=FALSE)

n.numeric<-res.uniroot\$root

res<-list(n=n.numeric, res.uniroot=res.uniroot)
if(verbose)
{ print(res\$n) }
invisible(res)
}

tmpSS.mediation.MacKinnon<-function(n, power, theta.1a, lambda.a,
sigma.x, sigma.m, sigma.e, sigma.epsilon,
alpha=0.05, verbose=FALSE)
{
tmppower<-powerMediation.MacKinnon(n=n, theta.1a=theta.1a,
lambda.a=lambda.a, sigma.x=sigma.x, sigma.m=sigma.m,
sigma.e=sigma.e, sigma.epsilon=sigma.epsilon,
alpha=alpha, verbose=verbose)
res<- tmppower\$power-power
return(res)
}

#####################################################
ssMediation.MacKinnon <- function(power, theta.1a, lambda.a,
sigma.x, sigma.m, sigma.e, sigma.epsilon, n.lower=1, n.upper=1e+30,
alpha = 0.05, verbose=TRUE)
{
rho2.mx = (theta.1a*sigma.x/sigma.m)^2
if(n.lower< 1)
{
stop("n.lower must be >= 1")
}
if(n.lower >= n.upper)
{
stop("n.lower must be < n.upper")
}
if(power<0 || power > 1)
{
stop("power must be in the range [0, 1]")
}
if(rho2.mx >= 1 || rho2.mx < 0)
{
stop("rho2.mx should be in the range [0, 1)")
}

res.uniroot<-uniroot(f=tmpSS.mediation.MacKinnon,
interval=c(n.lower, n.upper),
power=power, theta.1a=theta.1a, lambda.a=lambda.a,
sigma.x=sigma.x, sigma.m=sigma.m,
sigma.e=sigma.e, sigma.epsilon=sigma.epsilon,
alpha=alpha, verbose=FALSE)

n.numeric<-res.uniroot\$root

res<-list(n=n.numeric, res.uniroot=res.uniroot)
if(verbose)
{ print(res\$n) }
invisible(res)
}

#########################################################
# for simple linear regression y=theta + lambda * x + epsilon
#########################################################

power.SLR<-function(n, lambda.a, sigma.x, sigma.y, alpha=0.05, verbose=TRUE)
{
if(n <= 2)
{
stop("n must be > 2")
}
if(lambda.a > sigma.y/sigma.x || lambda.a < -sigma.y/sigma.x)
{
stop("lambda.a must be in the range [-sigma.y/sigma.x, sigma.y/sigma.x]")
}

s=sqrt(sigma.y^2-(lambda.a*sigma.x)^2)
delta = lambda.a*sigma.x*sqrt(n)/s

rho=lambda.a*sigma.x/sigma.y

alpha2<-alpha/2
n2<-n-2
# critical value of t
# i.e. upper 100*alpha/2% percentile
# Pr(t>t.cr)=alpha/2, t~t_{n-2}
t.cr<-qt(1-alpha2, df=n2)
part1<-1-pt(t.cr-delta, df=n2)
part2<-pt(-t.cr-delta, df=n2)

power<-part1+part2
res<-list(power=power, delta=delta, s=s, t.cr=t.cr, rho=rho)
if(verbose)
{ print(res\$power) }
invisible(res)
}

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

tmpPower.SLR<-function(n, power, lambda.a, sigma.x, sigma.y,
alpha=0.05, verbose=FALSE)
{
tmppower<-power.SLR(n=n, lambda.a=lambda.a, sigma.x=sigma.x,
sigma.y=sigma.y, alpha=alpha, verbose=verbose)
res<- tmppower\$power-power
return(res)
}

#####################################################
ss.SLR <- function(power, lambda.a, sigma.x, sigma.y,
n.lower = 2.01, n.upper = 1.0e+30, alpha=0.05, verbose=TRUE)
{
if(n.lower<= 2)
{
stop("n.lower must be > 2")
}
if(n.lower >= n.upper)
{
stop("n.lower must be < n.upper")
}
if(power<0 || power > 1)
{
stop("power must be in the range [0, 1]")
}
if(lambda.a > sigma.y/sigma.x || lambda.a < -sigma.y/sigma.x)
{
stop("lambda.a must be in the range [-sigma.y/sigma.x, sigma.y/sigma.x]")
}

res.uniroot<-uniroot(f=tmpPower.SLR,
interval=c(n.lower, n.upper),
power=power, lambda.a=lambda.a, sigma.x=sigma.x,
sigma.y=sigma.y, alpha=alpha, verbose=FALSE)

n.numeric<-res.uniroot\$root

res<-list(n=n.numeric, res.uniroot=res.uniroot)
if(verbose)
{ print(res\$n) }
invisible(res)
}

############################
tmp.minEffect.SLR<-function(lambda.a, n, power, sigma.x, sigma.y,
alpha=0.05, verbose=FALSE)
{

tmppower<-power.SLR(n=n, lambda.a=lambda.a, sigma.x=sigma.x,
sigma.y=sigma.y, alpha=alpha, verbose=verbose)
res<- tmppower\$power-power
return(res)
}

#####################################################
minEffect.SLR<-function(n, power, sigma.x, sigma.y,
alpha=0.05, verbose=TRUE)
{
if(n <= 2)
{
stop("n must be > 2")
}
if(power<0 || power > 1)
{
stop("power must be in the range [0, 1]")
}

res.uniroot<-uniroot(f=tmp.minEffect.SLR,
interval=c(0.0001, sigma.y/sigma.x-0.0001),
n=n, power=power, sigma.x=sigma.x,
sigma.y=sigma.y, alpha=alpha, verbose=FALSE)

lambda.a<-res.uniroot\$root

res<-list(lambda.a=lambda.a, res.uniroot=res.uniroot)
if(verbose)
{ print(res\$lambda.a) }
invisible(res)
}

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

#########################################################
# for simple linear regression y=theta + lambda * x + epsilon
#########################################################

power.SLR.rho<-function(n, rho2, alpha=0.05, verbose=TRUE)
{
if(n <= 2)
{
stop("n must be > 2")
}
if(rho2 > 1 || rho2 <= 0)
{
stop("rho2 must be in the range (0, 1]")
}

delta=sqrt(n)/sqrt(1/rho2-1)

alpha2<-alpha/2
n2<-n-2
# critical value of t
# i.e. upper 100*alpha/2% percentile
# Pr(t>t.cr)=alpha/2, t~t_{n-2}
t.cr<-qt(1-alpha2, df=n2)
part1<-1-pt(t.cr-delta, df=n2)
part2<-pt(-t.cr-delta, df=n2)

power<-part1+part2
res<-list(power=power, delta=delta)
if(verbose)
{ print(res\$power) }
invisible(res)
}

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

tmpPower.SLR.rho<-function(n, power, rho2,
alpha=0.05, verbose=FALSE)
{
tmppower<-power.SLR.rho(n=n, rho2=rho2, alpha=alpha, verbose=verbose)
res<- tmppower\$power-power
return(res)
}

#####################################################
ss.SLR.rho <- function(power, rho2,
n.lower = 2.01, n.upper = 1.0e+30, alpha=0.05, verbose=TRUE)
{
if(n.lower<= 2)
{
stop("n.lower must be > 2")
}
if(n.lower >= n.upper)
{
stop("n.lower must be < n.upper")
}
if(power<0 || power > 1)
{
stop("power must be in the range [0, 1]")
}
if(rho2 > 1 || rho2 <= 0)
{
stop("rho2 must be in the range (0, 1]")
}

res.uniroot<-uniroot(f=tmpPower.SLR.rho,
interval=c(n.lower, n.upper),
power=power, rho2=rho2, alpha=alpha, verbose=FALSE)

n.numeric<-res.uniroot\$root

res<-list(n=n.numeric, res.uniroot=res.uniroot)
if(verbose)
{ print(res\$n) }
invisible(res)
}
```

## Try the powerMediation package in your browser

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

powerMediation documentation built on March 24, 2021, 1:06 a.m.