OptOTrials constructs optimal one- and two-stage
randomised trial designs for ordered categorical outcomes. The whole
workflow runs through two functions:
rule() builds a decision rule;op() estimates its operating characteristics by
simulation.The test statistic and the monitoring scheme are chosen by argument, so you do not need to remember a different function name for each combination.
The score test (test = "S") is only appropriate when the
proportional odds assumption holds.
Proportional_odds_assumption() reports whether it does and
returns the implied common log odds ratio.
p1 <- c(0.075, 0.182, 0.319, 0.243, 0.015, 0.166)
p2 <- p2_fun(p1, log(3.06))
Proportional_odds_assumption(p1, p2)
#> The proportional odds assumption holds. Common log odds ratio = 1.11841 (odds ratio = 3.0600); the score test is appropriate.
#> [1] 1.11841When the assumption fails the function says so and warns, so the condition cannot pass unnoticed in a script:
q1 <- c(1/3, 1/3, 1/3)
q2 <- c(1/2, 1/3, 1/6)
res <- tryCatch(Proportional_odds_assumption(q1, q2), warning = function(w) NA)
#> The proportional odds assumption does not hold for the supplied probabilities. It is not advisable to use the score test; consider the Mann-Whitney-Wilcoxon test ("M") or the win odds test ("W").
res
#> [1] NAFor that scenario use the Mann-Whitney-Wilcoxon test
("M") or the win odds test ("W") instead.
p2_fun() is a helper in this package that constructs the
experimental-group category probabilities from the control probabilities
and a specified log odds ratio under the proportional odds model, as
used above.
d <- rule(alpha = 0.05, beta = 0.2, p1 = p1, p2 = p2,
test = "M", stopping = "F", criterion = 1)
d
#> Optimal design for an ordered categorical outcome
#> -------------------------------------------------
#> Test : Mann-Whitney-Wilcoxon test
#> Monitoring : two-stage, futility only
#> Criterion : 1 - minimise E(N | H0)
#> Error rates : alpha = 0.050, beta = 0.200 (power 80.0%)
#> Allocation : 1 : 1 (control : experimental)
#>
#> Decision rule
#> Stage 1: enrol 28 patients in total
#> stop for futility if T1 <= 0.0830
#> otherwise continue to stage 2
#> Stage 2: enrol up to 89 patients in total
#> declare superiority if T2 > 0.1960
#>
#> Boundaries on the effect-size scale
#> interim futility threshold 0.0830 -> win odds 1.181
#> final analysis threshold 0.1960 -> win odds 1.488The printed output names every quantity and additionally translates each boundary onto an effect-size scale, so the design can be communicated to non-statistician collaborators: for the score test the implied odds ratio, and for the rank-based tests the implied win odds.
Individual components are available by name, so no positional indexing is needed:
o <- op(d, nsim = 2000, seed = 1234)
o
#> Operating characteristics
#> -------------------------
#> Test : Mann-Whitney-Wilcoxon test
#> Monitoring : futility only
#> Criterion : 1
#> Replicates : 2000 (seed 1234)
#>
#> Scenario Reject H0 (MCSE) Stop fut. Stop super. Continue E(N)
#> H0 (type I) 0.038 (0.004) 0.653 0.000 0.347 49.14
#> Ha (power) 0.848 (0.008) 0.091 0.000 0.908 83.42
#>
#> EN0 = 49.14 ENa = 83.42 EN = 66.28Both hypotheses are evaluated in one call. Futility and superiority stopping probabilities are reported separately, and Monte Carlo standard errors are attached so that a deviation of the estimated type I error rate from the nominal level can be judged against simulation noise.
Supplying seed re-seeds immediately before each
scenario, so the result does not depend on how much randomness has
already been consumed.
design_table() sweeps a set of tests and criteria and
returns a data frame, re-seeding for every row so each row reproduces on
its own.
design_table(0.05, 0.2, p1, p2, tests = c("S", "M"), criteria = c(1, 3),
stopping = "F", nsim = 1000)
#> test criterion n1 t1f n2 t2 alpha se_alpha power se_power p_fut_H0
#> 1 S 1 29 0.654 91 4.378 0.048 0.0068 0.828 0.0119 0.652
#> 2 S 3 31 -0.510 67 3.753 0.049 0.0068 0.808 0.0125 0.370
#> 3 M 1 28 0.083 89 0.196 0.039 0.0061 0.848 0.0114 0.646
#> 4 M 3 30 -0.058 69 0.223 0.041 0.0063 0.826 0.0120 0.360
#> p_sup_H0 p_fut_Ha p_sup_Ha EN0 ENa EN degenerate
#> 1 0 0.106 0 50.576 84.428 67.5020 FALSE
#> 2 0 0.023 0 53.680 66.172 59.9260 FALSE
#> 3 0 0.092 0 49.594 83.388 66.4910 FALSE
#> 4 0 0.021 0 54.960 68.181 61.5705 FALSECriterion 2 minimises the expected sample size under the alternative. This can drive the optimum to a stage-1 size of 1 with an extreme interim boundary: the criterion rewards designs under which continuation to the final analysis is almost certain when the treatment works, so the optimiser commits essentially all information to the final analysis and the design becomes two-stage in name only.
rule() detects this. By default it warns and substitutes
the single-stage design:
d2 <- rule(0.05, 0.2, p1, p2, test = "M", stopping = "F", criterion = 2)
#> Warning: The optimum for test = "M", criterion = 2 has collapsed to a
#> degenerate two-stage design (stage-1 sample size n1 = 1 is below min_n1 = 10;
#> the interim futility boundary t1f = -4.715 is so extreme that the trial stops
#> early with probability 0.0000 under H0). Returning the single-stage design
#> instead; set on_degenerate = "none" to obtain the raw optimum.
d2$substituted
#> [1] TRUEUse on_degenerate = "none" to inspect the raw optimum,
or "error" to make it fail loudly in an automated
pipeline.
The designs rely on the asymptotic normality of the test statistics. This is adequate at the stage-1 sizes arising in typical applications, but the approximation degrades when the stage-1 size is small or when some outcome categories are rare. In those situations an exact approach is preferable; see the package reference manual for the relevant references.
Versions up to 1.0.2 provided a separate function for each
combination of test and stopping rule (ruleF(),
ruleFS(), op.F(), op.FS(),
op.1stage() and the Decision_rule_* family).
These are deprecated in this version: each still works and still returns
exactly what it always returned, but warns and names its replacement.
They will be removed in version 1.1.0. The correspondence is tabulated
in help("OptOTrials-deprecated"):
| Deprecated | Replacement |
|---|---|
ruleF(), ruleFS() |
rule(..., stopping = "F"/"FS") |
Decision_rule_{S,M,W}.F / .FS /
_1stage |
rule(..., test = ..., stopping = ...) |
op.F(), op.FS(),
op.1stage() |
op(design, nsim, seed) |
Two habits change. rule() returns a named object, so
res[3] becomes res$n1, and the object is
passed straight to op() instead of its elements being
transcribed. op() evaluates both hypotheses in one call, so
the old pattern of calling the simulator twice, once with
(p1, p2) and once with (p1, p1), is no longer
needed.