#--------------------------------------------------------------------------- # Sample size for partial and full replicate design and scaled ABE # via simulated (empirical) power # estimation method via intra-subject contrasts & BE decision via # linearized reference scale BE criterion, no cap # # Author: dlabes #--------------------------------------------------------------------------- sampleN.RSABE <- function(alpha=0.05, targetpower=0.8, theta0, theta1, theta2, CV, design=c("2x3x3", "2x2x4", "2x2x3"), regulator = c("FDA", "EMA"), nsims=1E5, nstart, imax=100, print=TRUE, details=TRUE, setseed=TRUE) { if (missing(theta1) & missing(theta2)) theta1 <- 0.8 if (missing(theta2)) theta2=1/theta1 if (missing(theta1)) theta1=1/theta2 # according to the two Laszlos' paper 0.9 (alas 1.10) should be considered if (missing(theta0)) theta0 <- 0.90 # theta0 in range if ( (theta0 < theta1) | abs(theta0-theta1) <1e-8 | (theta0 > theta2) | abs(theta0-theta2) <1e-8 ){ stop("True ratio ", theta0," not within margins ", theta1," ... ", theta2,"!", call.=FALSE) } if (missing(CV)) stop("CV(s) must be given!", call.=FALSE) # regulator here only FDA, EMA # other regulatory bodies ("HC", "ANVISA") use all the EMA regulatory constant if (missing(regulator)) regulator <- "FDA" reg <- reg_check(regulator, choices=c("FDA", "EMA")) CVswitch <- reg$CVswitch r_const <- reg$r_const pe_constr <- reg$pe_constr # CVcap doesn't apply to the FDA recommended method # but in Munoz et al. method= Howe-EMA CVcap <- reg$CVcap # 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 <- log(1.0 + CVwT^2) s2wR <- log(1.0 + CVwR^2) # 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") { seqs <- 3 bk <- 1.5 # needed for n0? # in case of the FDA we are using the 'robust' df's # due to the fact that the described analysis in the # progesterone guidance is based in the intrasubject contrasts # T-R and R-R with df=n-seqs dfe <- parse(text="n-3", srcfile=NULL) dfRRe <- parse(text="n-3", srcfile=NULL) # expectation of mse of the ANOVA of intra-subject contrasts #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 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 <- (s2D + (s2wT + s2wR)/2) } if (design=="2x2x3") { 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){ designs <- c("2x2x4", "2x2x3", "2x3x3") type <- c("4 period full replicate", "3 period full replicate", "partial replicate") # clear words cat("\n++++++++ Reference scaled ABE crit. +++++++++\n") cat(" Sample size estimation\n") cat("---------------------------------------------\n") cat("Study design: ") cat(paste0(design, " (", type[match(design, designs)], ")\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 constraints =",theta1,"...", theta2,"\n") if (details | reg$name=="USER") { reg$CVcap <- NULL # CVcap doesn't apply here print(reg) cat("\n") } else { cat("Regulatory settings:", reg$name,"\n") } } # ----------------------------------------------------------------- # nstart from sampleN0 with widened limits # doesn't fit really good if theta0>=1.2 or <=0.85! ways out? ltheta1 <- -r_const*sqrt(s2wR) ltheta2 <- -ltheta1 # this does not function in case of CVwR=0.3 for the original code # calculating s2wR and back-calculating CVwR from that # numerical problem? if (CVwR <= CVswitch){ ltheta1 <- log(theta1) ltheta2 <- log(theta2) } 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 (regulator=="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(n100 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 <- .power.RSABE(mlog, sdm, C3, Emse, df, s2wR, dfRR, nsims, CVswitch, r_const, pe_constr, CVcap, ln_lBEL=log(theta1), ln_uBEL=log(theta2), alpha=alpha) pwr <- as.numeric(p["BE"]); # do not print first step down 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 one up again. is done in the next loop } while (pwrimax) break } nlast <- n if (up & pwrtargetpower & 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") # debug print # cat("iter=",iter+1,"\n") if (print | details) return(invisible(res)) else return(res) } # end function