#---------------------------------------------------------------------------
# Sample size for partial and full replicate design and scaled ABE
# via simulated (empirical) power
# estimation method via intra-subject contrasts & BE decision via
# linearized reference scale BE criterion, no cap
#
# Author: dlabes
#---------------------------------------------------------------------------
sampleN.RSABE <- function(alpha=0.05, targetpower=0.8, theta0, theta1,
theta2, CV, design=c("2x3x3", "2x2x4", "2x2x3"),
regulator = c("FDA", "EMA"), nsims=1E5, nstart, imax=100,
print=TRUE, details=TRUE, setseed=TRUE) {
if (missing(theta1) & missing(theta2)) theta1 <- 0.8
if (missing(theta2)) theta2=1/theta1
if (missing(theta1)) theta1=1/theta2
# according to the two Laszlos' paper 0.9 (alas 1.10) should be considered
if (missing(theta0)) theta0 <- 0.90
# 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(s) must be given!", call.=FALSE)
# regulator here only FDA, EMA
# other regulatory bodies ("HC", "ANVISA") use all the EMA regulatory constant
if (missing(regulator)) regulator <- "FDA"
reg <- reg_check(regulator, choices=c("FDA", "EMA"))
CVswitch <- reg$CVswitch
r_const <- reg$r_const
pe_constr <- reg$pe_constr
# CVcap doesn't apply to the FDA recommended method
# but in Munoz et al. method= Howe-EMA
CVcap <- reg$CVcap
# for later enhancement taking into account the
# subject-by-formulation interaction
s2D <- 0
CVwT <- CV[1]
if (length(CV)==2) CVwR <- CV[2] else CVwR <- CVwT
s2wT <- log(1.0 + CVwT^2)
s2wR <- log(1.0 + CVwR^2)
# 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") {
seqs <- 3
bk <- 1.5 # needed for n0?
# in case of the FDA we are using the 'robust' df's
# due to the fact that the described analysis in the
# progesterone guidance is based in the intrasubject contrasts
# T-R and R-R with df=n-seqs
dfe <- parse(text="n-3", srcfile=NULL)
dfRRe <- parse(text="n-3", srcfile=NULL)
# expectation of mse of the ANOVA of intra-subject contrasts
#sd2 <- s2D + (s2wT + s2wR)/2 # used in v1.1-00 - v1.1-02
# according to McNally et al., verified via simulations:
Emse <- s2D + s2wT + s2wR/2
}
if (design == "2x2x4") {
seqs <- 2
bk <- 1.0 # needed for n0
dfe <- parse(text="n-2", srcfile=NULL)
dfRRe <- parse(text="n-2", srcfile=NULL)
# expectation of mse of the ANOVA of intra-subject contrasts
Emse <- (s2D + (s2wT + s2wR)/2)
}
if (design=="2x2x3") {
seqs <- 2
bk <- 1.5 # needed for n0?
dfe <- parse(text="n-2", srcfile=NULL)
dfRRe <- parse(text="n/2-1", srcfile=NULL) # for balanced designs
# expectation of mse of the ANOVA of intra-subject contrasts
Emse <- 1.5*(s2wT + s2wR)/2 # for balanced designs
}
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++++++++ Reference scaled ABE crit. +++++++++\n")
cat(" Sample size estimation\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 constraints =",theta1,"...", theta2,"\n")
if (details | reg$name=="USER") {
reg$CVcap <- NULL # CVcap doesn't apply here
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 or <=0.85! ways out?
ltheta1 <- -r_const*sqrt(s2wR)
ltheta2 <- -ltheta1
# this does not function in case of CVwR=0.3 for the original code
# calculating s2wR and back-calculating CVwR from that
# numerical problem?
if (CVwR <= CVswitch){
ltheta1 <- log(theta1)
ltheta2 <- log(theta2)
}
if (missing(nstart)){
# try to use the empirical alpha for start sample size, some sort of
al <- alpha
if (reg$name=="FDA") {
if(Emse/bk <= CV2mse(0.30001) & Emse/bk >= CV2mse(0.2975)) al=0.12
if(Emse/bk > CV2mse(0.30001)) al <- 0.035
}
if (reg$name=="EMA") {
#does not fit!
#if(Emse/bk <= CV2mse(0.321) & Emse/bk >= CV2mse(0.28)) al=0.065
}
# debug print
# cat(al,"\n")
# we use bk=1 here since our formula is sem=sqrt(Emse/n)
n01 <- .sampleN0(alpha=al, targetpower, ltheta1, ltheta2, diffm=mlog,
se=sqrt(Emse), steps=seqs, bk=1, diffmthreshold=0.01)
# empirical correction in the vicinity of CV=0.3, for ratios
# outside 0.86 ... 1/0.86
# if(Emse/bk <= CV2mse(0.305) & Emse/bk >= CV2mse(0.295) & abs(mlog)>log(1/0.865)) {
# if (regulator=="EMA") n01 <- 0.9*n01 else n01 <- 0.8*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=1)
# debug print
# cat(n01,n02,"\n")
n <- max(c(n01,n02))+seqs
} else n <- seqs*round(nstart/seqs)
nmin <- 6 # fits 2x3x3 and 2x2x4
if(n<nmin) n <- nmin
# we are simulating for balanced designs
C3 <- 1/n
# sd of the sample mean T-R (point estimator)
sdm <- sqrt(Emse*C3)
df <- eval(dfe)
dfRR <- eval(dfRRe)
if(setseed) set.seed(123456)
p <- .power.RSABE(mlog, sdm, C3, Emse, df, s2wR, dfRR, nsims, CVswitch,
r_const, pe_constr, CVcap,
ln_lBEL=log(theta1), ln_uBEL=log(theta2), alpha=alpha)
pwr <- as.numeric(p["BE"]);
pd <- max(4,round(log10(nsims),0)) # 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
C3 <- 1/n
# sd of the sample mean T-R (point estimator)
sdm <- sqrt(Emse*C3)
df <- eval(dfe)
dfRR <- eval(dfRRe)
if(setseed) set.seed(123456)
p <- .power.RSABE(mlog, sdm, C3, Emse, df, s2wR, dfRR, nsims, CVswitch,
r_const, pe_constr, CVcap, ln_lBEL=log(theta1),
ln_uBEL=log(theta2), alpha=alpha)
pwr <- as.numeric(p["BE"]);
# do not print first step down
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 one up again. is done in the next loop
}
while (pwr<targetpower) {
down <- FALSE; up <- TRUE
n <- n+seqs # step up
iter <- iter+1
C3 <- 1/n
sdm <- sqrt(Emse*C3)
df <- eval(dfe)
dfRR <- eval(dfRRe)
if(setseed) set.seed(123456)
p <- .power.RSABE(mlog, sdm, C3, Emse, df, s2wR, dfRR, nsims, CVswitch,
r_const, pe_constr, CVcap, 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")
# debug print
# cat("iter=",iter+1,"\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.