Simulating 3+3 Operating Characteristics

Why Simulate a Rule-Based Design?

The 3+3 design is a deterministic decision rule applied to outcomes that are random. Its operating characteristics are therefore probabilities over all possible trial paths: for example, the probability of selecting each dose, selecting a dose whose true toxicity rate is above a chosen target, or stopping without a recommendation. Simulation estimates these probabilities by repeatedly generating hypothetical trials under a specified set of dose-specific dose-limiting toxicity (DLT) probabilities.

The simulation does not predict what will happen in a particular real trial. It answers a conditional question: if the assumed DLT probabilities and model assumptions were true, how often would this protocol produce each outcome? The results depend on the scenario supplied, the exact protocol rules coded, and the number of simulated trials.

Define the Algorithm and the Scenario

This tutorial uses the conventional rule with complete DLT ascertainment before each decision, no unevaluable participants, and no dose skipping. The underlying participant outcomes are independent Bernoulli variables at each dose, with DLT probability \(p_j\); outcomes are identically distributed for participants treated at the same dose. The scenario probabilities are treated as fixed and known for the purpose of evaluating the design. A monotone scenario, \(p_1 \leq p_2 \leq \cdots \leq p_K\), is natural for ordered doses, although the algorithm can be run for nonmonotone probabilities too. In a nonmonotone scenario, dose-based interpretations such as “highest acceptable dose” require particular care.

For interpretation, define the true target dose \(j^*\) as the highest dose with \(p_j \leq p_T\), when at least one dose meets that condition. A selected dose above \(j^*\) is an overdose selection; one below \(j^*\) is an underdose selection. Selecting no dose is reported separately, not counted as either kind of selection. If no dose meets the target criterion, or if every dose meets it, describe that scenario explicitly rather than forcing an ambiguous underdose or overdose label.

Calculate a Branch Probability Before Simulating

A short calculation checks the logic of the first decision. At dose 1, the chance of 0 DLTs among three participants is \((1-p_1)^3\). The chance of exactly 1 DLT is \(3p_1(1-p_1)^2\). After exactly 1/3, escalation occurs only if the next three have 0 DLTs, with probability \((1-p_1)^3\), under the independence assumption. Thus, the probability of escalating after evaluating dose 1 is

The first term is the 0/3 path; the second is the 1/3 path followed by 0/3 in the added cohort. All other paths stop at dose 1. This identity is also a useful check on a simulation: escalation and stopping probabilities at the first dose must sum to 1.

A Reproducible Simulation Procedure

One simulated trial is generated by applying the protocol to newly generated participant outcomes. At any dose with 1/3 DLTs, the next cohort is generated at the same dose, not at the next dose. At a dose with 0/3, no additional cohort is enrolled there. Each simulated trial ends with a selected dose (possibly none), a dose at which toxicity stopped escalation or the trial reached its maximum dose, and a path of enrolled cohorts.

The following base R code gives an executable implementation. The output value 0 denotes no selected dose; positive values denote the selected dose index. The recorded stopping dose is the dose where a toxicity boundary was met, or the highest dose if the trial reached the top under the rule.

set.seed(552); sim1 <- function(p) { K <- length(p); selected <- 0L; stopdose <- K; for (j in seq_len(K)) { x <- rbinom(1, 3, p[j]); if (x == 0) { if (j == K) { selected <- j; break } } else if (x == 1) { y <- rbinom(1, 3, p[j]); if (x + y <= 1) { if (j == K) { selected <- j; break } } else { selected <- j - 1L; stopdose <- j; break } } else { selected <- j - 1L; stopdose <- j; break } }; c(selected = selected, stopdose = stopdose) }; p <- c(0.10, 0.25, 0.40); N <- 10000L; z <- t(replicate(N, sim1(p))); table(factor(z[, "selected"], levels = 0:3))

The seed makes the random-number stream reproducible within the same R environment and generator settings; it does not make the simulation estimates exact probabilities. To estimate stopping at each dose, tabulate the recorded `stopdose` values as well. In a different programming language, the same algorithm and a documented seed should be used. Check that probabilities supplied to the code lie between 0 and 1, that the vector length matches the intended dose grid, and that each cohort follows the protocol.

Monte Carlo Error and Basic Diagnostics

For a prespecified event, such as selecting an overdose, each simulated trial either has the event or does not. Under independent simulation runs, the event indicator is Bernoulli with probability \(\theta\). If \(\hat{\theta}\) is the proportion of \(N\) runs with the event, its estimated Monte Carlo standard error is

This formula follows because the sample mean of \(N\) independent Bernoulli indicators has variance \(\theta(1-\theta)/N\); replacing the unknown \(\theta\) by \(\hat{\theta}\) estimates that variance. It measures simulation noise, not uncertainty in the true \(p_j\) values and not sampling uncertainty in a future clinical trial outcome.

Useful diagnostics include verifying that all selection frequencies, including “none,” sum to 1; checking the first-dose branch calculation against a direct simulation; and examining whether escalation, additional cohorts, and stopping occur exactly as the written rule specifies. For \(N\) independent trials, the standard error for a very rare event may still be substantial relative to its estimated probability. If an event is never observed, the estimate is zero, not proof that its probability is zero. Increase \(N\), or use exact path enumeration for a small dose grid when exact calculation is practical.

Interpretation, Alternatives, and Limitations

Simulation and exact enumeration are complementary ways to calculate operating characteristics. Exact enumeration lists all possible outcome paths and sums their probabilities; it can avoid Monte Carlo error and is valuable for checking a small design. The number of paths can grow quickly as doses and cohorts increase, and coding every branch can itself be error-prone. Simulation is often simpler to extend and makes it straightforward to record quantities such as enrollment and stopping dose, but its estimates have Monte Carlo error. Neither method repairs a misspecified scenario or an incorrectly implemented protocol.

The simulation conditions on the chosen \(p_j\) values. It does not establish that those values are clinically plausible, account for uncertainty about them, or estimate toxicity from actual trial data. The simple Bernoulli model also omits delayed DLT ascertainment, dependence among participants, changing eligibility or follow-up, and other features unless they are explicitly incorporated. Results are specific to the exact escalation and selection conventions coded. For example, changing how the highest dose is handled changes the operating characteristics.

A common mistake is to report the proportion of trials that stopped at the top dose as the probability of selecting the top dose. These are not necessarily the same event: a trial can reach the top dose and then stop for toxicity, recommending the dose below. Correction: record “selected dose” and “stopping dose” separately, and define every reported event before running the simulation. Similarly, do not call every trial ending without reaching the top a “no-dose” outcome; the design may have selected a lower dose after a toxicity boundary.