Nothing
#------------------------------------------------------------------------------
# functions for power calculation
# Author: dlabes
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
# helper function to allow partial match of the method
.powerMethod <- function(method){
meth <- tolower(method[1])
if (method=="") method <- "exact"
methods <- c("exact","owenq","noncentral","nct","shifted","central","mvt")
# ^ = match at start
meth <- methods[grep(paste("^",meth,sep=""), methods)]
if (length(meth)==0) meth <- tolower(method)
if (length(meth)>1) meth <- "nct" # happens if only "n" is given
return(meth)
}
#------------------------------------------------------------------------------
# 'raw' power function without any error checks,
# does not vectorize propperly!
# to be used in sampleN.TOST avoiding overhead of redundant calculations
# in case of multiplicative model:
# diffm=log(null ratio), theta1=log(lower BE limit), theta2=log(upper BE limit)
# in case of additive model:
# diffm=1-null ratio, theta1=lower BE limit-1, theta2=upper BE limit -1
# Jan 2015: Interface changed to sem
# so call it with sem= se*sqrt(bk/n) if balanced or se*sqrt(bkni*sum(1/n))
.power.TOST <- function(alpha=0.05, ltheta1, ltheta2, diffm, sem, df)
{
stopifnot(length(alpha) <= 2)
tval <- qt(1 - alpha, df, lower.tail = TRUE)
dl <- length(tval)
# 0/0 -> NaN in case diffm=ltheta1 or diffm=ltheta2 and sem=0!
delta1 <- (diffm-ltheta1)/sem
delta2 <- (diffm-ltheta2)/sem
# is this correct?
delta1[is.nan(delta1)] <- 0
delta2[is.nan(delta2)] <- 0
# Define integration boundary R for OwensQ
R <- (delta1-delta2)*sqrt(df)/(tval[1]+tval[dl])
# in case of se=0 it results: delta1=Inf, delta2=inf if diffm>ltheta2
# Inf - Inf is NaN
R[is.nan(R)] <- 0
# If delta1 > delta2, or equivalently ltheta1 < ltheta2:
# R is negative if and only if tval[1] <= -tval[dl].
# An example would be equal alphas for the two hypotheses and alpha > 0.5
# (very unusual!) => the upper integration limit is lower then the
# lower limit! SAS OwenQ gives missings if b or a are negative!
# On the other hand SAS Proc Power gives values which are seemingly calculated
# with abs(R).
# Correct acc. to Fig. 1 given in K.Philips
# "Power for Testing Multiple Instances of the Two One-Sided Tests Procedure"
# The International Journal of Biostatistics: Vol. 5: Iss. 1, Article 15.
# DOI: 10.2202/1557-4679.1169
# should be R=Inf, i.e. unlimited integration with respect to sigma.
# This gives the same values (within certain precision) as Ben's power.1TOST
# aka power.TOST(..., method="mvt").
# Can also be checked via function power.TOST.sim().
# Also gives the same results as Power code (base.powerfun.TOST) from
# Maurer et al (2018)
R[R<=0] <- Inf
# to check SAS Proc power values comment above out and write
# R <- abs(R)
# to avoid numerical errors in OwensQ implementation
if (min(df)>10000){
# 'shifted' normal approximation Jan 2015
# former Julious formula (57)/(58) doesn't work
tval <- qnorm(1-alpha)
p1 <- pnorm(tval[1]-delta1)
p2 <- pnorm(-tval[dl]-delta2)
# may give negative values
# thus set to zero
pwr <- p2-p1
pwr[pwr<0] <- 0
return(pwr)
}
if (min(df)>=5000 & min(df<=10000)) {
# approximation via non-central t-distribution
return(.approx.power.TOST(alpha, ltheta1, ltheta2, diffm, sem, df))
}
# attempt to vectorize (it vectorizes properly if diffm is a vector
# OR se OR n,df are vectors)
p1 <- vector("numeric", length(delta1))
p2 <- vector("numeric", length(delta1))
for (i in seq_along(delta1)) {
# get correct values for p1 and p2:
# p1 is for left hypothesis (right-tailed)
# -> critical value is always first component
# p2 for right hypothesis (left-tailed)
# -> critical value is first component if alpha is 1-dim,
# second component if alpha is 2-dim
p1[i] <- OwensQ(df[1], tval[1], delta1[i], 0, R[i])
p2[i] <- OwensQ(df[1], -tval[dl], delta2[i], 0, R[i])
}
pwr <- p2-p1
# due to numeric inaccuracies power < 0
# paranoia
pwr[pwr<0] <- 0
return( pwr )
}
#------------------------------------------------------------------------------
# Ben's implementation of power via integration of the bivariate t-distribution
# with correlation == 1, also exact
# does'nt vectorize in any respect!
.power.1TOST <- function(alpha, ltheta1, ltheta2, diffm, sem, df, setseed = TRUE)
{
stopifnot(length(alpha) <= 2)
if (setseed) set.seed(123456)
corr <- matrix(1, ncol = 2, nrow = 2)
tval <- qt(1 - alpha, df)
lower <- c(tval[1], -Inf)
upper <- c(Inf, -tval[length(tval)])
delta1 <- (diffm - ltheta1) / sem
delta2 <- (diffm - ltheta2) / sem
pow <- rep(0, times=length(delta1))
# attempt to vectorize if ltheta0 OR se is a vector
for(i in seq_along(delta1)){
delta <- c(delta1[i], delta2[i])
prob <- pmvt(lower = lower, upper = upper, delta = delta, df = df, corr = corr,
algorithm = GenzBretz(maxpts=100000, abseps = 1e-05))#[1]
# abseps=1e-6 gives often "Completion with error > abseps"
# give a warning if attr(prob,"msg") not equal "Normal completion"?
if(attr(prob, which="msg")!="Normal Completion")
warning("pmvt returned message ", attr(prob, which="msg"), call.=FALSE)
pow[i] <- prob[1]
}
pow
}
#------------------------------------------------------------------------------
# 'raw' approximate power function without any error checks,
# approximation based on non-central t
# this vectorizes ok
.approx.power.TOST <- function(alpha=0.05, ltheta1, ltheta2, diffm, sem, df)
{
stopifnot(length(alpha) <= 2)
tval <- qt(1 - alpha, df, lower.tail = TRUE, log.p = FALSE)
# 0/0 -> NaN in case diffm=ltheta1 or diffm=ltheta2
# and sem=0!
delta1 <- (diffm-ltheta1)/sem
delta2 <- (diffm-ltheta2)/sem
# is this correct?
delta1[is.nan(delta1)] <- 0
delta2[is.nan(delta2)] <- 0
# suppress warnings with regard to insufficient precision of nct
pow <- suppressWarnings(pt(-tval[length(tval)], df, ncp=delta2) -
pt(tval[1], df, ncp=delta1))
pow[pow<0] <- 0 # this is to avoid neg. power due to approx. (vector form)
return(pow)
}
#------------------------------------------------------------------------------
# 'raw' power function without any error checks,
# approximation based on central 'shifted' central t distribution
# according to Chow, Liu "Design and Analysis of Bioavailability ..."
# Chapter 9.6 and implemented in PASS 2008
# where does this all come from?
.approx2.power.TOST <- function(alpha=0.05, ltheta1, ltheta2, diffm, sem, df)
{
tval <- qt(1 - alpha, df, lower.tail = TRUE)
# 0/0 -> NaN in case diffm=ltheta1 or diffm=ltheta2
# and se=0!
delta1 <- (diffm-ltheta1)/sem
delta2 <- (diffm-ltheta2)/sem
# is this correct?
delta1[is.nan(delta1)] <- 0
delta2[is.nan(delta2)] <- 0
pow <- pt(-tval[length(tval)]-delta2, df) - pt(tval[1]-delta1, df)
pow[pow<0] <- 0 # this is to avoid neg. power due to approx. (vector form)
return(pow)
}
#------------------------------------------------------------------------------
# function for merging the various power calculations
.calc.power <- function(alpha=0.05, ltheta1, ltheta2, diffm, sem, df, method="exact")
{
pow <- switch(
method,
exact=.power.TOST(alpha, ltheta1, ltheta2, diffm, sem, df),
owenq=.power.TOST(alpha, ltheta1, ltheta2, diffm, sem, df),
mvt= .power.1TOST(alpha, ltheta1, ltheta2, diffm, sem, df),
nct= .approx.power.TOST(alpha, ltheta1, ltheta2, diffm, sem, df),
noncentral=.approx.power.TOST(alpha, ltheta1, ltheta2, diffm, sem, df),
shifted=.approx2.power.TOST(alpha, ltheta1, ltheta2, diffm, sem, df),
central=.approx2.power.TOST(alpha, ltheta1, ltheta2, diffm, sem, df),
stop("Method '", method, "' unknown!\n", call.=TRUE)
)
return(pow)
}
#------------------------------------------------------------------------------
# Power of two-one-sided-t-tests
# (this is a wrapper to the various power calculation methods)
# In case of logscale=TRUE give theta0, theata1 and theta2 as ratios
# f.i. theta0=0.95, theta1=0.8, theta2=1.25
# In case of logscale=FALSE give theta0, theata1 and theta2 as difference
# to 1 f.i. theta0=0.05 (5% difference),
# theata1=-0.2, theta2=0.2 20% equivalence margins)
# CV is always the coefficient of variation but as ratio, not %
# leave upper BE margin (ltheta2) empty and the function will use -lower
# in case of additive model or 1/lower if logscale=TRUE
power.TOST <- function(alpha=0.05, logscale=TRUE, theta1, theta2, theta0,
CV, n, design="2x2", method="exact", robust=FALSE)
{
if (missing(CV)) stop("CV must be given!")
if (missing(n)) stop("Number of subjects n must be given!")
if (length(alpha) != 1) stop("alpha must be a scalar!")
# check if design is implemented
d.no <- .design.no(design)
if (is.na(d.no)) stop("Design ",design, " unknown!", call.=FALSE)
# design characteristics
ades <- .design.props(d.no)
# degrees of freedom as expression
dfe <- .design.df(ades, robust=robust)
# design const.
# we use always bkni
#bk <- ades$bk
# step size = no of sequences
steps <- ades$steps
# handle n = ntotal if scalar else n's of the sequence groups
if (length(n) == 1) {
# total n given
# for unbalanced designs we divide the ns by ourself
# to have only small imbalance (function nvec() from Helper_dp.R)
if(is.finite(n)) n <- nvec(n=n, grps=ades$steps) else n <- rep(Inf, times=steps)
if (n[1]!=n[length(n)]){
message("Unbalanced design. n(i)=", paste(n, collapse="/"), " assumed.")
}
} else {
if (length(n) != ades$steps) {
stop("Length of n vector must be ", ades$steps, "!")
}
if (any(n<1)) stop("All n(i) have to be >0.")
}
nc <- sum(1/n)
n <- sum(n)
se.fac <- sqrt(ades$bkni * nc)
# regularize the method giving
method <- .powerMethod(method)
# handle log-transformation
if (logscale) {
if (missing(theta0)) theta0 <- 0.95
if (missing(theta1) & missing(theta2)) theta1 <- 0.80
if (missing(theta1)) theta1 = 1/theta2
if (missing(theta2)) theta2 = 1/theta1
ltheta1 <- log(theta1)
ltheta2 <- log(theta2)
ldiff <- log(theta0)
se <- CV2se(CV)
} else { # untransformed
if (missing(theta0)) theta0 <- 0.05
if (missing(theta1) & missing(theta2)) theta1 <- -0.2
if (missing(theta1)) theta1 <- -theta2
if (missing(theta2)) theta2 <- -theta1
ltheta1 <- theta1
ltheta2 <- theta2
ldiff <- theta0
se <- CV
}
if (length(CV) > 1 & length(theta0) > 1)
stop("Only CV or theta0 can be vectors, not both of them")
df <- eval(dfe)
if (any(df<1)) stop("n too small. Degrees of freedom <1!")
pow <- .calc.power(alpha, ltheta1, ltheta2, ldiff, sem=se*se.fac, df,
method=method)
return( pow )
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.