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

posted by d_labes  – Berlin, Germany, 2009-11-23 15:57 (5261 d 03:27 ago) – Posting: # 4378
Views: 10,623

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 :cool:.
—
Regards,

Detlew

Complete thread:

UA Flag
Activity
 Admin contact
22,988 posts in 4,825 threads, 1,661 registered users;
95 visitors (1 registered, 94 guests [including 5 identified bots]).
Forum time: 20:25 CEST (Europe/Vienna)

The only way to comprehend what mathematicians mean by Infinity
is to contemplate the extent of human stupidity.    Voltaire

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