Another one - R code for power of TROST [🇷 for BE/BA]

posted by d_labes  – Berlin, Germany, 2009-10-09 15:50 (5730 d 22:50 ago) – Posting: # 4335
Views: 34,387

Dear ElMaestro, dear All!

Also suffering from senile "Bettflucht" (insomnia) here is my quick shoot with Owen's Q:
(forgive me, I'm a very beginning adept of R)

# The functions of normal-, t-distributions and integrate()
require(stats)
# -----------------------------------------------------------------------------
# Owen's Q-function
OwensQ <- function (nu, t, delta, a, b)
{
   # .Machine$double.eps^.5 = 1.490116e-08 on my machine
   # What method integrate() uses is not clear to me. But seems sophisticated
   # to assure enough accuracy
   Qintegral <- integrate(Q.integrand, lower = a, upper = b,
                nu=nu, t=t, delta = delta,
                subdivisions = 1000,
                rel.tol = .Machine$double.eps^.5, abs.tol = .Machine$double.eps^.5,
                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 ?     
Q.integrand <- function(x, nu, t, delta)
{
   # assure a vector answer, copied from Power.density.md MBESS
   # look for better implementation! is loop necessary?
   pdim <- length(x)
   dens <- c(1:pdim)
   # constants for the loop
   sqrtnu   <- sqrt(nu)
   lnQconst <- 0.5*log(2*pi)-lgamma(nu/2)-((nu-2)/2)*log(2)

   for (ix in 1:pdim) {
     xx <- x[ix]
     if (xx > 0) {
       d <- pnorm(t*xx/sqrtnu - delta, mean = 0, sd = 1, log = TRUE) +
            (nu-1)*log(xx) + dnorm(xx, mean = 0, sd = 1, log = TRUE)
       dens[ix] <- exp(lnQconst + d)
     } else {
       # do we need  < 0 ? or must xx>=0? for our power it is.
       dens[ix] <-0
     }
   }

   return(dens)
}

# -----------------------------------------------------------------------------
# Power of two-one-sided-t-tests using OwensQ
# In case of logscale=TRUE give ldiff, ltheata1 and ltheta2 as ratios
# f.i. ldiff=0.95, ltheta1=0.8, ltheta2=1.25
# in case of logscale=FALSE give ldiff, ltheata1 and ltheta2 as difference to 1
# f.i. ldiff=0.05 (5% difference) ltheata1=-0.2, ltheta2=0.2 (20% equivalence margins)
# CV is always the coefficient of variation but as ratio, not %

Power.TOST <- function(alpha, logscale, ltheta1, ltheta2, ldiff, CV, n, df)
{
    theta1 <- ltheta1
    theta2 <- ltheta2
    diff   <- ldiff
    s      <- CV
    if (logscale) {
        theta1 <- log(ltheta1)
        theta2 <- log(ltheta2)
        diff   <- log(ldiff)
        s <- sqrt(log(1+CV*CV))
    }
    tval <- qt(1 - alpha, df, lower.tail = TRUE, log.p = FALSE)
    delta1 <- (diff-theta1)/(s*sqrt(2/n))
    delta2 <- (diff-theta2)/(s*sqrt(2/n))
    R <- (delta1-delta2)*sqrt(df)/(2*tval)
    p1 <- OwensQ(df,  tval, delta1, 0, R)
    p2 <- OwensQ(df, -tval, delta2, 0, R)

    power=p2-p1

    return(power)
}

Some spare 'validation' results (Sample sizes from Diletti et. al. for target power 70%, 80%, 90%):
CV=0.3, BE margins 0.8-1.25, Null ratio=0.95, alpha=0.05
   n power(Q) Apfelstrudel
1 32 0.7177182  0.7177182
2 40 0.8158453  0.8158453
3 52 0.9019652  0.9019652

PS: I think it's enough R power in this thread to have it in the R category now.

Edit: Category changed. require(stats) not needed, because loaded by default. [Helmut]
—
Regards,

Detlew

Complete thread:

UA Flag
Activity
 Admin contact
23,424 posts in 4,927 threads, 1,671 registered users;
21 visitors (0 registered, 21 guests [including 9 identified bots]).
Forum time: 14:40 CEST (Europe/Vienna)

Truth and clarity are complementary.    Niels Bohr

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