Exact power and sample size, Part IV [🇷 for BE/BA]

Dear All!

Ok, ok, I will quit with this "Eierlegende Wollmilchsau".
 # -------------------------------------------------------------------------- # Sample size for a desired power: see known.designs() # for covered experimental designs sampleN.TOST <- function(alpha=0.05, targetpower=0.8, logscale=TRUE,                          ltheta1=0.8, ltheta2, ldiff=0.95, CV,                          design="2x2", exact=TRUE, print=TRUE, details=FALSE) {   #number of the design   d.no <- .design.no(design)   if (is.na(d.no)) stop("Err: design not known")   # design charcteristics   d.name <- .design.name(d.no)   # get the df for the design as an unevaluated expression   dfe    <- .design.df(d.no,"Ns")   # get stepsize for search   steps  <- .design.step(d.no)   # get design constant   bk     <- .design.bk(d.no)   if (missing(CV)) stop("Err: CV must be given!")   # print the configuration:   if (print) {     cat("\n++++++++ Equivalence test - TOST ++++++++\n")     cat("         Sample size estimatio\n\n")     cat("-----------------------------------------\n")     cat("Study design: ",d.name,"\n")     if (details) {       cat("Design characteristics:\n")       cat("df = ",gsub("Ns","n",gsub(" ","",as.character(dfe))),           ", design const = ",bk,", step = ",steps,"\n\n", sep="")     }        }   # handle the log transformation   if (logscale) {     if (missing(ltheta2)) ltheta2=1/ltheta1     if ( (ldiff<=ltheta1) | (ldiff>=ltheta2) ) {       cat("Err: ratio not between margins!\n")       return(NaN)     }     theta1 <- log(ltheta1)     theta2 <- log(ltheta2)     diffm  <- log(ldiff)     se     <- sqrt(log(1.+CV^2))     if (print) cat("log-transformed data (multiplicative model)\n\n")   } else {     if (missing(ltheta2)) ltheta2=-ltheta1     if ( (ldiff<=ltheta1) | (ldiff>=ltheta2) ) {       cat("Err: diff not between margins!\n")       return(NaN)     }     theta1 <- ltheta1     theta2 <- ltheta2     diffm  <- ldiff     se     <- CV     if (print) cat("untransformed data (additive model)\n\n")   }   if (print) {     cat("alpha = ",alpha,", target power = ", targetpower,"\n", sep="")     cat("BE margins        =",ltheta1,"...", ltheta2,"\n")     if (logscale)       cat("Null (true) ratio = ",ldiff,",  CV = ",CV,"\n", sep="")     else       cat("Null (true) diff. = ",ldiff,",  CV = ",CV,"\n", sep="")   }   #start value from large sample approx. (hidden func.)   Ns  <- .sampleN0(alpha, targetpower, theta1, theta2, diffm, se, steps, bk)   df  <- eval(dfe)   if (exact) pow <- .power.TOST(alpha, theta1, theta2, diffm, se, n=Ns, df, bk)     else  pow <- .approx.power.TOST(alpha, theta1, theta2, diffm, se, n=Ns, df, bk)   if (details) {     cat("\nSample size search\n")     if (d.no == 0) cat("(n is sample size per group)\n")     cat(" n     power\n")     cat( Ns," ", formatC(pow, digits=6, format="f"),"\n")   }   # --- loop until power >= target power   iter <- 0   # iter>50 is emergency brake   # this is eventually not necessary, depends on quality of sampleN0   # in experimentation I have seen max of six steps   if (pow<targetpower) {     while (pow<targetpower) {       Ns   <- Ns+steps       iter <- iter+1       df   <- eval(dfe)       if (exact)         pow  <- .power.TOST(alpha, theta1, theta2, diffm, se, n=Ns, df, bk)       else         pow  <- .approx.power.TOST(alpha, theta1, theta2, diffm, se, n=Ns, df, bk)       if (details) cat( Ns," ", formatC(pow, digits=6, format="f"),"\n")       if (iter>50) break     }   } else {     while (pow>targetpower) {       if (Ns<=4) break     # min. sample size       Ns   <- Ns-steps     # step down if start power is to high       iter <- iter+1       df   <- eval(dfe)       if (exact)         pow  <- .power.TOST(alpha, theta1, theta2, diffm, se, n=Ns, df, bk)       else         pow  <- .approx.power.TOST(alpha, theta1, theta2, diffm, se, n=Ns, df, bk)       if (details) cat( Ns," ", formatC(pow, digits=6),"\n")       if (iter>50) break      }     if ( pow<targetpower ) {       Ns   <- Ns+steps     #step up once to have n with pow>=target       df   <- eval(dfe)       if (exact)         pow  <- .power.TOST(alpha, theta1, theta2, diffm, se, n=Ns, df, bk)       else         pow  <- .approx.power.TOST(alpha, theta1, theta2, diffm, se, n=Ns, df, bk)     }   }   if (print && !details) {     cat("\nSample size\n")     if (d.no == 0) cat("(n is sample size per group)\n")     cat(" n     power\n")     cat( Ns," ", formatC(pow, digits=6, format="f"),"\n")   }   if (details && print) {     if (exact) cat("\nExact power calculation with\nOwen's Q functions.\n")       else cat("\nApproximate power calculation with\nnon-central t-distribution.\n")   }   # always print if approx.   if (print && !exact)     cat("\nApproximate power calculation with\nnon-central t-distribution.\n")   if (print) cat("\n")   res <- data.frame(design=design, alpha=alpha, CV=CV, theta0=ldiff,                     theta1=ltheta1, theta2=ltheta2, n=Ns, power=pow,                     targetpower=targetpower)   names(res) <-c("Design","alpha","CV","theta0","theta1","theta2",                  "Sample size", "Achieved power", "Target power")   if (print)     return(invisible(res))   else     return(res) } 
The function .approx.power.TOST() not shown here is the code using non-central t-approximation for the power.

The whole code is released under the copy left (aka GPL).
Eventually there is someone out there to tying together all this snippets into a package, for christmas .

Regards,

Detlew