# ciBinomial roxy [sinew] ----
#' @seealso
#' \code{\link[stats]{Normal}},\code{\link[stats]{uniroot}}
#' @rdname varBinomial
#' @export
#' @importFrom stats qnorm uniroot
# ciBinomial function [sinew] ----
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)
< stats::qnorm(alpha / 2)) {
lower <- -1
}
else {
lower <- stats::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)
> -stats::qnorm(alpha / 2)) {
upper <- 1
}
else {
upper <- stats::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 <- stats::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 <- stats::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 <- stats::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 <- stats::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 roxy [sinew] ----
#' @rdname varBinomial
#' @export
#' @importFrom stats qnorm pnorm
# nBinomial function [sinew] ----
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 <- rep(p1, len)
}
if (length(p2) == 1) {
p2 <- rep(p2, len)
}
if (length(alpha) == 1) {
alpha <- rep(alpha, len)
}
if (length(beta) == 1) {
beta <- rep(beta, len)
}
if (length(delta0) == 1) {
delta0 <- rep(delta0, len)
}
if (length(ratio) == 1) {
ratio <- rep(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 <- stats::qnorm(1 - beta)
sided[sided != 2] <- 1
z.alpha <- stats::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 <- stats::pnorm(-(stats::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 <- stats::pnorm(-(stats::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 <- stats::pnorm(-(stats::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 roxy [sinew] ----
#' @rdname varBinomial
#' @export
#' @importFrom stats rbinom
# simBinomial function [sinew] ----
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 <- stats::rbinom(prob = p1, size = n1, n = nsim)
x2 <- stats::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 roxy [sinew] ----
#' @rdname varBinomial
#' @export
#' @importFrom stats rbinom
# testBinomial function [sinew] ----
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 <- rep(x1, len)
if (length(x2) == 1) p2 <- rep(x2, len)
if (length(chisq) == 1) alpha <- rep(chisq, len)
if (length(adj) == 1) beta <- rep(adj, len)
if (length(delta0) == 1) delta0 <- rep(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 <- rep(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 <- rep(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
}
# bpdiff function [sinew] ----
bpdiff <- function(delta, x1, x2, n1, n2, alpha = .05, adj = 0, scale = "Difference", lower.tail = FALSE) {
stats::pnorm(testBinomial(x1, x2, n1, n2, delta0 = delta, adj = adj, scale = scale),
lower.tail = lower.tail
) - alpha / 2
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.