R/sampsiz_n0.R

In Detlew/Power2Stage: Power and Sample-Size Distribution of 2-Stage Bioequivalence Studies

# ----- helper functions ---------------------------------------------------
# sample size search start values
#
# Sample size for a desired power, large sample approx.
# original 'large sample' sample size formula if diffmthreshold=0
# author D. Labes
# vectorizes properly if diffm is scalar and se a vector or both are vectors
# does not give the correct answer if one acceptance limit goes to Inf
.sampleN00 <- function(alpha=0.05, targetpower=0.8, ltheta1, ltheta2, diffm,
                       se, steps=2, bk=2, diffmthreshold=0.04)
{
  # handle unsymmetrical BE acceptance ranges
  dn <- pmin(abs(diffm-ltheta1), abs(diffm-ltheta2))
  # return Inf if diffm outside, dn=0 accomplishes that
  dn <- ifelse((diffm-ltheta1)<1.25e-5 | (ltheta2-diffm)<1.25e-5, 0, dn)

  z1 <- qnorm(1-alpha)
  # transform to limits symmetric around zero (if they are not)
  locc    <- (ltheta1+ltheta2)/2
  diffm   <- diffm - locc
  ltheta1 <- ltheta1 - locc
  ltheta2 <- -ltheta1
  # within +-diffmthreshold diffm will be regarded as 0
  # value diffmthreshold=0.04 corresponds roughly to log(0.96)
  # with lower values there are many steps in sample size around between 0.95 and 1.
  # no longer used. Zhangs formula instead
  # vectorized solution
  z2    <- ifelse(abs(diffm)>diffmthreshold, qnorm(targetpower),
                  qnorm(1-(1-targetpower)/2))
  diffm <- ifelse(abs(diffm)>diffmthreshold, diffm, 0)
  # 1-beta/2 seems only correct if ltheta1 = -ltheta2 ?

  n0 <- (bk/2)*((z1+z2)*(se*sqrt(2)/dn))^2

  # round up to next even >=
  # seems Golkowski has used simple round
  n0 <- steps*ceiling(n0/steps)

  return(n0)

}

# ----------------------------------------------------------------------------
# 'large sample' sample size with one additional step via t-distribution
# author D. Labes
# bk = design constant, see known.designs()
# does not vectorize if se is a vector
.sampleN0 <- function(alpha=0.05, targetpower=0.8, ltheta1, ltheta2, diffm,
                      se, steps=2, bk=2, diffmthreshold=0.04)
{
  # paranoia? swap
  if (ltheta2<ltheta1) {h <- ltheta2; ltheta2 <- ltheta1; ltheta1 <-h}

  z1 <- qnorm(1-alpha)
  # value diffmthreshold=0.04 corresponds roughly to log(0.96)
  # with lower values or 0 there are many steps around between 0.95 and 1
  # in sampleN.TOST, no longer used but Zhang's formula instead
  # vectorized solution
  z2    <- ifelse(abs(diffm)>diffmthreshold, qnorm(targetpower), qnorm(1-(1-targetpower)/2))
  diffm <- ifelse(abs(diffm)>diffmthreshold, diffm, 0)

  #dn <- ifelse(diffm<0, diffm-ltheta1, diffm-ltheta2)
  # handle unsymmetrical BE acceptance ranges, one may also be Inf, -Inf
  dn <- pmin(abs(diffm-ltheta1), abs(diffm-ltheta2))
  # return Inf if diffm outside, dn=0 accomplishes that
  dn <- ifelse((diffm-ltheta1)<1.25e-5 | (ltheta2-diffm)<1.25e-5, 0, dn)
  n0 <- (bk/2)*((z1+z2)*(se*sqrt(2)/dn))^2

  if (is.finite(n0) & n0>2){
    # make another step with t-distri
    n0 <- ceiling(n0)
    z1 <- qt(1-alpha, df=n0-2)
    z2    <- ifelse(abs(diffm)>diffmthreshold, qt(targetpower, df=n0-2),
                    qt(1-(1-targetpower)/2, df=n0-2))

    n0 <- (bk/2)*((z1+z2)*(se*sqrt(2)/dn))^2
  }
  # make an even multiple of step (=2 in case of 2x2 cross-over)
  n0 <- steps*trunc(n0/steps)

  # minimum sample size will be checked outside
  return(n0)
}

# ----------------------------------------------------------------------------
# Paul Zhang (2003)
# A Simple Formula for Sample Size Calculation in Equivalence Studies
# Journal of Biopharmaceutical Statistics, 13:3, 529-538
.sampleN0_2 <- function(alpha=0.05, targetpower=0.8, ltheta1, ltheta2, diffm,
                        se, steps=2, bk=2)
{
  # handle unsymmetric limits, Zhang's c0
  c0 <- 0.5*exp(-7.06*(ltheta1+ltheta2)/(ltheta1-ltheta2))
  # Zhang's formula, large sample
  beta <- 1-targetpower
  z1 <- qnorm(1-alpha)
  fz <- ifelse(diffm<0, c0*exp(-7.06*diffm/ltheta1), c0*exp(-7.06*diffm/ltheta2))
  # results in NaN if one of the acceptance range is Inf/-Inf
  fz[is.nan(fz)|fz>=1] <- 0
  z2 <- abs(qnorm((1-fz)*beta))

  #dn <- ifelse(diffm<0, diffm-ltheta1, diffm-ltheta2)
  dn <- pmin(abs(diffm-ltheta1), abs(diffm-ltheta2))
  # return Inf if diffm outside, dn=0 accomplishes that
  dn <- ifelse((diffm-ltheta1)<1.25e-5 | (ltheta2-diffm)<1.25e-5, 0, dn)
  n0 <- (bk/2)*((z1+z2)*(se*sqrt(2)/dn))^2

  # round up to next even >=
  n0 <- steps*ceiling(n0/steps)

  return(n0)
}

# -----------------------------------------------------------------------------
# variant of Zhangs 'large sample' formula with 'smooth' transition from beta/2
# to beta and a threshold
# vectorizes properly if diffm is scalar and se a vector or both are vectors
# needs finite ltheta1, ltheta2!
.sampleN0_3 <- function(alpha=0.05, targetpower=0.8, ltheta1, ltheta2, diffm,
                        se, steps=2, bk=2)
{
  # difference for denominator
  #dn <- ifelse(diffm<0, diffm-ltheta1, diffm-ltheta2)
  dn <- pmin(abs(diffm-ltheta1), abs(diffm-ltheta2))
  # return Inf if diffm outside, dn=0 accomplishes that
  dn <- ifelse((diffm-ltheta1)<1.25e-5 | (ltheta2-diffm)<1.25e-5, 0, dn)

  # transform to limits symmetric around zero (if they are not)
  locc    <- (ltheta1+ltheta2)/2
  diffm   <- diffm - locc
  ltheta1 <- ltheta1 - locc
  ltheta2 <- -ltheta1
  delta   <- abs((ltheta2-ltheta1)/2)

  z1   <- qnorm(1-alpha)
  beta <- 1-targetpower

  c  <- ifelse(is.finite(diffm), abs(diffm/delta), 1)
  # probability for second normal quantil
  # if c<0.2 is in general a good choice needs to be tested
  p2 <- ifelse(c<0.2, 1-(1-0.5*exp(-7.06*c))*beta, 1-beta)
  z2 <- qnorm(p2)

  n0 <- (bk/2)*((z1+z2)*(se*sqrt(2)/dn))^2
  # round up to next even >=
  n0 <- steps*ceiling(n0/steps)

  return(n0)

}