R/power_VSMc_poisson.R

In powerMediation: Power/Sample Size Calculation for Mediation Analysis

# implmement sample size calculation methods proposed in
#  Vittinghoff E., Sen, S., and McCulloch, C.E. (2009). Sample size calculation for evaluating mediation. Statistics in Medicine. 28:541-557

# for poisson regression model
# b2 is the regression coefficient of mediator
# sigma.m is the standard deviation of mediator
# corr.xm is the correlation between predictor and mediator
# EY is the marginal mean of the outcome
ssMediation.VSMc.poisson<-function(power, b2,
  sigma.m, EY, corr.xm, n.lower=1, n.upper=1e+30,
  alpha = 0.05, verbose=TRUE)
{
  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(corr.xm < -1  || corr.xm > 1)
  {
    stop("corr.xm must be in the range [-1, 1]")
  }

  res.uniroot<-uniroot(f=tmpSS.mediation.VSMc.poisson,
      interval=c(n.lower, n.upper),
      power=power, b2=b2,
      sigma.m=sigma.m,
      EY=EY, corr.xm=corr.xm, 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)

}

# formula (9) in Vittinghoff, Sen, and McCulloch (2009). Statist. Med. 28:541-557
# alpha - Type I error rate for one-sided test (half of Type I error rate
#    for two-sided test)
ss.VSMc.poisson<-function(power, b2,
  sigma.m, EY, corr.xm, alpha = 0.025, verbose=TRUE)
{

  if(power<0 || power > 1)
  {
    stop("power must be in the range [0, 1]")
  }
  if(corr.xm < -1  || corr.xm > 1)
  {
    stop("corr.xm must be in the range [-1, 1]")
  }
  za<-qnorm(1-alpha)
  zg<-qnorm(power)
  numer<-(za+zg)^2
  denom<-(b2*sigma.m)^2*EY*(1-corr.xm^2)
  n.numeric<-numer/denom

  if(verbose)
  { print(n.numeric) }
  invisible(n.numeric)

}


powerMediation.VSMc.poisson<-function(n, b2, sigma.m, EY, corr.xm,
  alpha=0.05, verbose=TRUE)
{
  if(n < 1)
  {
    stop("n must be >= 1")
  }
  if(corr.xm < -1  || corr.xm > 1)
  {
    stop("corr.xm must be in the range [-1, 1]")
  }

  alpha2<-alpha/2

  za2<-qnorm(1-alpha2)
  delta<-b2*sigma.m*sqrt(n*(1-corr.xm^2))*sqrt(EY)

  power<-2-pnorm(za2-delta)-pnorm(za2+delta)

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


#########################
tmpSS.mediation.VSMc.poisson<-function(n, power, b2,
   sigma.m, EY, corr.xm,
   alpha=0.05, verbose=FALSE)
{
  tmppower<-powerMediation.VSMc.poisson(n=n, b2,
    sigma.m=sigma.m, EY=EY, corr.xm=corr.xm,
    alpha=alpha, verbose=verbose)
  res<- tmppower$power-power
  return(res)
}


############################
tmp.minEffect.VSMc.poisson<-function(b2, n, power, sigma.m, EY, 
   corr.xm, alpha=0.05, verbose=FALSE)
{

  tmppower<-powerMediation.VSMc.poisson(n=n, b2=b2, sigma.m=sigma.m, 
    EY=EY, corr.xm=corr.xm, alpha=alpha, verbose=verbose)
  res<- tmppower$power-power
  return(res)
}

#####################################################
minEffect.VSMc.poisson<-function(n, power, sigma.m, EY, corr.xm, 
  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.VSMc.poisson, 
      interval=c(0.0001, 1.0e+30),
      n=n, power=power, sigma.m=sigma.m, 
      EY=EY, corr.xm=corr.xm, alpha=alpha, verbose=FALSE)

  b2<-res.uniroot$root

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

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