Nothing
## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#"
)
## ----setup--------------------------------------------------------------------
library(PowerTOST) # attach the library
## ----example1a----------------------------------------------------------------
sampleN.TOST(CV = 0.30)
## ----example1b----------------------------------------------------------------
sampleN.TOST(CV = 0.30, details = FALSE, print = FALSE)[["Sample size"]]
## ----example1c----------------------------------------------------------------
power.TOST(CV = 0.30, n = 39)
## ----example1vect-------------------------------------------------------------
sampleN.TOST.vectorized <- function(CVs, theta0s, ...) {
n <- power <- matrix(ncol = length(CVs), nrow = length(theta0s))
for (i in seq_along(theta0s)) {
for (j in seq_along(CVs)) {
tmp <- sampleN.TOST(CV = CVs[j], theta0 = theta0s[i], ...)
n[i, j] <- tmp[["Sample size"]]
power[i, j] <- tmp[["Achieved power"]]
}
}
DecPlaces <- function(x) match(TRUE, round(x, 1:15) == x)
fmt.col <- paste0("CV %.", max(sapply(CVs, FUN = DecPlaces),
na.rm = TRUE), "f")
fmt.row <- paste0("theta %.", max(sapply(theta0s, FUN = DecPlaces),
na.rm = TRUE), "f")
colnames(power) <- colnames(n) <- sprintf(fmt.col, CVs)
rownames(power) <- rownames(n) <- sprintf(fmt.row, theta0s)
res <- list(n = n, power = power)
return(res)
}
CVs <- seq(0.20, 0.40, 0.05)
theta0s <- seq(0.90, 0.95, 0.01)
x <- sampleN.TOST.vectorized(CV = CVs, theta0 = theta0s,
details = FALSE, print = FALSE)
cat("Sample size\n"); print(x$n); cat("Achieved power\n"); print(signif(x$power, digits = 5))
## ----example1d----------------------------------------------------------------
designs <- c("2x2x2", "2x2x3", "2x3x3", "2x2x4")
# data.frame of results
res <- data.frame(design = designs, n = NA_integer_,
power = NA_real_, n.do = NA_integer_,
power.do = NA_real_,
stringsAsFactors = FALSE) # this line for R <4.0.0
for (i in 1:4) {
# print = FALSE suppresses output to the console
# we are only interested in columns 7-8
# let's also calculate power for one dropout
res[i, 2:3] <- sampleN.TOST(CV = 0.30, design = res$design[i],
print = FALSE)[7:8]
res[i, 4] <- res[i, 2] - 1
res[i, 5] <- suppressMessages(
power.TOST(CV = 0.30, design = res$design[i],
n = res[i, 4]))
}
print(res, row.names = FALSE)
## ----example1e, echo = FALSE--------------------------------------------------
designs <- c("2x2x2", "2x2x3", "2x3x3", "2x2x4")
res <- data.frame(design = rep(NA_character_, 4), name = NA_character_,
n = NA_integer_, formula = NA_character_,
df = NA_real_, t.value = NA_real_,
stringsAsFactors = FALSE) # this line for R <4.0.0
res[, c(2:1, 4)] <- known.designs()[which(known.designs()[, 2] %in% designs), c(9, 2:3)]
for (i in 1:4) {
res$n[i] <- sampleN.TOST(CV = 0.30, design = res$design[i],
print = FALSE)[["Sample size"]]
e <- parse(text = res[i, 4], srcfile = NULL)
n <- res$n[i]
res[i, 5] <- eval(e)
res$t.value[i] <- signif(qt(1 - 0.05, df = res[i, 5]), 4)
}
res <- res[with(res, order(-n, design)), ]
print(res, row.names = FALSE)
## ----example1f----------------------------------------------------------------
grouping <- function(capacity, n) {
# split sample size into >=2 groups based on capacity
if (n <= capacity) { # make equal groups
ngrp <- rep(ceiling(n / 2), 2)
} else { # at least one = capacity
ngrp <- rep(0, ceiling(n / capacity))
grps <- length(ngrp)
ngrp[1] <- capacity
for (j in 2:grps) {
n.tot <- sum(ngrp) # what we have so far
if (n.tot + capacity <= n) {
ngrp[j] <- capacity
} else {
ngrp[j] <- n - n.tot
}
}
}
return(ngrp = list(grps = length(ngrp), ngrp = ngrp))
}
CV <- 0.30
capacity <- 24 # clinical capacity
res <- data.frame(n = NA_integer_, grps = NA_integer_,
n.grp = NA_integer_,
m.1 = NA_real_, m.2 = NA_real_)
x <- sampleN.TOST(CV = CV, print = FALSE, details = FALSE)
res$n <- x[["Sample size"]]
res$m.1 <- x[["Achieved power"]]
x <- grouping(capacity = capacity, n = res$n)
res$grps <- x[["grps"]]
ngrp <- x[["ngrp"]]
res$n.grp <- paste(ngrp, collapse = "|")
res$m.2 <- power.TOST.sds(CV = CV, n = res$n, grps = res$grps,
ngrp = ngrp, gmodel = 2, progress = FALSE)
res$loss <- 100*(res$m.2 - res$m.1) / res$m.1
names(res)[2:6] <- c("groups", "n/group", "pooled model",
"group model", "loss (%)")
res[1, 4:6] <- sprintf("%6.4f", res[1, 4:6])
cat("Achieved power and relative loss\n"); print(res, row.names = FALSE)
## ----example2-----------------------------------------------------------------
sampleN.RatioF(CV = 0.20, CVb = 0.40)
## ----example3a----------------------------------------------------------------
sampleN.TOST(CV = 0.20, theta0 = 0.92)
## ----example3b----------------------------------------------------------------
df <- 16 - 2 # degrees of freedom of the 2x2x2 crossover pilot
CVCL(CV = 0.20, df = df, side = "upper", alpha = 0.20)[["upper CL"]]
## ----example3c----------------------------------------------------------------
CL.upper <- CVCL(CV = 0.20, df = 16 - 2, side = "upper",
alpha = 0.20)[["upper CL"]]
res <- sampleN.TOST(CV = CL.upper, theta0 = 0.92, print = FALSE)
print(res[7:8], row.names = FALSE)
## ----example3d----------------------------------------------------------------
CL.upper <- CVCL(CV = 0.20, df = 16 - 2, side = "upper",
alpha = 0.20)[["upper CL"]]
power.TOST(CV = CL.upper, theta0 = 0.92, n = 28)
## ----example3e----------------------------------------------------------------
power.TOST(CV = 0.22, theta0 = 0.90, n = 40)
## ----example3f----------------------------------------------------------------
expsampleN.TOST(CV = 0.20, theta0 = 0.92, prior.type = "CV",
prior.parm = list(m = 16, design = "2x2x2"))
## ----example3g----------------------------------------------------------------
expsampleN.TOST(CV = 0.20, theta0 = 0.92, prior.type = "theta0",
prior.parm = list(m = 16, design = "2x2x2"))
## ----example3h----------------------------------------------------------------
expsampleN.TOST(CV = 0.20, theta0 = 0.92, prior.type = "both",
prior.parm = list(m = 16, design = "2x2x2"),
details = FALSE)
## ----example4-----------------------------------------------------------------
CV <- 0.214
res <- data.frame(target = c(rep(0.8, 5), rep(0.9, 5)),
theta0 = rep(c(1, seq(0.95, 0.92, -0.01)), 2),
n.1 = NA_integer_, power = NA_real_,
sigma.u = rep(c(0.0005, seq(0.05, 0.08, 0.01)), 2),
n.2 = NA_integer_, assurance = NA_real_)
for (i in 1:nrow(res)) {
res[i, 3:4] <- sampleN.TOST(CV = CV, targetpower = res$target[i],
theta0 = res$theta0[i], print = FALSE)[7:8]
res[i, 6:7] <- expsampleN.TOST(CV = CV, targetpower = res$target[i],
theta0 = 1, # mandatory!
prior.type = "theta0",
prior.parm = list(sem = res$sigma.u[i]),
print = FALSE)[9:10]
}
res <- signif(res, 3)
res[, 5] <- sprintf("%.2f", res[, 5])
names(res)[c(3, 6)] <- "n"
print(res, row.names = FALSE)
## ----example4a----------------------------------------------------------------
planned <- "2x2x2"
logscale <- FALSE
theta0 <- -5
theta1 <- -15
theta2 <- +15 # if not given, -theta1 is used
SD.resid <- 20 # residual standard deviation
sampleN.TOST(CV = SD.resid, theta0 = theta0, theta1 = theta1,
theta2 = theta2, logscale = logscale, design = planned)
## ----example4b----------------------------------------------------------------
known <- known.designs()[, c(2, 6)] # extract relevant information
bk <- known[known$design == planned, "bk"] # retrieve design constant
txt <- paste0("The design constant for design \"",
planned, "\" is ", bk)
SD.delta <- 28 # standard deviation of the difference
SD.resid <- SD.delta / sqrt(bk) # convert to residual SD
cat(txt); sampleN.TOST(CV = SD.resid, theta0 = theta0,
theta1 = theta1, theta2 = theta2,
logscale = logscale, design = planned)
## ----higher_order1------------------------------------------------------------
CV <- 0.20
res <- data.frame(design = c("3x6x3", "2x2x2"), n = NA_integer_,
power = NA_real_, stringsAsFactors = FALSE)
for (i in 1:2) {
res[i, 2:3] <- sampleN.TOST(CV = CV, design = res$design[i],
print = FALSE)[7:8]
}
print(res, row.names = FALSE)
## ----higher_order2------------------------------------------------------------
CV <- 0.20
res <- data.frame(design = c("4x4", "2x2x2"), n = NA_integer_,
power = NA_real_, stringsAsFactors = FALSE)
for (i in 1:2) {
res[i, 2:3] <- sampleN.TOST(CV = CV, design = res$design[i],
print = FALSE)[7:8]
}
print(res, row.names = FALSE)
## ----power1-------------------------------------------------------------------
n <- c(14 - 1, 14 - 2) # 14 dosed in each sequence
round(100*CI.BE(pe = 0.90, CV = 0.25, n = n), 2)
## ----power2-------------------------------------------------------------------
power.TOST(CV = 0.25, theta0 = 0.90, n = c(13, 12)) # observed values
## ----power3-------------------------------------------------------------------
power.TOST(CV = 0.20, theta0 = 0.92, n = 28) # assumed in planning
power.TOST(CV = 0.25, theta0 = 1, n = c(13, 12)) # observed CV but wrong T/R-ratio
## ----pooling------------------------------------------------------------------
CVs <- (" CV | n | design | study
0.20 | 16 | 2x2x2 | pilot
0.25 | 25 | 2x2x2 | pivotal")
txtcon <- textConnection(CVs)
data <- read.table(txtcon, header = TRUE, sep = "|",
strip.white = TRUE, as.is = TRUE)
close(txtcon)
print(CVpooled(data, alpha = 0.20), digits = 4, verbose = TRUE)
## ----PA-----------------------------------------------------------------------
pa.ABE(CV = 0.20, theta0 = 0.92)
## ----DO-----------------------------------------------------------------------
balance <- function(x, y) {
return(y * (x %/% y + as.logical(x %% y)))
}
do <- 0.15 # anticipated dropout-rate 15%
seqs <- 3
n <- seq(12L, 96L, 12L)
res <- data.frame(n = n,
adj1 = balance(n / (1 - do), seqs), # correct
elig1 = NA_integer_, diff1 = NA_integer_,
adj2 = balance(n * (1 + do), seqs), # wrong
elig2 = NA_integer_, diff2 = NA_integer_)
res$elig1 <- floor(res$adj1 * (1 - do))
res$diff1 <- sprintf("%+i", res$elig1 - n)
res$elig2 <- floor(res$adj2 * (1 - do))
res$diff2 <- sprintf("%+i", res$elig2 - n)
invisible(
ifelse(res$elig1 - n >=0,
res$optim <- res$elig1,
res$optim <- res$elig2))
res$diff <- sprintf("%+i", res$optim - n)
names(res)[c(2, 5)] <- c("n'1", "n'2")
res$diff1[which(res$diff1 == "+0")] <- "\u00B10"
res$diff2[which(res$diff2 == "+0")] <- "\u00B10"
res$diff[which(res$diff == "+0")] <- "\u00B10"
print(res, row.names = FALSE)
## ----lit1---------------------------------------------------------------------
CVfromCI(lower = 0.8323, upper = 1.0392, design = "2x2x4", n = 26)
## ----lit2---------------------------------------------------------------------
n <- 26
CV.est <- CVfromCI(lower = 0.8323, upper = 1.0392,
design = "2x2x4", n = 26)
n.est <- sampleN.TOST(CV = CV.est, design = "2x2x4",
print = FALSE)[["Sample size"]]
n1 <- balance(seq(n, 18, -1), 2) / 2
n2 <- n - n1
nseqs <- unique(data.frame(n1 = n1, n2 = n2, n = n))
res <- data.frame(n1 = nseqs$n1, n2 = nseqs$n2,
CV.true = NA_real_,
CV.est = CV.est, n.true = NA_integer_,
n.est = n.est)
for (i in 1:nrow(res)) {
res$CV.true[i] <- CVfromCI(lower = 0.8323, upper = 1.0392,
design = "2x2x4",
n = c(res$n1[i], res$n2[i]))
res$n.true[i] <- sampleN.TOST(CV = res$CV.true[i], design = "2x2x4",
print = FALSE)[["Sample size"]]
res$n.est[i] <- sampleN.TOST(CV = CV.est, design = "2x2x4",
print = FALSE)[["Sample size"]]
}
print(signif(res, 5), row.names = FALSE)
## ----dev----------------------------------------------------------------------
CV <- 0.21
d <- 0.05 # delta 5%, direction unknown
n <- sampleN.TOST(CV = CV, theta0 = 1 - d, print = FALSE,
details = FALSE)[["Sample size"]]
res1 <- data.frame(CV = CV, theta0 = c(1 - d, 1 / (1 - d)),
n = n, power = NA_real_)
for (i in 1:nrow(res1)) {
res1$power[i] <- power.TOST(CV = CV, theta0 = res1$theta0[i], n = n)
}
n <- sampleN.TOST(CV = CV, theta0 = 1 + d,
print = FALSE)[["Sample size"]]
res2 <- data.frame(CV = CV, theta0 = c(1 + d, 1 / (1 + d)),
n = n, power = NA_real_)
for (i in 1:nrow(res1)) {
res2$power[i] <- power.TOST(CV = CV, theta0 = res2$theta0[i], n = n)
}
res <- rbind(res1, res2)
print(signif(res[order(res$n, res$theta0), ], 4), row.names = FALSE)
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.