Nothing
#---------------------------------------------------------------------------
# Sample size for partial and full replicate design with BE decision
# via ABE + sratio <=2.5 according to the FDA Dabigatran / rivaroxaban guidances
# via simulated (empirical) power
#
# Author: dlabes
#---------------------------------------------------------------------------
sampleN.HVNTID <- function(alpha=0.05, targetpower=0.8, theta0, theta1, theta2,
CV, design=c("2x2x4", "2x2x3"), nsims=1E5, nstart,
imax=100, print=TRUE, details=TRUE, setseed=TRUE)
{
if (missing(theta0)) theta0 <- 0.95
if (missing(theta1) & missing(theta2)) theta1 <- 0.8
if (missing(theta2)) theta2 <- 1/theta1
if (missing(theta1)) theta1 <- 1/theta2
# 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)
# CV's for Test/Reference
CVwT <- CV[1]
if (length(CV)>1) CVwR <- CV[2] else CVwR <- CVwT
if (length(CV)>2) warning("Only first 2 entries from CV vector used.")
s2wT <- CV2mse(CVwT)
s2wR <- CV2mse(CVwR)
# check design argument
design <- match.arg(design)
# we are treating only balanced designs
# thus we use here bk - the design constant for ntotal
# expressions for the df's
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)
dfTTe <- parse(text="n-2", srcfile=NULL)
# expectation of mse of the ANOVA of intra-subject contrasts T-R
Emse <- (s2wT + s2wR)/2
desi <- "2x2x4 (TRTR|RTRT)"
}
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) # balanced only
dfTTe <- parse(text="n/2-1", srcfile=NULL) # balanced only
# expectation of mse of the ANOVA of intra-subject contrasts T-R
Emse <- 1.5*(s2wT + s2wR)/2 # balanced only
desi <- "2x2x3 (TRT|RTR)"
}
mlog <- log(theta0)
ltheta1 <- log(theta1)
ltheta2 <- log(theta2)
# we have eventually to adapt this empirically for the n0 start value
bkk <- bk
if (print){
cat("\n+++++++++ FDA method for HV NTIDs ++++++++++++\n")
cat(" Sample size estimation\n")
cat("----------------------------------------------\n")
cat("Study design:",desi,"\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")
cat("ABE limits =", theta1, "...", theta2,"\n")
# if (details) {
# cat("Implied scABEL =", formatC(exp(ltheta1), format="f", digits=4),
# "...", formatC(exp(ltheta2), format="f", digits=4),"\n")
# }
# cat("Regulatory settings:",regulator,"\n")
# if (details) {
# cat("- Regulatory const. =",r_const,"\n")
# # CVcap?
# cat("- 'CVcap' =", formatC(CVcap, format="f", digits=4),"\n")
# }
}
# -----------------------------------------------------------------
# nstart? from sampleN0 with shrunken/widened limits
# does'nt fit always really good
if (missing(nstart)){
n <- .sampleN0(alpha=alpha, targetpower, ltheta1, ltheta2, diffm=mlog,
se=sqrt(Emse), steps=seqs, bk=bkk, diffmthreshold=0.01)
#cat("n0=",n,"\n")
} else n <- seqs*round(nstart/seqs)
nmin <- 6
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)
dfTT <- eval(dfTTe)
if(setseed) set.seed(123456)
p <- .power.HVNTID(mlog, sdm, C3, Emse, df, s2wR, dfRR, s2wT, dfTT, nsims,
ln_lBEL=log(theta1),ln_uBEL=log(theta2), alpha)
pwr <- as.numeric(p["BE"]);
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=6, format="f"),"\n")
}
iter <- 0
# iter>100 is emergency brake
# --- loop until power <= target power, step-down
down <- FALSE; up <- FALSE
while (pwr>targetpower) {
if (n<=nmin) {
if (details & iter==0) cat( n," ", formatC(pwr, digits=6, format="f"),"\n")
break
}
down <- TRUE
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)
dfTT <- eval(dfTTe)
if(setseed) set.seed(123456)
p <- .power.HVNTID(mlog, sdm, C3, Emse, df, s2wR, dfRR, s2wT, dfTT, nsims,
ln_lBEL=log(theta1),ln_uBEL=log(theta2), alpha)
pwr <- as.numeric(p["BE"]);
# do not print first step down
if (details) cat( n," ", formatC(pwr, digits=6, 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) {
up <- TRUE ; down <- FALSE
n <- n+seqs # step up
iter <- iter+1
C3 <- 1/n
sdm <- sqrt(Emse*C3)
df <- eval(dfe)
dfRR <- eval(dfRRe)
dfTT <- eval(dfTTe)
if(setseed) set.seed(123456)
p <- .power.HVNTID(mlog, sdm, C3, Emse, df, s2wR, dfRR, s2wT, dfTT, nsims,
ln_lBEL=log(theta1),ln_uBEL=log(theta2), alpha)
pwr <- as.numeric(p["BE"]);
if (details) cat( n," ", formatC(pwr, digits=6, 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=6, 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")
if (print | details) return(invisible(res)) else return(res)
} # end function
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.