# ----------------------------------------------------------------------------
# simulate replicate design subject data and evaluate via EMA ABEL method
# ANOVA & average BE with expanding limits
#
# Author D. Labes with suggestions by Ben
# This is the version trying to optimize more things
# using Ben's function to simulate multiple right-hand sides
# Trying to made more self-contained, df and C2 from a first, single lm() call
# design definition now alternatively via data.frame 'design_dta'
# a data.frame with columns subject, sequence, period and tmt and optional logval
# sequence and tmt in R/T notation
# ----------------------------------------------------------------------------
power.scABEL.sdsims <- function(alpha=0.05, theta1, theta2, theta0, CV, n,
design=c("2x3x3", "2x2x4", "2x2x3"),
design_dta=NULL, regulator,
nsims=1E5, details=FALSE, setseed=TRUE, progress)
{
# check theta0 and ABE limits
if (missing(theta0)) theta0 <- 0.90
if (length(theta0)>1) {
theta0 <- theta0[2]
warning(paste0("theta0 has to be scalar. theta0 = ",
theta0, " used."), call. = FALSE)
}
if (missing(theta1) & missing(theta2)) theta1 <- 0.8
if (missing(theta2)) theta2 <- 1/theta1
if (missing(theta1)) theta1 <- 1/theta2
# check regulator
if (missing(regulator)) regulator <- "EMA"
reg <- reg_check(regulator)
if (reg$est_method=="ISC") stop("ISC evaluation not allowed in this function.")
# Check CV
if (missing(CV)) stop("CV must be given!")
# CV scalar or vector
CVwT <- CV[1]
if (length(CV)==2) CVwR <- CV[2] else CVwR <- CVwT
# intra-subject variabilities from CV
s2wT <- CV2mse(CVwT)
s2wR <- CV2mse(CVwR)
if (is.null(design_dta)){
# check design
desi <- match.arg(design)
# degrees of freedom no longer set here
# will be determined in the working horse via a first, single call of lm
if(desi=="2x3x3"){
bk <- 1.5
bkni <- 1/6
seqs <- c("TRR", "RTR", "RRT")
}
if(desi=="2x2x3"){
bk <- 1.5
bkni <- 3/8
seqs <- c("TRT", "RTR")
}
if(desi=="2x2x4"){
bk <- 1
bkni <- 1/4
seqs <- c("TRTR", "RTRT")
}
seqn <- length(seqs)
# check n
if (missing(n)) stop("Number of subjects n must be given!")
# check n: vector or scalar
if (length(n)==1){
# for unbalanced designs we divide the ns by ourself
# to have 'smallest' imbalance
nv <- nvec(n=n, grps=seqn)
if (nv[1]!=nv[length(nv)]){
message("Unbalanced design. n(i)=", paste(nv, collapse="/"), " assumed.")
}
C2 <- sum(1/nv)*bkni
n <- sum(nv)
} else {
# check length
if (length(n)!=seqn) stop("n must be a vector of length=",seqn,"!")
C2 <- sum(1/n)*bkni
nv <- n
n <- sum(n)
}
design_dta <- prep_data2(seqs, nseq=nv, muR=log(10), ldiff=log(theta0),
s2wT=s2wT, s2wR=s2wR)
} else {
# check data.frame design_dta
# TODO
# delete NA in logval if logval is defined in data.frame
if ("logval" %in% names(design_dta)){
design_dta <- design_dta[!is.na(design_dta$logval),]
}
seqs <- unique(design_dta$sequence)
seqn <- length(seqs)
nv <- as.numeric(table(unique(design_dta[,c("subject", "sequence")])[,"sequence"]))
n <- sum(nv)
# C2 is calculated in the working horse
# above for definition via design and n it is given to the working horse
# for comparative purposes only
C2 <- NULL
}
# progressbar or not
if(missing(progress)) {
progress <- FALSE
if(nsims>=5E5) progress <- TRUE
if(nsims>=1E5 & n>72) progress <- TRUE
}
# after introducing multiple right-hand sides .lm.fit seems to have no longer
# an advantage in run-time, fitmethod removed May 2018
# call the working horse
pwr <- .pwr.ABEL.sds(design_dta=design_dta, ldiff=log(theta0),
s2WR=s2wR, s2WT=s2wT, C2=C2, nsims=nsims, regulator=reg,
ln_lBEL=log(theta1), ln_uBEL=log(theta2), alpha=alpha,
setseed=setseed, details=details, progress=progress)
pwr
}
# alias
power.scABEL.sds <- power.scABEL.sdsims
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.