# -------------------------------------------------------------------------
# sample size function without all the overhead
# for 2x2 crossover or 2-group parallel
# power via nct or shifted t approximation or exact
#
# Version with vectorize sample size w.r.t mse and/or ltheta0
#
# author D. Labes Jan 2015
# Aug 2017: df changed to n-3 via new argument dfc (df as character string)
# but only for Potvin like TSD's
# -------------------------------------------------------------------------
# is also used for parallel groups with bk=4
.sampleN2 <- function(alpha=0.05, targetpower=0.8, ltheta0, ltheta1=log(0.8),
ltheta2=log(1.25), mse, method="nct", bk=2, dfc="n-3")
{
# se and ltheta0/diffm must have the same length to vectorize propperly!
if (length(mse)==1) mse <- rep(mse, times=length(ltheta0))
if (length(ltheta0)==1) ltheta0 <- rep(ltheta0, times=length(mse))
# return 'Inf' if ltheta0 not between or very near to ltheta1, ltheta2
ns <- ifelse((ltheta0-ltheta1)<1.25e-5 | (ltheta2-ltheta0)<1.25e-5, Inf, 0)
# design characteristics for 2-group parallel and 2x2 crossover design
steps <- 2 # stepsize for sample size search (balanced designs)
nmin <- 4 # minimum n
se <- sqrt(mse[is.finite(ns)])
diffm <- ltheta0[is.finite(ns)]
# start value from large sample approx. (hidden func.)
# Jan 2015 changed to modified Zhang's formula
# gives at least for 2x2 the best estimate (max diff to n: +-4)
.sampleN0_3 <- NULL
n <- .sampleN0_3(alpha, targetpower, ltheta1, ltheta2, diffm, se, steps, bk)
n <- ifelse(n<nmin, nmin, n)
# method=="ls" is not used yet, in power.2stage.ssr() the 'original' ls approx
# is used exactly as given in the paper of Golkowski et al.
if(method=="ls") return(n)
# degrees of freedom as expression
# usual n-2 for 2x2 crossover and 2-group parallel design
# dfe <- parse(text="n-2", srcfile=NULL)
# or should that read n-3? see Kieser/Rauch for inclusion of a stage term
# the jury has decided so
dfe <- parse(text=dfc, srcfile=NULL)
df <- eval(dfe)
pow <- .calc.power(alpha, ltheta1, ltheta2, diffm, sem=se*sqrt(bk/n), df, method)
iter <- rep(0, times=length(se))
imax <- 50
# imax =max. iter >50 is emergency brake
# this is eventually not necessary, depends on quality of sampleN0
# in experimentation I have seen max. of 2-3 steps
# reformulation with only one loop does not shorten the code considerable
# --- loop until power <= target power, step-down
index <- pow>targetpower & n>nmin
while (any(index)) {
n[index] <- n[index]-steps # step down if start power is to high
iter[index] <- iter[index]+1
df <- eval(dfe)
pow[index] <- .calc.power(alpha, ltheta1, ltheta2, diffm[index],
sem=se[index]*sqrt(bk/n[index]), df[index],
method)
index <- pow>targetpower & n>nmin & iter<=imax
# loop results in n with power too low
# must step one up again. is done in the next loop
}
# --- loop until power >= target power
index <- pow<targetpower
while (any(index)) {
n[index] <- n[index]+steps
iter[index] <- iter[index]+1
df <- eval(dfe)
pow[index] <- .calc.power(alpha, ltheta1, ltheta2, diffm[index],
sem=se[index]*sqrt(bk/n[index]), df[index],
method)
index <- pow<targetpower & iter<=imax
}
# nlast <- n
# if ((up & pow<targetpower) | (down & pow>targetpower) ) {
# n <- NA
# }
# combine the Inf and n
ns[ns==0] <- n
n <- ns
# return only n here
return(n)
} # end of function
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.