Nothing
#require(PowerTOST)
#-------------------------------------------------------------------------
# power analysis of a sample size plan for scaled ABE (EMA & FDA)
# Input area at the end of the code.
# coded originally by Helmut Schuetz
# adapted to PowerTOST infra structure by D. Labes
#-------------------------------------------------------------------------
pa.scABE <- function(CV, theta0=0.9, targetpower=0.8, minpower=0.7,
design=c("2x3x3", "2x2x4", "2x2x3"),
regulator=c("EMA", "HC", "FDA", "GCC"), ...)
{
# Rversion must be >=3.1.0 for the uniroot call with argument extendInt
Rver <- paste0(R.Version()$major, ".", R.Version()$minor)
if (length(CV)>1) stop("Only CVwT=CVwR=CV implemented. CV has to be a scalar.")
if (CV<0) {
CV <- -CV
message("Negative CV changed to ",CV,".")
}
if (targetpower>=1 | minpower>=1 | targetpower<=0 | minpower<=0 )
stop("Power values have to be within 0...1")
if (minpower>=targetpower) stop("Minimum acceptable power must be < than target")
if (Rver<"3.1.0"){
if (minpower < 0.5) stop("Minimum acceptable power must be >=0.5.")
if (targetpower < 0.5) stop("Target power must be >=0.5.")
} else {
if (minpower <= 0.5)
message("Note: Minimum acceptable power <=0.5 doesn't make much sense.")
if (targetpower <= 0.5)
message("Note: Target power <=0.5 doesn't make much sense.")
}
# check regulator, check design
regulator <- toupper(regulator)
reg <- match.arg(regulator)
design <- match.arg(design)
# to avoid recoding of Helmut's code
GMR <- theta0
# functions to be used with uniroot
# reg is visible as long as these functions are defined within pa.scABE()
pwrCV <- function(x, ...) {
if (reg=="FDA") {
power.RSABE(CV=x, ...) - minpower
} else {
power.scABEL(CV=x, ...)- minpower
}
}
pwrGMR <- function(x, ...) {
if (reg == "FDA"){
power.RSABE(theta0=x, ...) - minpower
} else {
power.scABEL(theta0=x, ...) - minpower
}
}
if (reg == "FDA") {
res <- sampleN.RSABE(CV=CV, theta0=GMR, targetpower=targetpower,
design=design, print=FALSE, details=FALSE, ...)
} else {
res <- sampleN.scABEL(CV=CV, theta0=GMR, targetpower=targetpower,
design=design, regulator=reg,
print=FALSE, details=FALSE, ...)
}
n.est <- res[1, "Sample size" ]
pwr.est <- res[1, "Achieved power"]
# various minimum sample size conditions
# generally 12
# EMA: 24 if 2x2x3
# FDA: 24
incr <- FALSE
if (reg == "EMA" & design == "2x2x3") {
if (n.est < 24) {
incr <- TRUE
n.est <- 24
pwr.est <- power.scABEL(CV=CV, n=n.est, theta0=GMR, design=design,
regulator="EMA", ...)
}
}
if (reg == "FDA") {
if (n.est < 24) {
incr <- TRUE
n.est <- 24
pwr.est <- power.RSABE(CV=CV, n=n.est, theta0=GMR, design=design, ...)
}
} else {
if (n.est < 12) {
incr <- TRUE
n.ext <- 12
pwr.est <- power.scABEL(CV=CV, n=n.est, theta0=GMR, design=design,
regulator=reg, ...)
}
}
res[1, "Sample size"] <- n.est
res[1, "Achieved power"] <- pwr.est
########################################
# max. CV for minimum acceptable power #
########################################
# TODO: test many cases. the form of the curves suggest that uniroot may fail!
if (Rver<"3.1.0"){
CV.max <- uniroot(pwrCV, c(CV, 30*CV), tol=1e-7,
n=n.est, design=design, theta0=GMR, regulator=reg, ...)$root
} else {
CV.max <- uniroot(pwrCV, c(CV, 30*CV), tol=1e-7, extendInt ="downX",
n=n.est, design=design, theta0=GMR, regulator=reg, ...)$root
}
# points for plotting
seg <- 75;
# what did he do here? D.L.
# amateur! [Helmuts comment for himself :-))]
CVs <- c(seq(CV*2/3, CV, length.out=seg*.25), seq(CV, CV.max, length.out=seg*.75))
CVs <- c(head(CVs, floor(seg*0.25)), tail(CVs, seg*0.75))
pBECV <- vector("numeric", length=length(CVs))
for (j in seq_along(CVs)) {
if (reg == "FDA") {
pBECV[j] <- power.RSABE(CV=CVs[j], n=n.est, theta0=GMR, design=design, ...)
# } else if (reg == "HC"){
# pBECV[j] <- power.scABEL2(CV=CVs[j], n=n.est, theta0=GMR, design=design,
# regulator=reg, ...)
} else {
pBECV[j] <- power.scABEL(CV=CVs[j], n=n.est, theta0=GMR, design=design,
regulator=reg, ...)
}
}
##################################################################
# min. GMR for minimum accept. power, may also be max if GMR>1! #
##################################################################
# if (reg == "EMA") {
# # scABEL (rounded switch acc. to BE-GL and Q&A document)
# ifelse(CV <= 0.5, UL <- exp(0.76*CV2se(CV)), UL <- exp(0.76*CV2se(0.5)))
# if (CV <= 0.3) UL <- 1.25
# } else {
# # RSABE (exact switching cond.: 0.8925742...)
# ifelse(CV > 0.3, UL <- exp(log(1.25)/0.25*CV2se(CV)), UL <- 1.25)
# }
UL <- scABEL(CV, regulator=reg)[2]
ifelse(GMR <= 1, interval <- c(1/UL, 1), interval <- c(GMR, UL)) # scaled borders!
# extending interval for extrem cases
updown <- ifelse(GMR <= 1, "upX", "downX")
if (Rver<"3.1.0"){
GMR.min <- uniroot(pwrGMR, interval, tol=1e-7,
n=n.est, CV=CV, design=design, regulator=reg, ...)$root
} else {
GMR.min <- uniroot(pwrGMR, interval, tol=1e-7, extendInt=updown,
n=n.est, CV=CV, design=design, regulator=reg, ...)$root
}
seg <- 50 # save some speed compared to Helmuts always 75
GMRs <- seq(GMR.min, GMR, length.out=seg)
pBEGMR <- vector("numeric", length=length(GMRs))
for (j in seq_along(GMRs)) {
if (reg == "FDA") {
pBEGMR[j] <- power.RSABE(CV=CV, n=n.est, theta0=GMRs[j],
design=design, ...)
} else {
pBEGMR[j] <- power.scABEL(CV=CV, n=n.est, theta0=GMRs[j],
design=design, regulator=reg, ...)
}
}
####################################
# min. n for minimum accept. power #
# workaround, since uniroot() does #
# not accept two vectors as limits #
####################################
#Ns <- seq(n.est, 12, by=-1) # don't drop below 12 subjects (original)
Ns <- seq(n.est, 12)
if (n.est==12) Ns <- seq(n.est, 6)
pwrN <- pwr.est
n.min <- NULL; pBEn <- NULL
#######################################################
# THX to Detlew Labes for this part of the code! #
# See: http://forum.bebac.at/forum_entry.php?id=13375 #
#######################################################
# get # of sequences
# dont need this here since unbalancedness is handled
# by power.scABEL() or power.RSABE() itself
# seqs <- designs[designs$design==design,"steps"]
# n <- vector("numeric", length=seqs)
# ni <- 1:seqs
# may it be that j grows greater than length(Ns)? paranoia
nNs <- length(Ns)
j <- 0
while(pwrN >= minpower & j<nNs){
j <- j+1
# distribute total Ns to the sequence groups
# n[-seqs] is n[1:(seqs-1)]
# don't need this here, is handled in the power functions
# n[-seqs] <- diff(floor(Ns[j]*ni/seqs))
# n[seqs] <- Ns[j] -sum(n[-seqs])
# suppress messages regarding unbalanced designs
suppressMessages(
if (reg == "FDA") {
pwrN <- power.RSABE(CV=CV, n=Ns[j], theta0=GMR, design=design, ...)
} else {
pwrN <- power.scABEL(CV=CV, n=Ns[j], theta0=GMR, design=design,
regulator=reg, ...)
}
)
if (pwrN >= minpower) {
n.min <- c(n.min, Ns[j])
pBEn <- c(pBEn, pwrN)
} else {
break
}
}
# plots are contained in the S3 method plot
# return what?
ret <-list(plan=res,
paCV=data.frame(CV=CVs, pwr=pBECV),
paGMR=data.frame(theta0=GMRs, pwr=pBEGMR),
paN=data.frame(N=n.min, pwr=pBEn),
minpower=minpower,
method="scABE", regulator=reg,
incr=incr
)
class(ret) <- c("pwrA", class(ret))
return(ret)
}
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.