Designing two-stage trials with ordered categorical outcomes

Overview

OptOTrials constructs optimal one- and two-stage randomised trial designs for ordered categorical outcomes. The whole workflow runs through two functions:

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.

Step 1: check the proportional odds assumption

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.11841

When 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] NA

For that scenario use the Mann-Whitney-Wilcoxon test ("M") or the win odds test ("W") instead.

Step 2: build a design

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.488

The 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:

d$n1
#> [1] 28
d$t1f
#> [1] 0.083
d$effect
#>   boundary            stage threshold    scale effect
#> 1      t1f interim futility     0.083 win odds  1.181
#> 2       t2   final analysis     0.196 win odds  1.488

Step 3: operating characteristics

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.28

Both 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.

Step 4: comparing criteria

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      FALSE

Degenerate optima

Criterion 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] TRUE

Use on_degenerate = "none" to inspect the raw optimum, or "error" to make it fail loudly in an automated pipeline.

Choosing between asymptotic and exact approaches

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.

Moving from the earlier interface

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.