Nothing
##################################################################################
# Binomial functionality for the gsDesign package
#
# Exported Functions:
#
# ciBinomial
# nBinomial
# simBinomial
# testBinomial
#
# Hidden Functions:
#
# bpdiff
#
# Author(s): Keaven Anderson, PhD.
#
# Reviewer(s): REvolution Computing 19DEC2008 v.2.0 - William Constantine, Ph.D., Kellie Wills
#
# R Version: 2.7.2
#
##################################################################################
###
# Exported Functions
###
"ciBinomial" <- function(x1, x2, n1, n2, alpha=.05, adj=0, scale="Difference")
{
# check input arguments
checkScalar(n1, "integer", c(1, Inf))
checkScalar(n2, "integer", c(1, Inf))
checkScalar(x1, "integer", c(0, n1))
checkScalar(x2, "integer", c(0, n2))
checkScalar(alpha, "numeric", c(0, 1), c(FALSE, FALSE))
checkScalar(adj, "integer", c(0, 1))
checkScalar(scale, "character")
scale <- match.arg(tolower(scale), c("difference", "rr", "or"))
checkLengths(x1, x2)
if (scale == "difference")
{
delta <- x1 / n1 - x2 / n2
if (delta == -1)
{
lower <- -1
}
else if (testBinomial(delta0=-.9999, x1=x1, x2=x2, n1=n1, n2=n2, adj=adj)
< qnorm(alpha / 2))
{
lower <- -1
}
else
{
lower <- uniroot(bpdiff, interval=c(-.9999, delta), x1=x1, x2=x2,
n1=n1, n2=n2, adj=adj, alpha=alpha)$root
}
if (delta == 1)
{
upper <- 1
}
else if (testBinomial(delta0=.9999, x1=x1, x2=x2, n1=n1, n2=n2, adj=adj)
> -qnorm(alpha/2))
{
upper<-1
}
else
{
upper <- uniroot(bpdiff, interval=c(delta,.9999), lower.tail=TRUE,
x1=x1, x2=x2, n1=n1, n2=n2, adj=adj,
alpha=alpha)$root
}
}
else if (scale=="rr")
{
if (x1 == 0) lower <- 0
else
{
lower <- uniroot(bpdiff, interval=c(-20, 20), x1=x1, x2=x2,
n1=n1, n2=n2, adj=adj, scale="RR", alpha=alpha)$root
lower <- exp(lower)
}
if (x2 == 0) upper <- Inf
else
{
upper <- uniroot(bpdiff, interval=c(-20, 20), lower.tail=TRUE,
x1=x1, x2=x2, n1=n1, n2=n2, adj=adj, scale="RR",
alpha=alpha)$root
upper <- exp(upper)
}
}
else
{
if (x1 == 0 || x2 == n2) lower <- -Inf
else
{
lower <- uniroot(bpdiff, interval=c(-10, 10), x1=x1, x2=x2,
n1=n1, n2=n2, adj=adj, scale=scale, alpha=alpha)$root
lower <- exp(lower)
}
if (x2 == 0 || x1 == n1) upper <- Inf
else
{
upper <- uniroot(bpdiff, interval=c(-10, 10), lower.tail=TRUE,
x1=x1, x2=x2, n1=n1, n2=n2, adj=adj, scale=scale,
alpha=alpha)$root
upper <- exp(upper)
}
}
data.frame(lower=lower,upper=upper)
}
"nBinomial" <- function(p1, p2, alpha = 0.025, beta = 0.1, delta0 = 0, ratio = 1,
sided = 1, outtype = 1, scale = "Difference", n = NULL)
{
checkVector(p1, "numeric", c(0, 1), c(FALSE, FALSE))
checkVector(p2, "numeric", c(0, 1), c(FALSE, FALSE))
checkScalar(sided, "integer", c(1, 2))
checkScalar(alpha, "numeric", c(0, 1/sided), c(FALSE, FALSE))
checkVector(beta, "numeric", c(0, 1 - alpha/sided), c(FALSE,
FALSE))
checkVector(delta0, "numeric")
checkVector(ratio, "numeric", c(0, Inf), c(FALSE, FALSE))
checkScalar(outtype, "integer", c(1, 3))
checkScalar(scale, "character")
if (!is.null(n))
checkVector(n, "numeric")
scale <- match.arg(tolower(scale), c("difference", "rr",
"or", "lnor"))
if (is.null(n))
checkLengths(p1, p2, beta, delta0, ratio, allowSingle = TRUE)
else checkLengths(n, p1, p2, delta0, ratio, allowSingle = TRUE)
len <- max(sapply(list(p1, p2, beta, delta0, ratio), length))
if (len > 1) {
if (length(p1) == 1)
p1 <- array(p1, len)
if (length(p2) == 1)
p2 <- array(p2, len)
if (length(alpha) == 1)
alpha <- array(alpha, len)
if (length(beta) == 1)
beta <- array(beta, len)
if (length(delta0) == 1)
delta0 <- array(delta0, len)
if (length(ratio) == 1)
ratio <- array(ratio, len)
}
if (max(delta0 == 0) > 0 && max(p1[delta0 == 0] == p2[delta0 ==
0]) > 0) {
stop("p1 may not equal p2 when delta0 is zero")
}
z.beta <- qnorm(1 - beta)
sided[sided != 2] <- 1
z.alpha <- qnorm(1 - alpha/sided)
d0 <- (delta0 == 0)
if (scale == "difference") {
if (min(abs(p1 - p2 - delta0)) < 1e-11) {
stop("p1 - p2 may not equal delta0 when scale is \"Difference\"")
}
a <- 1 + ratio
b <- -(a + p1 + ratio * p2 + delta0 * (ratio + 2))
c <- delta0^2 + delta0 * (2 * p1 + a) + p1 + ratio *
p2
d <- -p1 * delta0 * (1 + delta0)
v <- (b/(3 * a))^3 - b * c/6/a^2 + d/2/a
u <- (sign(v) + (v == 0)) * sqrt((b/3/a)^2 - c/3/a)
w <- (pi + acos(v/u^3))/3
p10 <- 2 * u * cos(w) - b/3/a
p20 <- p10 - delta0
p10[d0] <- (p1[d0] + ratio[d0] * p2[d0])/(1 + ratio[d0])
p20[d0] <- p10[d0]
sigma0 <- sqrt((p10 * (1 - p10) + p20 * (1 - p20)/ratio) *
(ratio + 1))
sigma1 <- sqrt((p1 * (1 - p1) + p2 * (1 - p2)/ratio) *
(ratio + 1))
if (is.null(n)) {
n <- ((z.alpha * sigma0 + z.beta * sigma1)/(p1 - p2 - delta0))^2
if (outtype == 2) {
return(data.frame(cbind(n1 = n/(ratio + 1), n2 = ratio *
n/(ratio + 1))))
}
else if (outtype == 3) {
return(data.frame(cbind(n = n, n1 = n/(ratio + 1),
n2 = ratio * n/(ratio + 1), alpha = alpha,
sided = sided, beta = beta, Power = 1 - beta,
sigma0 = sigma0, sigma1 = sigma1, p1 = p1,
p2 = p2, delta0 = delta0, p10 = p10, p20 = p20)))
}
else return(n = n)
}
else {
pwr <- pnorm(-(qnorm(1 - alpha/sided) - sqrt(n) *
((p1 - p2 - delta0)/sigma0)) * sigma0/sigma1)
if (outtype == 2) {
return(data.frame(cbind(n1 = n/(ratio + 1), n2 = ratio *
n/(ratio + 1), Power = pwr)))
}
else if (outtype == 3) {
return(data.frame(cbind(n = n, n1 = n/(ratio + 1),
n2 = ratio * n/(ratio + 1), alpha = alpha,
sided = sided, beta = 1 - pwr, Power = pwr,
sigma0 = sigma0, sigma1 = sigma1, p1 = p1,
p2 = p2, delta0 = delta0, p10 = p10, p20 = p20)))
}
else return(Power = pwr)
}
}
else if (scale == "rr") {
RR <- exp(delta0)
if (min(abs(p1/p2 - RR)) < 1e-07) {
stop("p1/p2 may not equal exp(delta0) when scale=\"RR\"")
}
a <- (1 + ratio)
b <- -(RR * (1 + ratio * p2) + ratio + p1)
c <- RR * (p1 + ratio * p2)
p10 <- (-b - sqrt(b^2 - 4 * a * c))/2/a
p20 <- p10/RR
p10[d0] <- (p1[d0] + ratio[d0] * p2[d0])/(1 + ratio[d0])
p20[d0] <- p10[d0]
sigma0 <- sqrt((ratio + 1) * (p10 * (1 - p10) + RR^2 *
p20 * (1 - p20)/ratio))
sigma1 <- sqrt((ratio + 1) * (p1 * (1 - p1) + RR^2 *
p2 * (1 - p2)/ratio))
if (is.null(n)) {
n <- ((z.alpha * sigma0 + z.beta * sigma1)/(p1 - p2 * RR))^2
if (outtype == 2) {
return(data.frame(cbind(n1 = n/(ratio + 1),
n2 = ratio * n/(ratio + 1))))
}
else if (outtype == 3) {
return(data.frame(cbind(n = n, n1 = n/(ratio + 1),
n2 = ratio * n/(ratio + 1), alpha = alpha,
sided = sided, beta = beta, Power = 1-beta,
sigma0 = sigma0, sigma1 = sigma1, p1 = p1,
p2 = p2, delta0 = delta0, p10 = p10, p20 = p20)))
}
else return(n = n)
}
else {
pwr <- pnorm(-(qnorm(1 - alpha/sided) - sqrt(n) *
((p1 - p2 * RR)/sigma0)) * sigma0/sigma1)
if (outtype == 2) {
return(data.frame(cbind(n1 = n/(ratio + 1), n2 = ratio *
n/(ratio + 1), Power = pwr)))
}
else if (outtype == 3) {
return(data.frame(cbind(n = n, n1 = n/(ratio + 1),
n2 = ratio * n/(ratio + 1), alpha = alpha,
sided = sided, beta = 1 - pwr, Power = pwr,
sigma0 = sigma0, sigma1 = sigma1, p1 = p1,
p2 = p2, delta0 = delta0, p10 = p10, p20 = p20)))
}
else return(Power = pwr)
}
}
else {
OR <- exp(-delta0)
if (min(abs(p1/(1 - p1)/p2 * (1 - p2) * OR) - 1) < 1e-07) {
stop("p1/(1-p1)/p2*(1-p2) may not equal exp(delta0) when scale=\"OR\"")
}
a <- OR - 1
b <- 1 + ratio * OR + (1 - OR) * (ratio * p2 + p1)
c <- -(ratio * p2 + p1)
p10 <- (-b + sqrt(b^2 - 4 * a * c))/2/a
p20 <- OR * p10/(1 + p10 * (OR - 1))
p10[d0] <- (p1[d0] + ratio[d0] * p2[d0])/(1 + ratio[d0])
p20[d0] <- p10[d0]
sigma0 <- sqrt((ratio + 1) * (1/p10/(1 - p10) + 1/p20/(1 - p20)/ratio))
sigma1 <- sqrt((ratio + 1) * (1/p1/(1 - p1) + 1/p2/(1 - p2)/ratio))
if (is.null(n)) {
n <- ((z.alpha * sigma0 + z.beta * sigma1)/log(OR/p2 * (1 - p2) * p1/(1 - p1)))^2
if (outtype == 2) {
return(data.frame(cbind(n1 = n/(ratio + 1), n2 = ratio *
n/(ratio + 1))))
}
else if (outtype == 3) {
return(data.frame(cbind(n = n, n1 = n/(ratio + 1),
n2 = ratio * n/(ratio + 1), alpha = alpha,
sided = sided, beta = beta, Power = 1-beta,
sigma0 = sigma0, sigma1 = sigma1, p1 = p1,
p2 = p2, delta0 = delta0, p10 = p10, p20 = p20)))
}
else {
return(n = n)
}
}
else {
pwr <- pnorm(-(qnorm(1 - alpha/sided) - sqrt(n) *
(log(OR/p2 * (1 - p2) * p1/(1 - p1))/sigma0)) *
sigma0/sigma1)
if (outtype == 2) {
return(data.frame(cbind(n1 = n/(ratio + 1), n2 = ratio *
n/(ratio + 1), Power = pwr)))
}
else if (outtype == 3) {
return(data.frame(cbind(n = n, n1 = n/(ratio + 1),
n2 = ratio * n/(ratio + 1), alpha = alpha,
sided = sided, beta = 1 - pwr, Power = pwr,
sigma0 = sigma0, sigma1 = sigma1, p1 = p1,
p2 = p2, delta0 = delta0, p10 = p10, p20 = p20)))
}
else return(Power = pwr)
}
}
}
"simBinomial" <- function(p1, p2, n1, n2, delta0=0, nsim=10000, chisq=0, adj=0,
scale="Difference")
{
# check input arguments
# Note: delta0, chisq, adj, and scale checked by testBinomial() function
checkVector(p1, "numeric", c(0, 1), c(FALSE, FALSE))
checkVector(p2, "numeric", c(0, 1), c(FALSE, FALSE))
checkScalar(n1, "integer", c(1, Inf))
checkScalar(n2, "integer", c(1, Inf))
checkScalar(nsim, "integer", c(1, Inf))
checkLengths(p1, p2)
x1 <- rbinom(prob=p1, size=n1, n=nsim)
x2 <- rbinom(prob=p2, size=n2, n=nsim)
scale <- match.arg(tolower(scale), c("difference", "rr", "or", "lnor"))
testBinomial(x1=x1, x2=x2, n1=n1, n2=n2, delta0=delta0, adj=adj,
chisq=chisq, scale=scale)
}
"testBinomial" <- function(x1, x2, n1, n2, delta0=0, chisq=0, adj=0,
scale="Difference", tol=.1e-10)
{
# check input arguments
checkVector(n1, "integer", c(1, Inf))
checkVector(n2, "integer", c(1, Inf))
checkVector(x1, "integer", c(0, Inf))
checkVector(x2, "integer", c(0, Inf))
checkVector(chisq, "integer", c(0, 1))
checkVector(adj, "integer", c(0, 1))
checkScalar(scale, "character")
scale <- match.arg(tolower(scale), c("difference", "rr", "or", "lnor"))
checkScalar(tol, "numeric", c(0, Inf), c(FALSE, TRUE))
checkLengths(n1, n2, x1, x2, delta0, chisq, adj, allowSingle=TRUE)
checkVector(n1 - x1, "integer", c(0, Inf))
checkVector(n2 - x2, "integer", c(0, Inf))
# make all vector arguments the same length (don't extend n1, n2)
len<-max(sapply(list(n1, n2, x1, x2, chisq, adj, delta0),length))
if (len > 1)
{ if (length(x1) == 1) p1<-array(x1,len)
if (length(x2) == 1) p2<-array(x2,len)
if (length(chisq) == 1) alpha<-array(chisq,len)
if (length(adj) == 1) beta<-array(adj,len)
if (length(delta0) == 1) delta0<-array(delta0,len)
}
ntot <- n1 + n2
xtot <- x1 + x2
# risk difference test - from Miettinen and Nurminen eqn (9)
if (scale == "difference")
{
checkVector(delta0, "numeric", c(-1, 1), c(FALSE, FALSE))
L2 <- (n1 + 2 * n2) * delta0 - ntot - xtot
L1 <- (n2 * delta0 - ntot - 2 * x2) * delta0 + xtot
L0 <- x2 * delta0 * (1 - delta0)
q <- (L2 / (3 * ntot)) ^ 3 - L1 * L2 / 6 / ntot ^ 2 + L0 / 2 / ntot
p <- (sign(q) + (q==0)) * sqrt((L2 / (3 * ntot)) ^ 2 - L1 / (3 * ntot))
a <- q / p ^ 3
a[a > 1] <- 1
a <- (pi + acos(a)) / 3
R0 <- 2 * p * cos(a) - L2 / 3 / ntot
R1 <- R0 + delta0
V <- (R1 * (1 - R1) / n1 + R0 * (1 - R0) / n2)
# V=0 only if no or all events and delta0=0
# in which case test statistic will be 0
V[V <= tol] <- 1
# compute z-statistic
z <- x1 / n1 - x2 / n2 - delta0
}
# relative risk test - from Miettinen and Nurminen eqn (10)
else if (scale == "rr")
{
checkVector(delta0, "numeric", c(-Inf, Inf), c(FALSE, FALSE))
delta0 <- exp(delta0) # value of 0 input represents equal rates
A <- delta0 * ntot
B <- -(n1 * delta0 + x1 + n2 + x2 * delta0)
C <- xtot
R0 <- ( - B - sqrt(B ^ 2 - 4 * A * C)) / 2 / A
R1 <- R0 * delta0
V <- R1 * (1 - R1) / n1 + delta0 ^ 2 * R0 * (1 - R0) / n2
# V=0 only if no events or (all events and delta0=1)
# in which case test statistic will be 0
V[V <= 0 | is.na(V)] <- 1
z <- x1 / n1 - x2 / n2 * delta0
}
# odds-ratio and log-odds-ratio
else
{
checkVector(delta0, "numeric", c(-Inf, Inf), c(FALSE, FALSE))
delta0 <- exp(delta0) # change from log odds-ratio to odds-ratio
A <- n2 * (delta0 - 1)
B <- n1 * delta0 + n2 - xtot * (delta0 - 1)
C <- -xtot
R0 <- ( - B + sqrt(B ^ 2 - 4 * A * C)) / 2 / A
R1 <- R0 * delta0 / (1 + R0 * (delta0 - 1))
# next 2 lines deal with length(delta0)>length(xtot) or length(ntot)
xtem <- (delta0 == delta0) * xtot
ntem<-(delta0==delta0)*ntot
R0[delta0 == 1] <- xtem[delta0 == 1] / ntem[delta0 == 1]
R1[delta0 == 1] <- R0[delta0 == 1]
# odds-ratio test - from Miettinen and Nurminen eqn (13)
if (scale == "or")
{
V <- 1 / (1 / n1 / R1 / (1-R1) + 1 / n2 / R0 / (1 - R0))
V[xtot == 0 | xtot == ntot] <- 1
z <- n1 * (x1 / n1 - R1)
}
# log-odds ratio - based on asymptotic distribution of log-odds
# see vignette
else if (scale == "lnor")
{
V <- 1 / n1 / R1 / (1-R1) + 1 / n2 / R0 / (1-R0)
V[xtot == 0 | xtot == ntot] <- 1
z <- log(x1 / (n1 - x1) / x2 * (n2 - x2) / delta0)
z[xtot == 0 | xtot == ntot] <- 0
z[xtot > 0 & xtot < ntot] <- z[xtot > 0 & xtot < ntot]
}
}
# do continuity correction, where required
one <- array(TRUE, max(length(V), length(adj)))
adj <- (adj==1) & one
V[adj] <- V[adj] * ntot / (ntot - 1)
# square z where required to get chi-square
one <- array(TRUE, max(length(z), length(chisq)))
chisq <- chisq & one
if (length(z) == 1) z <- z * one
z[chisq] <- z[chisq] ^ 2 / V[chisq]
z[ ! chisq] <- z[ ! chisq] / sqrt(V[ ! chisq])
z
}
###
# Hidden Functions
###
"bpdiff" <- function(delta, x1, x2, n1, n2, alpha=.05, adj=0, scale="Difference",
lower.tail=FALSE)
{
pnorm(testBinomial(x1, x2, n1, n2, delta0=delta, adj=adj, scale=scale),
lower.tail=lower.tail) - alpha / 2
}
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.