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

posted by d_labes  – Berlin, Germany, 2009-11-23 14:57 (4220 d 01:18 ago) – Posting: # 4378
Views: 9,607

(edited by d_labes on 2009-11-23 15:48)

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:

Activity
 Admin contact
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
BEBAC Ing. Helmut Schütz
HTML5