Bioequivalence and Bioavailability Forum

Dear All!

In the era of 'Appelssstrudel' and Owen's Q-function as implemented in R (see here) R-users don't have the necessity to use approximations but should calculate always exact power and sample size!

Since other designs (parallel group, cross-over with more than 2 formulations, replicate designs) used occasionally in BE studies also rely on the two-one-sided t-tests procedure (TOST) here comes my R-super-code:

# -- The functions of normal-, t-distributions and integrate() --------------
require(stats)
# ---------------------------------------------------------------------------
# Owen's Q-function
OwensQ <- function (nu, t, delta, a, b)
{
   # Observation: for 'really' large df (nu>2000) and large delta/b the 
   # density function is zero over nearly all its range. Q is than 
   # returned sometimes as =0! See documentation (?integrate) for that.
   # example: OwensQ(3000,1.64,-10,0,300) = 1
   #          OwensQ(3000,1.64,-10,0,350) = 5.614476e-12 !
   # Idea: adapt upper and/or lower to account for that
   low <- a; up<- b
   if (nu>=1000){
     # try to shorten the range over that is integrated
     x <- seq(a, b, by=(b-a)/400)
     dens <- .Q.integrand(x, nu, t, delta)
     r <- data.frame(x=x, dens=dens)
     r <- r[r$dens>0,]      if (length(r)>0) {# if any >0
       up <- max(r$x)+(b-a)/600  # only the upper x+step is used
       rm(r,x,dens)
     }
   }
   Qintegral <- integrate(.Q.integrand, lower = low, upper = up, nu=nu, t=t,
                delta = delta, subdivisions = 10000,
                rel.tol = 1.e-10, abs.tol=1.e-10,
                stop.on.error = TRUE, keep.xy = FALSE, aux = NULL)
   Q <- Qintegral[[1]]

   return(Q) 
}
# --------------------------------------------------------------------------
# Integrand of the definit integral in Owen's Q. Used in the call of integrate()
# Not useful alone, I think ? Leading . hides this function    
.Q.integrand <- function(x, nu, t, delta)
{ #version without for - loop
   lnQconst <- -((nu/2.0)-1.0)*log(2.0) - lgamma(nu/2.)
   dens <- pnorm( t*x/sqrt(nu) - delta, mean = 0, sd = 1, log = FALSE) *
                 exp( (nu-1)*log(x) - 0.5*x^2 + lnQconst )

   return(dens)
}
# --------------------------------------------------------------------------
# 'raw' power function to be used in sampleN.TOST avoiding overhead
# bk is the so called design constant (see later)
.power.TOST <- function(alpha=0.05, theta1, theta2, diffm, se, n, df, bk=2)
{
    tval   <- qt(1 - alpha, df, lower.tail = TRUE, log.p = FALSE)
    delta1 <- (diffm-theta1)/(se*sqrt(bk/n))
    delta2 <- (diffm-theta2)/(se*sqrt(bk/n))
    R      <- (delta1-delta2)*sqrt(df)/(2.*tval)
   
    # to avoid numerical errors in OwensQ implementation
    if (df>10000) {
      # Joulious formula (57) or (67), normal approximation
      p1 <- pnorm( (abs(delta1)-tval), lower.tail = TRUE)
      p2 <- pnorm( (abs(delta2)-tval), lower.tail = TRUE)
      return (p1 + p2 -1)
    }
   
    p1 <- OwensQ(df,  tval, delta1, 0, R)
    p2 <- OwensQ(df, -tval, delta2, 0, R)

    return( p2-p1 )
}

Continuation follows.

Use it on your own risk! The author does not take any responsibility for damage of your computer, your career or of the world .