R/power_RSABE2L_isc.R

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

#---------------------------------------------------------------------------
# Simulate partial or full replicate designs and scaled ABE power
# estimation method via intra-subject contrasts of T-R and R-R
# BE decision via "exact" method of the 2L, ABEL or hyslop's method
#  
# no cap if regulator = "FDA", else CVcap is taken into account if defined
# in regulator argument of class 'regSet'
# 
# Author: dlabes
#---------------------------------------------------------------------------

# degrees of freedom for the TR/RR  analysis: 
# Using the intrasubject contrasts T-R and R-R and analyze them  
# by sequence groups the df's are = n-seq.
# 2x3x3  dfRR = n-3
# 2x2x4  dfRR = n-2
# 2x2x3  dfRR = n/2 - 2

power.RSABE2L.isc <- function(alpha=0.05, theta1, theta2, theta0, CV, n,   
                              design=c("2x3x3", "2x2x4", "2x2x3"), regulator,
                              SABE_test="exact", k_est=TRUE,
                              nsims=1E5, details=FALSE, setseed=TRUE)
{
  if (missing(CV)) stop("CV must be given!")
  if (missing(n))  stop("Number of subjects n must be given!")

  if (missing(theta0)) theta0 <- 0.90
  if (missing(theta1) & missing(theta2)) theta1 <- 0.8
  if (missing(theta2)) theta2 <- 1/theta1
  if (missing(theta1)) theta1 <- 1/theta2
  
  ptm <- proc.time()
  
  # for later enhancement taking into account the 
  # subject-by-formulation interaction
  s2D  <- 0 
  CVwT <- CV[1]
  if (length(CV)==2) CVwR <- CV[2] else CVwR <- CVwT
  s2wT <- CV2mse(CVwT)
  s2wR <- CV2mse(CVwR)

  # regulator here only FDA, EMA
  # other regulatory bodies ("HC", "ANVISA") use all the EMA regulatory constant
  if (missing(regulator)) regulator <- "FDA"
  rc <- reg_check(regulator, choices=c("FDA", "EMA"))
  CVswitch  <- rc$CVswitch
  r_const   <- rc$r_const
  pe_constr <- rc$pe_constr
  # CVcap doesn't apply to the FDA recommended method
  # but in Munoz et al. method = Howe-EMA
  CVcap     <- rc$CVcap

  # check SABE_test
  SABE_test <- tolower(SABE_test)
  SABE_test <- match.arg(SABE_test, choices=c("exact", "abel", "hyslop", "fda"))
  
  # check design and give the design characteristics
  design <- match.arg(design)
  if (design=="2x3x3") {
    seqs  <- 3
    dfe   <- parse(text="n-3", srcfile=NULL)
    dfRRe <- parse(text="n-3", srcfile=NULL)
    #sd2  <- s2D + (s2wT + s2wR)/2 # used in v1.1-00 - v1.1-02
    # according to McNally et al.
    # verified via simulations:
    Emse  <- s2D + s2wT + s2wR/2
  }
  if (design=="2x2x4") {
    seqs  <- 2
    dfe   <- parse(text="n-2", srcfile=NULL)
    dfRRe <- parse(text="n-2", srcfile=NULL)
    # sd^2 of the differences T-R from their components
    Emse  <- (s2D + (s2wT + s2wR)/2) 
  }
  if (design=="2x2x3") {
    seqs  <- 2
    dfe   <- parse(text="n-2", srcfile=NULL)
    # next was pre-V1.2-08
#     dfRRe <- parse(text="n/2-2", srcfile=NULL) # for balanced designs
#     dfTTe <- parse(text="n/2-2", srcfile=NULL) # for balanced designs
    # correct should be (only 1 sequence for each, f.i. RR from RTR):
    dfRRe <- parse(text="n/2-1", srcfile=NULL) # for balanced designs
    dfTTe <- parse(text="n/2-1", srcfile=NULL) # for balanced designs
    # sd^2 of the differences T-R from their components
    Emse  <- 1.5*(s2wT + s2wR)/2               # for balanced design 
  }
  
  if (length(n)==1){
    # then we assume n=ntotal
    # for unbalanced designs we divide the ns by ourself
    # to have only small imbalance
    nv <- nvec(n=n, grps=seqs)
    if (nv[1]!=nv[length(nv)]){
      message("Unbalanced design. n(i)=", paste(nv, collapse="/"), " assumed.")
    } 
    C3 <- sum(1/nv)/seqs^2
    n  <- sum(nv)
  } else {
    # then we assume n = vector of n's in sequences
    # check length
    if (length(n)!=seqs) stop("n must be a vector of length=",seqs,"!")
    
    C3 <- sum(1/n)/seqs^2
    nv <- n
    n  <- sum(n)
  }
  
  df <- eval(dfe)
  if (design=="2x2x3"){
    dfTT <- nv[1]-1
    dfRR <- nv[2]-1
    # where does this next came from?
    Emse <- (dfRR*(s2wT + s2wR/2)+dfTT*(s2wT/2 + s2wR))/(dfRR+dfTT)
    # warning in case of unbalanced design and heteroscdasticity
    if (abs(s2wT - s2wR)>1e-5 & abs(dfRR-dfTT)>2){
      warning(paste("Heteroscedasticity in unbalanced 2x2x3 design may lead", 
              "to poor accuracy of power!"), call.=FALSE)
    }
  } else {
    dfRR <- eval(dfRRe)
  }
  #cat("dfRR=", dfRR," dfTT=",dfTT," E(s2I)=", Emse, "\n")
  # sd of the mean T-R (point estimator)
  sdm  <- sqrt(Emse*C3)
  mlog <- log(theta0)
  
  if(setseed) set.seed(123456)
  p <- .pwr.SABE.isc(mlog, sdm, C3, Emse, df, s2wT, s2wR, dfRR, k_est, nsims, 
                     CVswitch, r_const, pe_constr, CVcap, SABE_test=SABE_test,
                     ln_lBEL=log(theta1),ln_uBEL=log(theta2), alpha=alpha)
    
  if (details) {
    ptm <- summary(proc.time()-ptm)
    message(nsims," sims. Time elapsed (sec): ", 
            formatC(ptm["elapsed"], digits=2), "\n")
    # return the components
    names(p) <- c("p(BE)", "p(BE-RSABE)", "p(BE-pe)", "p(BE-ABE)")
    if (SABE_test=="abel") names(p) <- c("p(BE)", "p(BE-ABEL)", "p(BE-pe)", "p(BE-ABE)")
    if (!pe_constr) p <- p[-3] # without pe constraint
    p
  } else {
    # return only the 'power'
    as.numeric(p["BE"])
  }
}

# ----------------------------------------------------------------------------
# working horse of SABE, estimation via ISC
.pwr.SABE.isc <- function(mlog, sdm, C3, Emse, df, s2wT, s2wR, dfRR, k_est, 
                          nsims, CVswitch, r_const, pe_constr, CVcap, 
                          SABE_test="exact", ln_lBEL, ln_uBEL, alpha=0.05)
{
  tval     <- qt(1-alpha,df)
  chisqval <- qchisq(1-alpha, dfRR)
  r2const  <- r_const^2
  s2switch <- log(CVswitch^2+1) 
  s2cap    <- log(CVcap^2+1)
  if(is.finite(CVcap)) wABEL_cap <- r_const*CV2se(CVcap) else wABEL_cap <- Inf
  
  counts <- rep.int(0, times=4)
  names(counts) <- c("BE", "BEul", "BEpe", "BEabe")
  # to avoid memory problems for high number of sims we work in chunks
  chunks <- 1
  nsi    <- nsims
  if (nsims>1E7) {
    chunks <- ceiling(nsims/1E7)
    nsi    <- 1E7
  } 
  for (iter in 1:chunks) {
    # if chunks*1E7 >nsims correct nsi to given nsims
    if(iter==chunks) nsi <- nsims-(chunks-1)*nsi
    # simulate sample mean (pe) of T-R via its normal distribution
    pes  <- rnorm(nsi, mean=mlog, sd=sdm)
    # simulate sample sd2s via chi-square distri
    sd2s   <- Emse*C3*rchisq(nsi, df)/df
    # simulate sample value s2wRs via chi-square distri
    s2wRs  <- s2wR*rchisq(nsi, dfRR)/dfRR
    SEs <- sqrt(sd2s)
    # conventional (1-2*alpha) CI's for T-R
    hw  <- tval*SEs
    lCL <- pes - hw 
    uCL <- pes + hw
    # conventional ABE, only for comparative purposes
    BE_ABE <- ((ln_lBEL<=lCL) & (uCL<=ln_uBEL))
    
    if (SABE_test=="exact"){
      # step 1: compute k
      if(k_est){
        # eqn (12), but is valid only for the population values?
        # is this valid for ANOVA only?
        k <- SEs/sqrt(s2wRs)
        # try to empirical correct the alpha overshot
        #k <- 1.04*k
        # geometric mean of 'estimated' and constant value
        #k <- sqrt(k * sdm/sqrt(s2wR))
        # geometric mean
        #k <- exp(mean(log(k)))
      } else {
        # use constant k based on population values
        #k <- sqrt((Emse/s2wR)*C3)
        k <- sdm/sqrt(s2wR)
        # maybe replaced by median(k) or mean(k)
      }
      #browser()
      # Hedges correction factor
      Hf <- 1-3/(4*dfRR-1)
      # step 2: compute L/U using eqn. (26)
      # attention! in the 2016 paper the non-centrality parm is defined
      # different, also the effect size
      # see f.i. eqn (17a, 17b) of the 2016 paper
      #
      # df for non-central t-distri; Which one?
      # here dfRR equals df, except for TRT|RTR
        #Ltheta <- qt(p=1-alpha, df=dfRR, ncp=-(Hf/k)*r_const)
        # suppress warnings wrt to full precision in pnt
        Ltheta <- suppressWarnings(-qt(p=alpha, df=dfRR, ncp=(Hf/k)*r_const))
        # using non-central f distribution
        # doesn't give the same values!!!
        #Ltheta <- -sqrt(qf(p=alpha, df1=1, df2=dfRR, ncp=(r_const*Hf/k)^2))
        #browser()
        #Utheta <- qt(alpha, dfRR, +(Hf/k)*r_const)
      # 2016 paper  
        #Ltheta <- qt(1-alpha, dfRR, -r_const/k)
        #Utheta <- qt(alpha, dfRR, +r_const/k)
        #Ltheta <- -qt(alpha, dfRR, +r_const/k)
        Utheta <- -Ltheta # is this in all cases correct?
      # effect size
      es <- (pes/sqrt(s2wRs))/k
      # 2016 paper
      #es <- pes/sqrt(s2wRs)/k/Hf
      # RSABE ("exact") decision
      BE_RSABE  <- (Ltheta < es) & (es < Utheta)
      #browser()
    } else if (SABE_test=="hyslop" | SABE_test=="howe" | SABE_test=="fda") {
      # linearized SABE criterion + 95%CI
      # with -SEs^2 the 'unknown' x from the progesterone guidance
      # and this gives the same values as power.RSABE()
      Em <- pes^2 
      if (SABE_test=="fda") Em <- Em - SEs^2 # bias corr.  
      Es <- r2const*s2wRs
      Cm <- (abs(pes) + hw)^2
      Cs <- Es*dfRR/chisqval    
      SABEc95 <- Em - Es + sqrt((Cm-Em)^2 + (Cs-Es)^2)
      # save memory
      rm(SEs, hw, Em, Es, Cm, Cs)
      BE_RSABE <- SABEc95<=0
    } else {
      # ABEL method
      # 'widened' acceptance limits or conventional limits below s2switch
      wABEL  <- ifelse(s2wRs<=s2switch, ln_uBEL, r_const*sqrt(s2wRs))
      # cap on widening
      wABEL  <- ifelse(s2wRs>s2cap, wABEL_cap, wABEL)
      # scaled ABE (ABEL) decision
      BE_RSABE  <- (lCL >= -wABEL) & (uCL <= wABEL)
    }
    # mixed procedure
    BE <- ifelse(s2wRs>s2switch, BE_RSABE, BE_ABE)
    # use capped acceptance limits if CVwR > CVcap
    # already incorporated with "abel"
    if (is.finite(CVcap) & SABE_test!="abel"){
      s2cap <- CV2mse(CVcap)
      # calculate the capped widened acceptance limits in log domain
      uprABEL <- r_const*sqrt(s2cap)
      lwrABEL <- -uprABEL
      BE <- ifelse(s2wRs>=s2cap, ((lwrABEL<=lCL) & (uCL<=uprABEL)), BE)
    }
    
    # point est. constraint true?
    BEpe  <- ( pes>=ln_lBEL & pes<=ln_uBEL )

    counts["BEabe"] <- counts["BEabe"] + sum(BE_ABE)
    counts["BEpe"]  <- counts["BEpe"]  + sum(BEpe)
    counts["BEul"]  <- counts["BEul"]  + sum(BE)
    if(pe_constr) {
      counts["BE"]    <- counts["BE"]    + sum(BE & BEpe)
    } else {
      counts["BE"]    <- counts["BE"]    + sum(BE) # no pe constraint
    }
    
  } # end over chunks
  # return the pBEs
  counts/nsims
}

# ----------------------------------------------------------------------------
# working horse for RSABE (FDA) based on .pwr.SABE.isc
.power.RSABE <- function(mlog, sdm, C3, Emse, df, s2wR, dfRR, nsims, 
                         CVswitch, r_const, pe_constr, CVcap,
                         ln_lBEL, ln_uBEL, alpha=0.05)
{
  pwr <- .pwr.SABE.isc(mlog, sdm, C3, Emse, df, s2wT=NULL, s2wR, dfRR, k_est=NULL, 
                       nsims, CVswitch, r_const, pe_constr, CVcap, 
                       SABE_test="fda", ln_lBEL, ln_uBEL, alpha=alpha)
  pwr
}

# ----------------------------------------------------------------------------
# working horse for ABEL with ISC est. based on .pwr.SABE.isc
# k_est is only used with SABE_test="exact"
.pwr.ABEL.ISC <- function(mlog, sdm, C3, Emse, df, s2wR, dfRR, nsims, 
                            CVswitch, r_const, pe_constr, CVcap,
                            ln_lBEL, ln_uBEL, alpha=0.05)
{
  pwr <- .pwr.SABE.isc(mlog, sdm, C3, Emse, df, s2wT=NULL, s2wR, dfRR, k_est=NULL, 
                       nsims, CVswitch, r_const, pe_constr, CVcap, 
                       SABE_test="abel", ln_lBEL, ln_uBEL, alpha=alpha)
  pwr
}