## Exact power and sample size, Part IV [R 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

21,518 posts in 4,498 threads, 1,523 registered users;
online 10 (0 registered, 10 guests [including 6 identified bots]).
Forum time: Sunday 17:16 CEST (Europe/Vienna)

A refund for defective software might be nice,
except it would bankrupt the entire software industry
in the first year.    Andrew S. Tanenbaum

The Bioequivalence and Bioavailability Forum is hosted by
Ing. Helmut Schütz