Skip to content

Commit d0c31a5

Browse files
nanxstatskeaven
nanxstats
and
keaven
authored
Add vignette for binomial two-arm trial design and analysis (#209)
* Adding binomial vignette and summary functions * ongoing updates * Fix unit tests and remove object files * Near final? * Updated documentation * binomialPowerTable added * Minor updates * Minor updates * Increased accuracy in uniroot calls to tol = 1e-10 * range checks added and documented * reversion * minor code fix * revert to master * plot scale fix * range documentation * updates * minor edits to documentation * Remove spurious new rules * Sort and align BibTeX items * Reformat nbinomial input tests * Run styler * Improve document style for `binomialPowerTable()` * Improve binomial two sample vignette formatting * Add dplyr variable from `binomialPowerTable()` --------- Co-authored-by: keaven <[email protected]> Co-authored-by: keaven <[email protected]>
1 parent 4f717c8 commit d0c31a5

File tree

9 files changed

+782
-78
lines changed

9 files changed

+782
-78
lines changed

‎NAMESPACE‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,7 @@ export(Power.ssrCP)
2424
export(as_gt)
2525
export(as_rtf)
2626
export(as_table)
27+
export(binomialPowerTable)
2728
export(binomialSPRT)
2829
export(checkLengths)
2930
export(checkRange)

‎R/binomialPowerTable.R‎

Lines changed: 156 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,156 @@
1+
#' Power Table for Binomial Tests
2+
#'
3+
#' Creates a power table for binomial tests with various control group response rates and treatment effects.
4+
#' The function can compute power and Type I error either analytically or through simulation.
5+
#' With large simulations, the function is still fast and can produce exact power values to within
6+
#' simulation error.
7+
#'
8+
#' @param pC Vector of control group response rates.
9+
#' @param delta Vector of treatment effects (differences in response rates).
10+
#' @param n Total sample size.
11+
#' @param ratio Ratio of experimental to control sample size.
12+
#' @param alpha Type I error rate.
13+
#' @param delta0 Non-inferiority margin.
14+
#' @param scale Scale for the test
15+
#' (\code{"Difference"}, \code{"RR"}, or \code{"OR"}).
16+
#' @param failureEndpoint Logical indicating if the endpoint is a
17+
#' failure (\code{TRUE}) or success (\code{FALSE}).
18+
#' @param simulation Logical indicating whether to use simulation (\code{TRUE})
19+
#' or analytical (\code{FALSE}) power calculation.
20+
#' @param nsim Number of simulations to run when \code{simulation = TRUE}.
21+
#' @param adj Use continuity correction for the testing (default is 0;
22+
#' only used if \code{simulation = TRUE}).
23+
#' @param chisq Chi-squared value for the test (default is 0;
24+
#' only used if \code{simulation = TRUE}).
25+
#'
26+
#' @return A data frame containing:
27+
#' \describe{
28+
#' \item{\code{pC}}{Control group response or failure rate.}
29+
#' \item{\code{delta}}{Treatment effect.}
30+
#' \item{\code{pE}}{Experimental group response or failure rate.}
31+
#' \item{\code{Power}}{Power for the test (asymptotic or simulated).}
32+
#' }
33+
#'
34+
#' @details
35+
#' The function \code{binomialPowerTable()} creates a grid of all combinations of control group response rates and treatment effects.
36+
#' All out of range values (i.e., where the experimental group response rate is not between 0 and 1) are filtered out.
37+
#' For each combination, it computes the power either analytically using \code{nBinomial()} or through
38+
#' simulation using \code{simBinomial()}.
39+
#' When using simulation, the \code{simPowerBinomial()} function (not exported) is called
40+
#' internally to perform the simulations.
41+
#' Assuming \eqn{p} is the true probability of a positive test, the simulation standard error is
42+
#' \deqn{\text{SE} = \sqrt{p(1 - p) / \text{nsim}}.}
43+
#' For example, when approximating an underlying Type I error rate of 0.025, the simulation standard error is
44+
#' 0.000156 with 1000000 simulations and the approximated power 95% confidence interval
45+
#' is 0.025 +/- 1.96 * SE = 0.025 +/- 0.000306.
46+
#'
47+
#' @seealso \code{\link{nBinomial}}, \code{\link{simBinomial}}
48+
#'
49+
#' @rdname binomialPowerTable
50+
#'
51+
#' @export
52+
#'
53+
#' @examples
54+
#' # Create a power table with analytical power calculation
55+
#' power_table <- binomialPowerTable(
56+
#' pC = c(0.8, 0.9),
57+
#' delta = seq(-0.05, 0.05, 0.025),
58+
#' n = 70
59+
#' )
60+
#'
61+
#' # Create a power table with simulation
62+
#' power_table_sim <- binomialPowerTable(
63+
#' pC = c(0.8, 0.9),
64+
#' delta = seq(-0.05, 0.05, 0.025),
65+
#' n = 70,
66+
#' simulation = TRUE,
67+
#' nsim = 10000
68+
#' )
69+
binomialPowerTable <- function(
70+
pC = c(.8, .9, .95),
71+
delta = seq(-.05, .05, .025),
72+
n = 70,
73+
ratio = 1,
74+
alpha = .025,
75+
delta0 = 0,
76+
scale = "Difference",
77+
failureEndpoint = TRUE,
78+
simulation = FALSE,
79+
nsim = 1000000,
80+
adj = 0,
81+
chisq = 0) {
82+
# Create a grid of all combinations of pC and delta
83+
pC_grid <- expand.grid(pC = pC, delta = delta)
84+
# Compute the experimental group response rate
85+
pC_grid <- pC_grid %>%
86+
mutate(pE = pC + delta) %>%
87+
filter(pE < 1 & pE > 0 & pC < 1 & pC > 0) %>% # Filter out of range values
88+
select(pC, delta, pE)
89+
90+
# Compute the probability of non-inferiority using the nBinomial function
91+
if (!simulation) {
92+
pC_grid$Power <- mapply(function(pC, pE, delta) {
93+
if (n > 0 && pE == pC && delta0 == 0) {
94+
return(alpha)
95+
}
96+
if (failureEndpoint) {
97+
gsDesign::nBinomial(
98+
p1 = pE, p2 = pC, delta0 = delta0,
99+
n = n, alpha = alpha, scale = scale,
100+
ratio = ratio, sided = 1
101+
)
102+
} else {
103+
gsDesign::nBinomial(
104+
p1 = pC, p2 = pE, delta0 = delta0,
105+
n = n, alpha = alpha, scale = scale,
106+
ratio = 1 / ratio, sided = 1
107+
)
108+
}
109+
}, pC_grid$pC, pC_grid$pE, pC_grid$delta)
110+
return(pC_grid)
111+
} else {
112+
simPowerBinomial(pC_grid, n, ratio, alpha, delta0, scale, failureEndpoint, nsim)
113+
}
114+
}
115+
116+
simPowerBinomial <- function(
117+
longTable,
118+
n = 70,
119+
ratio = 1,
120+
alpha = .025,
121+
delta0 = 0,
122+
scale = "Difference",
123+
failureEndpoint = TRUE,
124+
nsim = 10000) {
125+
# Get sample size per arm
126+
nC <- round(n / (1 + ratio))
127+
nE <- n - nC
128+
129+
# Initialize vector for simulation results
130+
Power <- numeric(nrow(longTable))
131+
132+
# Loop through each row in the table
133+
for (i in 1:nrow(longTable)) {
134+
pC <- longTable$pC[i]
135+
pE <- longTable$pE[i]
136+
if (failureEndpoint) {
137+
sim <- gsDesign::simBinomial(
138+
n1 = nE, n2 = nC, p1 = pE, p2 = pC,
139+
delta0 = delta0, nsim = nsim,
140+
scale = scale
141+
)
142+
Power[i] <- mean(sim >= qnorm(1 - alpha))
143+
} else {
144+
sim <- gsDesign::simBinomial(
145+
n1 = nC, n2 = nE, p1 = pC, p2 = pE,
146+
delta0 = delta0, nsim = nsim,
147+
scale = scale
148+
)
149+
Power[i] <- mean(sim >= qnorm(1 - alpha))
150+
}
151+
}
152+
153+
# Add simulation results to the table
154+
longTable$Power <- Power
155+
return(longTable)
156+
}

‎R/globals.R‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,9 @@ utils::globalVariables(
22
unique(
33
c(
44
# From `as_gt.gsBinomialExactTable()`
5-
c("Lower", "Upper", "en")
5+
c("Lower", "Upper", "en"),
6+
# From `binomialPowerTable()`
7+
c("pE")
68
)
79
)
810
)

‎R/gsBinomial.R‎

Lines changed: 24 additions & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -22,12 +22,10 @@ ciBinomial <- function(x1, x2, n1, n2, alpha = .05, adj = 0, scale = "Difference
2222

2323
if (delta == -1) {
2424
lower <- -1
25-
}
26-
else if (testBinomial(delta0 = -.9999, x1 = x1, x2 = x2, n1 = n1, n2 = n2, adj = adj)
25+
} else if (testBinomial(delta0 = -.9999, x1 = x1, x2 = x2, n1 = n1, n2 = n2, adj = adj)
2726
< stats::qnorm(alpha / 2)) {
2827
lower <- -1
29-
}
30-
else {
28+
} else {
3129
lower <- stats::uniroot(bpdiff,
3230
interval = c(-.9999, delta), x1 = x1, x2 = x2,
3331
n1 = n1, n2 = n2, adj = adj, alpha = alpha
@@ -36,20 +34,17 @@ ciBinomial <- function(x1, x2, n1, n2, alpha = .05, adj = 0, scale = "Difference
3634

3735
if (delta == 1) {
3836
upper <- 1
39-
}
40-
else if (testBinomial(delta0 = .9999, x1 = x1, x2 = x2, n1 = n1, n2 = n2, adj = adj)
37+
} else if (testBinomial(delta0 = .9999, x1 = x1, x2 = x2, n1 = n1, n2 = n2, adj = adj)
4138
> -stats::qnorm(alpha / 2)) {
4239
upper <- 1
43-
}
44-
else {
40+
} else {
4541
upper <- stats::uniroot(bpdiff,
4642
interval = c(delta, .9999), lower.tail = TRUE,
4743
x1 = x1, x2 = x2, n1 = n1, n2 = n2, adj = adj,
4844
alpha = alpha
4945
)$root
5046
}
51-
}
52-
else if (scale == "rr") {
47+
} else if (scale == "rr") {
5348
if (x1 == 0) {
5449
lower <- 0
5550
} else {
@@ -69,8 +64,7 @@ ciBinomial <- function(x1, x2, n1, n2, alpha = .05, adj = 0, scale = "Difference
6964
)$root
7065
upper <- exp(upper)
7166
}
72-
}
73-
else {
67+
} else {
7468
if (x1 == 0 || x2 == n2) {
7569
lower <- -Inf
7670
} else {
@@ -178,42 +172,36 @@ nBinomial <- function(p1, p2, alpha = 0.025, beta = 0.1, delta0 = 0, ratio = 1,
178172
if (outtype == 2) {
179173
return(data.frame(n1 = n / (ratio + 1), n2 = ratio *
180174
n / (ratio + 1)))
181-
}
182-
else if (outtype == 3) {
175+
} else if (outtype == 3) {
183176
return(data.frame(
184177
n = n, n1 = n / (ratio + 1),
185178
n2 = ratio * n / (ratio + 1), alpha = alpha,
186179
sided = sided, beta = beta, Power = 1 - beta,
187180
sigma0 = sigma0, sigma1 = sigma1, p1 = p1,
188181
p2 = p2, delta0 = delta0, p10 = p10, p20 = p20
189182
))
190-
}
191-
else {
183+
} else {
192184
return(n = n)
193185
}
194-
}
195-
else {
186+
} else {
196187
pwr <- stats::pnorm(-(stats::qnorm(1 - alpha / sided) - sqrt(n) *
197188
((p1 - p2 - delta0) / sigma0)) * sigma0 / sigma1)
198189
if (outtype == 2) {
199190
return(data.frame(n1 = n / (ratio + 1), n2 = ratio *
200191
n / (ratio + 1), Power = pwr))
201-
}
202-
else if (outtype == 3) {
192+
} else if (outtype == 3) {
203193
return(data.frame(
204194
n = n, n1 = n / (ratio + 1),
205195
n2 = ratio * n / (ratio + 1), alpha = alpha,
206196
sided = sided, beta = 1 - pwr, Power = pwr,
207197
sigma0 = sigma0, sigma1 = sigma1, p1 = p1,
208198
p2 = p2, delta0 = delta0, p10 = p10, p20 = p20
209199
))
210-
}
211-
else {
200+
} else {
212201
return(Power = pwr)
213202
}
214203
}
215-
}
216-
else if (scale == "rr") {
204+
} else if (scale == "rr") {
217205
RR <- exp(delta0)
218206
if (min(abs(p1 / p2 - RR)) < 1e-07) {
219207
stop("p1/p2 may not equal exp(delta0) when scale=\"RR\"")
@@ -236,42 +224,36 @@ nBinomial <- function(p1, p2, alpha = 0.025, beta = 0.1, delta0 = 0, ratio = 1,
236224
n1 = n / (ratio + 1),
237225
n2 = ratio * n / (ratio + 1)
238226
))
239-
}
240-
else if (outtype == 3) {
227+
} else if (outtype == 3) {
241228
return(data.frame(
242229
n = n, n1 = n / (ratio + 1),
243230
n2 = ratio * n / (ratio + 1), alpha = alpha,
244231
sided = sided, beta = beta, Power = 1 - beta,
245232
sigma0 = sigma0, sigma1 = sigma1, p1 = p1,
246233
p2 = p2, delta0 = delta0, p10 = p10, p20 = p20
247234
))
248-
}
249-
else {
235+
} else {
250236
return(n = n)
251237
}
252-
}
253-
else {
238+
} else {
254239
pwr <- stats::pnorm(-(stats::qnorm(1 - alpha / sided) - sqrt(n) *
255240
((p1 - p2 * RR) / sigma0)) * sigma0 / sigma1)
256241
if (outtype == 2) {
257242
return(data.frame(n1 = n / (ratio + 1), n2 = ratio *
258243
n / (ratio + 1), Power = pwr))
259-
}
260-
else if (outtype == 3) {
244+
} else if (outtype == 3) {
261245
return(data.frame(
262246
n = n, n1 = n / (ratio + 1),
263247
n2 = ratio * n / (ratio + 1), alpha = alpha,
264248
sided = sided, beta = 1 - pwr, Power = pwr,
265249
sigma0 = sigma0, sigma1 = sigma1, p1 = p1,
266250
p2 = p2, delta0 = delta0, p10 = p10, p20 = p20
267251
))
268-
}
269-
else {
252+
} else {
270253
return(Power = pwr)
271254
}
272255
}
273-
}
274-
else {
256+
} else {
275257
OR <- exp(-delta0)
276258
if (min(abs(p1 / (1 - p1) / p2 * (1 - p2) * OR - 1)) < 1e-07) {
277259
stop("p1/(1-p1)/p2*(1-p2) may not equal exp(delta0) when scale=\"OR\"")
@@ -290,38 +272,33 @@ nBinomial <- function(p1, p2, alpha = 0.025, beta = 0.1, delta0 = 0, ratio = 1,
290272
if (outtype == 2) {
291273
return(data.frame(n1 = n / (ratio + 1), n2 = ratio *
292274
n / (ratio + 1)))
293-
}
294-
else if (outtype == 3) {
275+
} else if (outtype == 3) {
295276
return(data.frame(
296277
n = n, n1 = n / (ratio + 1),
297278
n2 = ratio * n / (ratio + 1), alpha = alpha,
298279
sided = sided, beta = beta, Power = 1 - beta,
299280
sigma0 = sigma0, sigma1 = sigma1, p1 = p1,
300281
p2 = p2, delta0 = delta0, p10 = p10, p20 = p20
301282
))
302-
}
303-
else {
283+
} else {
304284
return(n = n)
305285
}
306-
}
307-
else {
286+
} else {
308287
pwr <- stats::pnorm(-(stats::qnorm(1 - alpha / sided) - sqrt(n) *
309288
(log(OR / p2 * (1 - p2) * p1 / (1 - p1)) / sigma0)) *
310289
sigma0 / sigma1)
311290
if (outtype == 2) {
312291
return(data.frame(n1 = n / (ratio + 1), n2 = ratio *
313292
n / (ratio + 1), Power = pwr))
314-
}
315-
else if (outtype == 3) {
293+
} else if (outtype == 3) {
316294
return(data.frame(
317295
n = n, n1 = n / (ratio + 1),
318296
n2 = ratio * n / (ratio + 1), alpha = alpha,
319297
sided = sided, beta = 1 - pwr, Power = pwr,
320298
sigma0 = sigma0, sigma1 = sigma1, p1 = p1,
321299
p2 = p2, delta0 = delta0, p10 = p10, p20 = p20
322300
))
323-
}
324-
else {
301+
} else {
325302
return(Power = pwr)
326303
}
327304
}
@@ -464,7 +441,7 @@ testBinomial <- function(x1, x2, n1, n2, delta0 = 0, chisq = 0, adj = 0, scale =
464441
chisq <- chisq & one
465442
if (length(z) == 1) z <- z * one
466443
z[chisq] <- z[chisq]^2 / V[chisq]
467-
z[ !chisq] <- z[ !chisq] / sqrt(V[ !chisq])
444+
z[!chisq] <- z[!chisq] / sqrt(V[!chisq])
468445

469446
z
470447
}

‎_pkgdown.yml‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -36,6 +36,7 @@ reference:
3636
- testBinomial
3737
- varBinomial
3838
- ciBinomial
39+
- binomialPowerTable
3940
- title: Time-to-Event Endpoint Design
4041
contents:
4142
- nEvents
@@ -152,4 +153,5 @@ articles:
152153
- hGraph
153154
- GraphicalMultiplicity
154155
- VaccineEfficacy
156+
- binomialTwoSample
155157
- binomialSPRTExample

0 commit comments

Comments
 (0)