R/samplesize_scABEL2.R

In PowerTOST: Power and Sample Size for (Bio)Equivalence Studies

#-----------------------------------------------------------------------------
# Sample size for partial and full replicate designs and scaled ABE 
# via simulated (empirical) power
# estimation method intra-subject contrasts and BE decision via ABEL 
# (Average bioequivalence with expanding limits)
# 
# Author: dlabes
#-----------------------------------------------------------------------------

sampleN.scABEL2 <- function(alpha=0.05, targetpower=0.8, theta0, theta1, 
                            theta2, CV, design=c("2x3x3", "2x2x4", "2x2x3"), 
                            regulator, nsims=1E5, nstart, imax=100, print=TRUE, 
                            details=TRUE, setseed=TRUE)
{
  if (missing(theta1) & missing(theta2)) theta1 <- 0.8
  if (missing(theta1)) theta1=1/theta2
  if (missing(theta2)) theta2=1/theta1
  
  # the two Laszlo's recommend theta0=0.9 for HVD's
  if (missing(theta0)) theta0 <- 0.9
  if ( (theta0<=theta1) | (theta0>=theta2) ) {
    stop("True ratio ",theta0," not between margins ",theta1," / ",
         theta2,"!", call.=FALSE)
  }
  if (missing(CV)) stop("CV must be given!", call.=FALSE)
  
  #if (!print) details <- FALSE # do not print anything
  
  CVwT <- CV[1]
  if (length(CV)==2) CVwR <- CV[2] else CVwR <- CVwT
  s2wT <- log(1.0 + CVwT^2)
  s2wR <- log(1.0 + CVwR^2)
  
  if(missing(regulator)) regulator <- "EMA"
  # check regulator and get 
  # constants acc. to regulatory bodies (function in scABEL.R)
  reg <- reg_check(regulator)
  CVcap     <- reg$CVcap
  CVswitch  <- reg$CVswitch
  r_const   <- reg$r_const
  pe_constr <- reg$pe_constr
  
  # check design
  design <- match.arg(design)
  # we are treating only balanced designs
  # thus we use here bk - design constant for ntotal
  # expressions for the df's
  if (design=="2x3x3") {
    desi <- "2x3x3 (partial replicate)"
    seqs <- 3
    bk   <- 1.5    # needed for n0
    # in case of the ISC we are using the 'robust' df's
    # T-R and R-R with df=n-seqs
    dfe   <- parse(text="n-3", srcfile=NULL)
    dfRRe <- parse(text="n-3", srcfile=NULL)
    # according to McNally et al., verified via simulations:
    Emse  <- s2wT + s2wR/2
  }
  if (design=="2x2x4") {
    desi  <- "2x2x4 (full replicate)"
    seqs  <- 2
    bk    <- 1.0    # needed for n0
    dfe   <- parse(text="n-2", srcfile=NULL)
    dfRRe <- parse(text="n-2", srcfile=NULL)
    # expectation of mse of the ANOVA of intra-subject contrasts
    Emse  <- (s2wT + s2wR)/2 
  }
  if (design=="2x2x3") {
    desi <- "2x2x3 (TRT|RTR)"
    seqs <- 2
    bk   <- 1.5    # needed for n0?
    dfe   <- parse(text="n-2", srcfile=NULL)
    dfRRe <- parse(text="n/2-1", srcfile=NULL) # for balanced designs
    # expectation of mse of the ANOVA of intra-subject contrasts
    Emse  <- 1.5*(s2wT + s2wR)/2               # for balanced designs
  }
  
  mlog <- log(theta0)
  
  if (print){
    cat("\n+++++++++++ scaled (widened) ABEL +++++++++++\n")
    cat("            Sample size estimation\n")
    cat("(simulations based on intra-subject contrasts)\n")
    cat("----------------------------------------------\n")
    cat("Study design: ",desi,"\n")
    cat("log-transformed data (multiplicative model)\n")
    cat(nsims,"studies for each step simulated.\n\n")
    cat("alpha  = ", alpha,", target power = ", targetpower,"\n", sep="")
    cat("CVw(T) = ", CVwT,"; CVw(R) = ", CVwR,"\n", sep="")
    cat("True ratio = ",theta0,"\n", sep="")
    cat("ABE limits / PE constraint =", theta1,"...", theta2,"\n")
    if (details | reg$name=="USER") { 
      print(reg)
      cat("\n")
    } else {
      cat("Regulatory settings:", reg$name,"\n")
    }
  }
  
  # -----------------------------------------------------------------
  # nstart? from sampleN0 with widened limits
  # does'nt fit really good if theta0>=1.2! ways out? see sampleN0.2
  ltheta1 <- -sqrt(s2wR)*r_const
  ltheta2 <- -ltheta1
  if (CVwR <= CVswitch){
    ltheta1 <- log(theta1)
    ltheta2 <- log(theta2)
  }
  if (CVwR > CVcap){
    ltheta1 <- -sqrt(log(1.0 + CVcap^2))*r_const
    ltheta2 <- -ltheta1
  }
  if (missing(nstart)){
    # try to use the empirical alpha for start sample size, some sort of
    al <- alpha
    if (reg$name=="FDA") {
      if(Emse/bk <= CV2mse(0.30001) & Emse/bk >= CV2mse(0.2975)) al=0.12
      if(Emse/bk > CV2mse(0.30001)) al <- 0.035   
    }
    if (reg$name=="EMA") {
      #does not fit!
      #if(Emse/bk <= CV2mse(0.321) & Emse/bk >= CV2mse(0.28)) al=0.065
    }
    # debug print
    # cat(al,"\n")
    # we use bk=1 here since our formula is sem=sqrt(Emse/n)
    n01 <- .sampleN0(alpha=al, targetpower, ltheta1, ltheta2, diffm=mlog, 
                     se=sqrt(Emse), steps=seqs, bk=1, diffmthreshold=0.01)
    # empirical correction in the vicinity of CV=0.3, for ratios 
    # outside 0.86 ... 1/0.86
    #     if(Emse/bk <= CV2mse(0.305) & Emse/bk >= CV2mse(0.295) & abs(mlog)>log(1/0.865)) {
    #       if (reg$name=="EMA") n01 <- 0.9*n01 else  n01 <- 0.8*n01
    #       n01 <- seqs*trunc(n01/seqs)
    #     }  
    # start from PE constraint sample size
    n02 <- .sampleN0.2(targetpower, ltheta2=log(theta2), diffm=mlog, 
                       se=sqrt(Emse), steps=seqs, bk=1)
    # debug print
    # cat(n01,n02,"\n")
    n <- max(c(n01,n02))+seqs
  } else n <- seqs*round(nstart/seqs)
  nmin <- 6 # fits 2x3x3 and 2x2x4
  if(n<nmin) n <- nmin
  # we are simulating for balanced designs
  C3 <- 1/n
  # sd of the sample mean T-R (point estimator)
  sdm  <- sqrt(Emse*C3)
  df   <- eval(dfe)
  dfRR <- eval(dfRRe)
  
  if(setseed) set.seed(123456)
  
  p <- .pwr.ABEL.ISC(mlog, sdm, C3, Emse, df, s2wR, dfRR, nsims, 
                     CVswitch=CVswitch, r_const=r_const, CVcap=CVcap, 
                     pe_constr=pe_constr,
                     ln_lBEL=log(theta1),ln_uBEL=log(theta2), alpha=alpha)
  pwr <- as.numeric(p["BE"]);
  pd <- max(4,round(log10(nsims),0)-1)  # digits for power
  if (details) {
    cat("\nSample size search\n")
    cat(" n     power\n")
    # do not print first too high
    # this is for cases with only one step-down and than step up
    if (pwr<=targetpower) cat( n," ", formatC(pwr, digits=pd, format="f"),"\n")
  }
  iter <- 0
  # iter>100 is emergency brake
  # --- loop until power <= target power, step-down
  down <- FALSE
  up   <- FALSE
  while (pwr>targetpower) {
    down <- TRUE
    if (n<=nmin) { 
      if (details & iter==0) cat( n," ", formatC(pwr, digits=pd, format="f"),"\n")
      break
    }
    n  <- n-seqs     # step down if start power is to high
    iter <- iter + 1
    C3 <- 1/n
    # sd of the sample mean T-R (point estimator)
    sdm  <- sqrt(Emse*C3)
    df   <- eval(dfe)
    dfRR <- eval(dfRRe)
    if(setseed) set.seed(123456)
    p <- .pwr.ABEL.ISC(mlog, sdm, C3, Emse, df, s2wR, dfRR, nsims, 
                       CVswitch=CVswitch, r_const=r_const, CVcap=CVcap, 
                       pe_constr=pe_constr,
                       ln_lBEL=log(theta1),ln_uBEL=log(theta2), alpha=alpha)
    pwr <- as.numeric(p["BE"]);
    
    if (details) cat( n," ", formatC(pwr, digits=pd, format="f"),"\n")
    if (iter>imax) break  
    # loop results in n with power too low
    # must step up again one step. is done in the next loop
  }
  while (pwr<targetpower) {
    up   <- TRUE; down <- FALSE
    n    <- n+seqs   # step-up
    iter <- iter+1
    C3 <- 1/n
    # sd of the sample mean T-R (point estimator)
    sdm  <- sqrt(Emse*C3)
    df   <- eval(dfe)
    dfRR <- eval(dfRRe)
    if(setseed) set.seed(123456)
    p <- .pwr.ABEL.ISC(mlog, sdm, C3, Emse, df, s2wR, dfRR, nsims, 
                       CVswitch=CVswitch, r_const=r_const, CVcap=CVcap, 
                       pe_constr=pe_constr,
                       ln_lBEL=log(theta1),ln_uBEL=log(theta2), alpha=alpha)
    pwr <- as.numeric(p["BE"]);
    
    if (details) cat( n," ", formatC(pwr, digits=pd, format="f"),"\n")
    if (iter>imax) break 
  }

  nlast <- n
  if (up & pwr<targetpower) {
    n <- NA
    if (details) cat("Sample size search failed!\n")
  }
  if (down & pwr>targetpower & n != nmin) {
    n <- NA
    if (details) cat("Sample size search failed!\n")
  }
  
  if (print && !details) {
    cat("\nSample size\n")
    cat(" n     power\n")
    cat( n," ", formatC(pwr, digits=pd, format="f"),"\n")
    if (is.na(n)) cat("Sample size search failed!\n")
  }
  if (print) cat("\n")
  
  #return results as data.frame
  res <- data.frame(design=design, alpha=alpha, CVwT=CVwT, CVwR=CVwR,
                    theta0=theta0, theta1=theta1, theta2=theta2, n=n, power=pwr, 
                    targetpower=targetpower,nlast=nlast)
  names(res) <-c("Design", "alpha", "CVwT", "CVwR", "theta0", "theta1", "theta2",
                 "Sample size", "Achieved power", "Target power", "nlast")

  #cat("iter=",iter,"\n")
  
  if (print | details) return(invisible(res)) else return(res)
  
} # end function