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

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

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,911 posts in 4,806 threads, 1,635 registered users;
35 visitors (0 registered, 35 guests [including 8 identified bots]).
Forum time: 17:00 CET (Europe/Vienna)

The history of statistics is like a telephone directory:
the plot is boring, full of numbers and the cast is endless.    Stephen Senn

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