R/powerMediation.R

Defines functions ss.SLR.rho tmpPower.SLR.rho power.SLR.rho minEffect.SLR tmp.minEffect.SLR ss.SLR tmpPower.SLR power.SLR ssMediation.MacKinnon tmpSS.mediation.MacKinnon ssMediation.Sobel tmpSS.mediation.Sobel qprod tmpCDFprod CDFprod powerMediation.Sobel powerMediation.MacKinnon testMediation.MacKinnon testMediation.Sobel

Documented in minEffect.SLR powerMediation.Sobel power.SLR power.SLR.rho ssMediation.Sobel ss.SLR ss.SLR.rho testMediation.Sobel

# 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),
         answer = as.double(0),
         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.
     resList <- list(prob = res$answer, 
         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 Feb. 9, 2018, 6:10 a.m.