#---------------------------------------------------------------------------
# unified function
# chooses the sample size function according to regulator$est_method
# former is now sampleN.scABEL1
#
# author dlabes
#---------------------------------------------------------------------------
sampleN.scABEL <- function(alpha=0.05, targetpower=0.8, theta0, theta1,
theta2, CV, design=c("2x3x3", "2x2x4", "2x2x3"),
regulator, nsims=1E5, nstart, imax=100, print=TRUE,
details=TRUE, setseed=TRUE)
{
# design must be checked outside
desi <- match.arg(design)
# check regulator
if (missing(regulator)) regulator <- "EMA"
reg <- reg_check(regulator)
ssfun <- "sampleN.scABEL1"
if (reg$est_method=="ISC") ssfun <- "sampleN.scABEL2"
# print(sys.call())
# next doesn't function if arguments are missing
# r <- do.call(ssfun,
# list(alpha, targetpower, theta0, theta1, theta2, CV,
# design=desi, regulator=reg, nsims, nstart, imax,
# print, details, setseed))
if (reg$est_method!="ISC"){
r <- sampleN.scABEL1(alpha, targetpower, theta0, theta1, theta2, CV,
design=desi, regulator=reg, nsims, nstart, imax,
print, details, setseed)
} else {
r <- sampleN.scABEL2(alpha, targetpower, theta0, theta1, theta2, CV,
design=desi, regulator=reg, nsims, nstart, imax,
print, details, setseed)
}
if(print | details) return(invisible(r)) else return(r)
}
#-----------------------------------------------------------------------------
# Sample size for partial and full replicate design and scaled ABE
# via simulated (empirical) power
# estimation method ANOVA and BE decision via ABEL (Average bioequivalence with
# expanding limits)
#
# Author: dlabes
#-----------------------------------------------------------------------------
# helper function: sample size for pe in a range?
# definition is more or less empirical (i.e. not understood by me)
.sampleN0.2 <- function(targetpower, ltheta2, diffm, se, bk, steps)
{
n <- qnorm(targetpower)^2*se^2*bk/(abs(diffm)-ltheta2)^2
#make an even multiple of step (=2 in case of 2x2 cross-over)
n <- steps*trunc(n/steps)
n
}
sampleN.scABEL1 <- function(alpha=0.05, targetpower=0.8, theta0, theta1,
theta2, CV, design=c("2x3x3", "2x2x4", "2x2x3"),
regulator, nsims=1E5, nstart, imax=100, print=TRUE,
details=TRUE, setseed=TRUE)
{
if (missing(theta1) & missing(theta2)) theta1 <- 0.8
if (missing(theta1)) theta1 <- 1/theta2
if (missing(theta2)) theta2 <- 1/theta1
# the two Laszlos recommend theta0=0.9 for HVDs
if (missing(theta0)) theta0 <- 0.9
# theta0 in range
if ( (theta0 < theta1) | abs(theta0-theta1) <1e-8
| (theta0 > theta2) | abs(theta0-theta2) <1e-8 ){
stop("True ratio ", theta0," not within margins ", theta1," ... ",
theta2,"!", call.=FALSE)
}
if (missing(CV)) stop("CV must be given!", call.=FALSE)
# subject-by-formulation interaction can't play a role here I think
# since the EMA model doesn't allow such term
CVwT <- CV[1]
# should we allow different variabilities in the EMA method at all?
if (length(CV)==2) CVwR <- CV[2] else CVwR <- CVwT
s2wT <- log(1.0 + CVwT^2)
s2wR <- log(1.0 + CVwR^2)
if(missing(regulator)) regulator <- "EMA"
# check regulator and get
# constants acc. to regulatory bodies (function in scABEL.R)
reg <- reg_check(regulator)
CVcap <- reg$CVcap
CVswitch <- reg$CVswitch
r_const <- reg$r_const
pe_constr <- reg$pe_constr
# check design
design <- match.arg(design)
# we are treating only balanced designs
# thus we use here bk - design constant for ntotal
# expressions for the df's
if (design=="2x3x3") {
bk <- 1.5; seqs <- 3
dfe <- parse(text="2*n-3", srcfile=NULL)
dfRRe <- parse(text="n-2", srcfile=NULL)
#sd2 <- (s2wT + s2wR)/2 # used in v1.1-00 - v1.1-02, wrong
# simulations with s2D=0 show:
Emse <- (s2wT + 2.0*s2wR)/3
}
if (design=="2x2x4") {
bk <- 1.0; seqs <- 2
# only EMA settings
dfe <- parse(text="3*n-4", srcfile=NULL)
dfRRe <- parse(text="n-2", srcfile=NULL)
# sd^2 (variance) of the differences T-R from their components
Emse <- (s2wT + s2wR)/2
}
if (design=="2x2x3") {
bk <- 1.5; seqs <- 2
# only EMA settings
dfe <- parse(text="2*n-3", srcfile=NULL)
dfRRe <- parse(text="n/2-1", srcfile=NULL)
# sd^2 (variance) of the differences T-R from their components
Emse <- (s2wT + s2wR)/2 # for balanced designs we use here
}
mlog <- log(theta0)
if (print){
designs <- c("2x2x4", "2x2x3", "2x3x3")
type <- c("4 period full replicate",
"3 period full replicate",
"partial replicate") # clear words
cat("\n+++++++++++ scaled (widened) ABEL +++++++++++\n")
cat(" Sample size estimation\n")
cat(" (simulation based on ANOVA evaluation)\n")
cat("---------------------------------------------\n")
cat("Study design: ")
cat(paste0(design, " (", type[match(design, designs)], ")\n"))
cat("log-transformed data (multiplicative model)\n")
cat(nsims,"studies for each step simulated.\n\n")
cat("alpha = ", alpha,", target power = ", targetpower,"\n", sep="")
cat("CVw(T) = ", CVwT,"; CVw(R) = ", CVwR,"\n", sep="")
cat("True ratio = ",theta0,"\n", sep="")
cat("ABE limits / PE constraint =", theta1, "...", theta2, "\n")
if (reg$name == "GCC" & CVwR > 0.3)
cat("Widened limits =", 0.75, "...", 1/0.75, "\n")
if (details | reg$name=="USER") {
print(reg)
cat("\n")
} else {
cat("Regulatory settings:", reg$name,"\n")
}
}
# -----------------------------------------------------------------
# nstart? from sampleN0 with widened limits
# doesn't fit really good if theta0>=1.2! ways out? see sampleN0.2
ltheta1 <- -sqrt(s2wR)*r_const
ltheta2 <- -ltheta1
if (CVwR <= CVswitch){
ltheta1 <- log(theta1)
ltheta2 <- log(theta2)
}
if (CVwR > CVcap){
ltheta1 <- -sqrt(log(1.0 + CVcap^2))*r_const
ltheta2 <- -ltheta1
}
if (missing(nstart)){
# start from ABE start with widened limits
n01 <- .sampleN0_2(alpha=alpha, targetpower, ltheta1, ltheta2, diffm=mlog,
se=sqrt(Emse), steps=seqs, bk=bk)
# empirical correction in the vicinity of CV=0.3 for ratios
# outside 0.86 ... 1/0.86
if(Emse < CV2mse(CVswitch+0.005) & Emse > CV2mse(CVswitch-0.005)
& abs(mlog)>log(1/0.865)) {
if (reg$name=="EMA") n01 <- 0.9*n01
if (reg$name=="FDA") n01 <- 0.65*n01
n01 <- seqs*trunc(n01/seqs)
}
# start from PE constraint sample size
n02 <- .sampleN0.2(targetpower, ltheta2=log(theta2), diffm=mlog,
se=sqrt(Emse), steps=seqs, bk=bk)
# debug print
# cat(n01,n02,"\n")
n <- max(c(n01,n02))
} else n <- seqs*round(nstart/seqs)
nmin <- 6
if (n<nmin) n <- nmin
# we are simulating for balanced designs
C2 <- bk/n
# sd of the sample mean T-R (point estimator)
sdm <- sqrt(Emse*C2)
df <- eval(dfe)
dfRR <- eval(dfRRe)
dfTT <- dfRR
if(setseed) set.seed(123456)
p <- .pwr.ABEL.ANOVA(mlog, sdm, C2, Emse, df, s2wR, dfRR, s2wT, dfTT,
design, nsims, CVswitch, r_const, CVcap, pe_constr,
ln_lBEL=log(theta1),ln_uBEL=log(theta2), alpha=alpha)
pwr <- as.numeric(p["BE"]);
pd <- max(4,round(log10(nsims),0)-1) # digits for power
if (details) {
cat("\nSample size search\n")
cat(" n power\n")
# do not print first too high
# this is for cases with only one step-down and than step up
if (pwr<=targetpower) cat( n," ", formatC(pwr, digits=pd, format="f"),"\n")
}
iter <- 0
# iter>100 is emergency brake
# --- loop until power <= target power, step-down
down <- FALSE
up <- FALSE
while (pwr>targetpower) {
down <- TRUE
if (n<=nmin) {
if (details & iter==0) cat( n," ", formatC(pwr, digits=pd, format="f"),"\n")
break
}
n <- n-seqs # step down if start power is to high
iter <- iter + 1
C2 <- bk/n
# sd of the sample mean T-R (point estimator)
sdm <- sqrt(Emse*C2)
df <- eval(dfe)
dfRR <- eval(dfRRe)
dfTT <- dfRR
if(setseed) set.seed(123456)
p <- .pwr.ABEL.ANOVA(mlog, sdm, C2, Emse, df, s2wR, dfRR, s2wT, dfTT,
design, nsims, CVswitch, r_const, CVcap, pe_constr,
ln_lBEL=log(theta1),ln_uBEL=log(theta2), alpha=alpha)
pwr <- as.numeric(p["BE"]);
if (details) cat( n," ", formatC(pwr, digits=pd, format="f"),"\n")
if (iter>imax) break
# loop results in n with power too low
# must step up again one step. is done in the next loop
}
while (pwr<targetpower) {
up <- TRUE; down <- FALSE
n <- n+seqs # step-up
iter <- iter+1
C2 <- bk/n
sdm <- sqrt(Emse*C2)
df <- eval(dfe)
dfRR <- eval(dfRRe)
dfTT <- dfRR
if(setseed) set.seed(123456)
p <- .pwr.ABEL.ANOVA(mlog, sdm, C2, Emse, df, s2wR, dfRR, s2wT, dfTT,
design, nsims, CVswitch, r_const, CVcap, pe_constr,
ln_lBEL=log(theta1),ln_uBEL=log(theta2), alpha=alpha)
pwr <- as.numeric(p["BE"]);
if (details) cat( n," ", formatC(pwr, digits=pd, format="f"),"\n")
if (iter>imax) break
}
nlast <- n
if (up & pwr<targetpower) {
n <- NA
if (details) cat("Sample size search failed!\n")
}
if (down & pwr>targetpower & n != nmin) {
n <- NA
if (details) cat("Sample size search failed!\n")
}
if (print && !details) {
cat("\nSample size\n")
cat(" n power\n")
cat( n," ", formatC(pwr, digits=pd, format="f"),"\n")
if (is.na(n)) cat("Sample size search failed!\n")
}
if (print) cat("\n")
#return results as data.frame
res <- data.frame(design=design, alpha=alpha, CVwT=CVwT, CVwR=CVwR,
theta0=theta0, theta1=theta1, theta2=theta2, n=n, power=pwr,
targetpower=targetpower, nlast=nlast)
names(res) <-c("Design", "alpha", "CVwT", "CVwR", "theta0", "theta1", "theta2",
"Sample size", "Achieved power", "Target power", "nlast")
#cat("iter=",iter,"\n")
if (print | details) return(invisible(res)) else return(res)
} # end function
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.