# function to check if function is called with FDA in its name
check4FDA <- function(fname)
{
if(grepl("FDA", fname))
message("Function ", fname, " is deprecated.\n",
"Use version without FDA in it's name.")
}
#---------------------------------------------------------------------------
# Simulate full replicate designs and calculate scaled ABE power
# according to FDA Warfarin guidance
#
# Author: dlabes
#---------------------------------------------------------------------------
power.NTID <- function(alpha=0.05, theta1, theta2, theta0, CV, n,
design=c("2x2x4", "2x2x3"), nsims=1E5, details=FALSE,
setseed=TRUE)
{
#check if function name contains FDA and warn function deprecated
check4FDA(fname=as.character(sys.call())[1])
if (missing(CV)) stop("CV must be given!", call.=FALSE)
if (missing(n)) stop("Number of subjects n must be given!", call.=FALSE)
if (missing(theta0)) theta0 <- 0.975 # tighter content limits for NTID
if (missing(theta1) & missing(theta2)) theta1 <- 0.8
if (missing(theta2)) theta2 <- 1/theta1
if (missing(theta1)) theta1 <- 1/theta2
design <- match.arg(design)
if(design=="2x2x4"){
seqs <- 2
dfe <- parse(text="n-2", srcfile=NULL)
dfRRe <- parse(text="n-2", srcfile=NULL)
dfTTe <- parse(text="n-2", srcfile=NULL)
}
if(design=="2x2x3"){
seqs <- 2
dfe <- parse(text="n-2", srcfile=NULL)
dfRRe <- parse(text="n/2-1", srcfile=NULL) # balanced only, not used here
dfTTe <- parse(text="n/2-1", srcfile=NULL) # balanced only, not used here
}
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)
# FDA constant
r_const <- -log(0.9)/0.10
if (length(n)==1){
# we assume n=ntotal
# for unbalanced designs we divide the ns by ourself
# in such a way that we have only small imbalance
nv <- nvec(n=n, grps=seqs)
if (nv[1]!=nv[length(nv)]){
message("Unbalanced design. n(i)=", paste(nv, collapse="/"), " assumed.")
}
C3 <- sum(1/nv)/seqs^2
n <- sum(nv)
} else {
# we assume n = vector of n's in sequence groups
# check length
if (length(n)!=seqs) stop("n must be a vector of length=",seqs,"!", call.=FALSE)
nv <- n
C3 <- sum(1/n)/seqs^2
n <- sum(n)
}
if(design=="2x2x4"){
dfRR <- eval(dfRRe)
dfTT <- dfRR
# expectation of mse of the ANOVA of intra-subject contrasts T-R
Emse <- (s2wT + s2wR)/2
}
if(design=="2x2x3"){
dfTT <- nv[1]-1
dfRR <- nv[2]-1
w1 <- dfRR/(dfRR+dfTT); w2 <- dfTT/(dfRR+dfTT)
# expectation of mse of the ANOVA of intra-subject contrasts T-R
# always via unbalanced formula
Emse <- w1*(s2wT+s2wR/2) + w2*(s2wT/2+s2wR)
}
# start time measurement
ptm <- proc.time()
df <- eval(dfe)
# sd of the mean T-R (point estimator)
sdm <- sqrt(Emse*C3)
mlog <- log(theta0)
if(setseed) set.seed(123456)
p <- .power.NTID(mlog, sdm, C3, Emse, df, s2wR, dfRR, s2wT, dfTT, nsims,
r_const, ln_lBEL=log(theta1),ln_uBEL=log(theta2), alpha=alpha)
if (details) {
ptm <- summary(proc.time()-ptm)
message(nsims,"sims. Time elapsed (sec): ",
formatC(ptm["elapsed"], digits=2), "\n")
#print(ptm)
# return power and components
names(p) <- c("p(BE)", "p(BE-sABEc)", "p(BE-ABE)", "p(BE-sratio)")
p
} else {
# return only the 'power'
as.numeric(p["BE"])
}
}
# working horse of RSABE for NTID's
.power.NTID <- function(mlog, sdm, C3, Emse, df, s2wR, dfRR, s2wT, dfTT, nsims,
r_const=-log(0.9)/0.1, ln_lBEL=log(0.8), ln_uBEL=log(1.25),
alpha=0.05)
{
tval <- qt(1-alpha,df)
chisqval <- qchisq(1-alpha, dfRR)
r2const <- r_const^2
Fval <- qf(1-alpha, dfTT, dfRR, lower.tail=FALSE)
counts <- rep.int(0, times=4)
names(counts) <- c("BE", "BEul", "BEabe", "BEsratio")
# to avoid memory problems for high number of sims
# we are working with chunks of 1e7
chunks <- 1
nsi <- nsims
if (nsims>1E7) {
chunks <- ceiling(nsims/1E7)
nsi <- 1E7
}
for (iter in 1:chunks) {
# if chunks*1E7 >nsims correct nsi to given nsims
if(iter==chunks) nsi <- nsims-(chunks-1)*nsi
# simulate sample mean via its normal distribution
means <- rnorm(nsi, mean=mlog, sd=sdm)
# simulate sample sd2s via chi-square distri
sd2s <- Emse*C3*rchisq(nsi, df)/df
# simulate sample value s2wRs via chi-square distri
s2wRs <- s2wR*rchisq(nsi, dfRR)/dfRR
# simulate sample value s2wTs via chi-square distri
s2wTs <- s2wT*rchisq(nsi, dfTT)/dfTT
SEs <- sqrt(sd2s)
# conventional (1-2*alpha) CI's for T-R
hw <- tval*SEs
lCL <- means - hw
uCL <- means + hw
# conventional ABE
BEABE <- ((ln_lBEL<=lCL) & (uCL<=ln_uBEL))
# upper 95% CI linearized SABE criterion
# with -SEs^2 the 'unknown' x from the warfarin guidance
Em <- means^2 - SEs^2
Es <- r2const*s2wRs
#Cm <- (abs(means) + hw)^2
Cm <- ifelse(abs(lCL)>abs(uCL),abs(lCL)^2,abs(uCL)^2)
Cs <- Es*dfRR/chisqval
SABEc95 <- Em - Es + sqrt((Cm-Em)^2 + (Cs-Es)^2)
# what if no scaling has to be applied?
BEscABE <- (SABEc95 <= 0)
# save memory
rm(SEs, hw, Em, Es, Cm, Cs)
# upper limit of ratio swT/swR
ul_sratio <- sqrt(s2wTs/s2wRs/Fval)
# upper limit <= 2.5?
BEsratio <- (ul_sratio <= 2.5)
counts["BEabe"] <- counts["BEabe"] + sum(BEABE)
counts["BEul"] <- counts["BEul"] + sum(BEscABE)
counts["BEsratio"] <- counts["BEsratio"] + sum(BEsratio)
# test without s-ratio
#counts["BE"] <- counts["BE"] + sum(BEscABE & BEABE)
counts["BE"] <- counts["BE"] + sum(BEscABE & BEABE & BEsratio)
} # end over chunks
# return the pBEs
counts/nsims
}
# alias 'power.NTIDFDA' removed since this evaluation is not only requested
# by FDA but also by the China CDE
# power.NTIDFDA <- power.NTID
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.