# Incertus.jl – The Julia Lego Blocks for Randomized Clinical Trial Designs

Yevgen Ryeznik
Department of Pharmacy
Uppsala University
Uppsala, Sweden
<EMAIL_ADDRESS>
& Oleksandr Sverdlov
Novartis Pharmaceuticals Corporation
East Hannover, NJ, USA
<EMAIL_ADDRESS>

###### Abstract

In this paper, we present Insertus.jl, the Julia package that can help the
user generate a randomization sequence of a given length for a multi-arm trial
with a pre-specified target allocation ratio and assess the operating
characteristics of the chosen randomization method through Monte Carlo
simulations. The developed package is computationally efficient, and it can be
invoked in R. Furthermore, the package is open-ended – it can flexibly
accommodate new randomization procedures and evaluate their statistical
properties via simulation. It may be also helpful for validating other
randomization methods for which software is not readily available. In summary,
Insertus.jl can be used as “Lego Blocks” to construct a fit-for-purpose
randomization procedure for a given clinical trial design.

_Keywords_ multi-arm randomized clinical trial $\cdot$ unequal allocation
$\cdot$ randomization software

##  1 Introduction

Randomization in the context of a randomized controlled trial (RCT) refers to
a process of sequentially assigning eligible participants to different
treatment groups through a process governed by chance (Wittes, 2005).
Randomization helps to minimize experimental bias, ensure the comparability of
treatment groups at baseline, and allow for a precise and valid treatment
comparison (Rosenberger and Lachin, 2015). While randomization is an important
and desirable element in studies across the entire research and development
continuum, it is essential in confirmatory RCTs where the goal is to test some
pre-specified clinical research hypothesis (Green, 2002).

A classical RCT design is a 1:1, parallel-group, randomized, double-blind,
placebo-controlled study (Friedman et al., 2015). However, modern RCTs
frequently utilize multi-arm parallel group designs with equal or unequal
randomization ratios. Some examples of such study designs are randomized phase
2 dose range finding studies (Bretz et al., 2008), phase 2/3 seamless designs
(Schmidli et al., 2006), platform trials to evaluate multiple experimental
treatments in a perpetual manner (The Adaptive Platform Trials Coalition,
2019), amongst others. The choice of a randomization method to implement the
target allocation ratio for a given sample size and other trial parameters
requires careful consideration. There is no shortage of randomization designs
with established statistical properties (Rosenberger and Lachin, 2015);
however, choosing one that is fit for purpose may be a challenging task. Some
authors developed systematic approaches for the choice of an adequate
randomization procedure in the context of a 1:1 RCT (Hilgers et al., 2017;
Berger et al., 2021; Sverdlov et al., 2024).

The evaluation of statistical properties of randomization designs can be done
both theoretically and using simulations. For 1:1 RCTs, one useful statistical
software package for randomization is randomizeR (Uschner et al., 2018). For
more complex settings, such as multi-arm RCTs, there is a dearth of
statistical software that would enable the implementation of a chosen
randomization method. Several recent references provide relevant simulation
evidence on randomization in three-arm or four-arm trials (Ryeznik and
Sverdlov, 2018; Ryeznik et al., 2018a, b). Given the increased popularity and
use of multi-arm clinical trial designs, there is a strong need for a flexible
software package that enables fast and efficient generation of treatment
randomization sequences for a multi-arm trial with a pre-specified target
allocation ratio and enables the investigation of statistical properties of
the chosen randomization method(s) through Monte Carlo simulations.

In this paper, we present Insertus.jl, the Julia package that can help
clinical trial researchers to address the aforementioned goals. The highlights
of Insertus.jl are:

  * •

It can be used to generate randomization sequences using a chosen method for a
two-arm or a multi-arm trial with equal or unequal allocation ratios.

  * •

It enables fast and efficient Monte Carlo simulation studies to compare
operating characteristics of various randomization designs for the chosen
trial parameters.

  * •

The package is open-ended – it can flexibly integrate new procedures and
evaluate their statistical properties via simulation. It can be also used as a
tool to validate other randomization methods for which software is not readily
available.

Section 2 provides the relevant statistical background and describes the
randomization methods implemented in Incertus.jl. Section 3 describes the
Julia functions and presents some illustrative examples. Section 4 provides a
summary and describes some potential extensions of the package.

##  2 Methods

Consider a parallel-group, randomized clinical trial to compare the effects of
$K\geq 2$ treatments. Let $w_{1}:w_{2}:\ldots:w_{K}$ denote the target
allocation ratio for the trial, where $w_{i}$’s are some positive, not
necessarily equal numbers. Frequently, $w_{i}$’s are chosen to be integers
with the common divisor of 1. For example, in a two-arm trial, $1{:}1$
allocation means that each treatment group should have an equal number of
patients, whereas $2{:}1$ allocation means the first group should have twice
as many patients as the second group. In some trials, an investigator may
pursue a non-integer (irrational-valued) allocation ratio, e.g., the square
root allocation $\sqrt{2}{:}1{:}1$ which is optimal for inference on pairwise
comparisons involving treatment 1 as the reference in a three-arm trial
(Fleiss, 1986).

Let $\rho_{k}=w_{k}/\sum_{i=1}^{K}w_{i}$ denote the target allocation
proportion for the $k$th treatment group. In our examples,
$\rho_{1}=\rho_{2}=1/2$ for 1:1 allocation; $\rho_{1}=2/3$, $\rho_{1}=1/3$ for
2:1 allocation; and $\rho_{1}=\frac{\sqrt{2}}{\sqrt{2}+2}$,
$\rho_{2}=\rho_{3}=\frac{1}{\sqrt{2}+2}$ for $\sqrt{2}{:}1{:}1$ allocation. We
will use the notation $\rho_{1}{:}\rho_{2}{:}\ldots{:}\rho_{K}$, where
$\rho_{i}\in(0,1)$ and $\sum_{i=1}^{K}\rho_{i}=1$ to refer to any fixed, pre-
specified target allocation proportions for a $K\geq 2$-arm trial.

In a clinical trial, $n$ eligible patients are enrolled sequentially and must
be randomly assigned to $K$ treatments in the desired allocation ratio to
achieve the sample sizes per group as $n_{k}\approx n\rho_{k}$ with
$\sum_{k=1}^{K}n_{k}=n$. For the $j$th patient in the sequence
($j=1,\ldots,n$), let
$\boldsymbol{\delta}_{j}=(\delta_{j1},\ldots,\delta_{jK})^{T}$ denote a
$K\times 1$ random vector of treatment assignment such that only one entry of
$\boldsymbol{\delta}_{j}$ is equal to 1, and the rest are zeros. For example,
$\delta_{jk}=1$ and $\delta_{jk^{\prime}}=0$ for $k^{\prime}\neq k$ indicates
that the $j$th patient is assigned to the $k$th treatment group. This can be
also written as $\boldsymbol{\delta}_{j}=\boldsymbol{e}_{k}$, where
$\boldsymbol{e}_{k}$ ($k=1,\ldots K$) is a unit vector of zeros with a 1 in
the $k$th position. Let
$\boldsymbol{\Delta}_{n}=(\boldsymbol{\delta}_{1},\boldsymbol{\delta}_{2},\ldots,\boldsymbol{\delta}_{n})^{T}$
denote the $n\times K$ matrix of treatment assignments. Note that for any
$j\geq 1$, $N_{k}(j)=\sum_{i=1}^{j}\delta_{ik}$ is the sample size for the
$k$th treatment group after $j$ allocation steps. In particular, $N_{k}(n)$ is
the final sample size for the $k$th group, with $\sum_{k=1}^{K}N_{k}(n)=n$.

A _restricted_ randomization procedure for a $K$-arm trial with the target
allocation ratio $\rho_{1}{:}\rho_{2}{:}\ldots{:}\rho_{K}$ is defined by
specifying the $j$th patient’s ($j=1,\ldots,n$) conditional randomization
probabilities, $\phi_{jk}=\Pr(\textrm{Patient}\;j\;\textrm{is assigned to
treatment}\;k|\boldsymbol{\delta}_{1},\ldots,\boldsymbol{\delta}_{j-1})$, as
follows:

$\phi_{1k}=\rho_{k}\quad\textrm{and}\quad\phi_{jk}=\Pr(\boldsymbol{\delta}_{j}=\boldsymbol{e}_{k}|\boldsymbol{\delta}_{1},\ldots,\boldsymbol{\delta}_{j-1}),\;j\geq
2\;\textrm{and}\;k=1,\ldots,K.$ (1)

Equation (1) implies that the first patient in the sequence is randomized to
treatments with probabilities $(\rho_{1},\ldots,\rho_{K})$, and thereafter the
$j$th patient ($j\geq 2$) is randomized to treatments with probabilities
($\phi_{j1},\ldots,\phi_{jK}$) where $0\leq\phi_{jk}\leq 1$,
$\sum_{k=1}^{K}\phi_{jk}=1$ and $\phi_{jk}$ is the conditional probability of
assigning the $j$th patient to treatment $k$ given the past assignments
$\boldsymbol{\delta}_{1},\ldots,\boldsymbol{\delta}_{j-1}$. Denote
$\boldsymbol{\phi}_{j}=(\phi_{j1},\ldots,\phi_{jK})^{T}$ the $K\times 1$
vector and
$\boldsymbol{\Phi}_{n}=(\boldsymbol{\phi}_{1},\boldsymbol{\phi}_{2},\ldots,\boldsymbol{\phi}_{K})^{T}$
the $n\times K$ matrix of conditional randomization probabilities for $n$
patients in the trial.

The choice of $\boldsymbol{\phi}_{j}$’s determines the probabilistic structure
of a randomization procedure and its statistical properties. Some important
and desirable features of a randomization procedure are:

  * (a)

It should have a high degree of randomness, with very few (if any)
deterministic assignments, to minimize the predictability and the potential
for selection bias.

  * (b)

It should have a high degree of balance, i.e., at any allocation step, the
sample treatment numbers should be close to their targeted values:
$N_{k}(j)\approx j\rho_{k}$ for $j=1,\ldots,n$ and $k=1,\ldots,K$.

  * (c)

The unconditional probability of assigning treatment $k$ and the $j$th
allocation step should be equal to the target allocation proportion for this
treatment (allocation ratio preserving (ARP) property)(Kuznetsova and
Tymofyeyev, 2012):

$E(\phi_{jk})=\rho_{k}\;\;\textrm{for}\;j=1,\ldots,n\;\textrm{and}\;k=1,\ldots,K.$
(2)

###  2.1 Two-arm equal randomization

Consider a $1{:}1$ RCT to compare the effects of an experimental treatment (E)
vs. control (C). In this case, $\boldsymbol{\delta}_{j}=(1,0)^{T}$ if the
$j^{\text{th}}$ patient is assigned to E, $\boldsymbol{\delta}_{j}=(0,1)^{T}$
if the $j^{\text{th}}$ patient is assigned to C, and it is sufficient to
consider a single indicator variable: $\delta_{j}=1$ (or $0$), if the
$j^{\text{th}}$ patient’s assigned treatment is E (or C). Let
$N_{E}(j)=\sum_{i=1}^{j}\delta_{i}$ and $N_{C}(j)=j-N_{E}(j)$ denote the
sample sizes for treatments E and C, respectively, after $j$ allocation steps,
and let $D(j)=N_{E}(j)-N_{C}(j)$ denote the treatment imbalance after $j$
steps. Note that for every $j=1,\ldots,n$, $D(j)$ is a random variable whose
probability distribution is determined by the implemented randomization
procedure.

####  2.1.1 Randomization procedures for 1:1 RCTs

The simplest randomization procedure for a $1{:}1$ RCT is _complete
randomization_ : each subject’s treatment assignment is made independently of
other assignments by a flip of a fair coin, i.e., $\Pr(\delta_{j}=1)=0.5$ for
$j=1,\ldots,n$ (Lachin, 1988). While such an approach may be reasonable for
large sample sizes, it has a non-negligible probability of treatment imbalance
when the sample size is small.

In practice, most of $1{:}1$ RCTs utilize _restricted_ randomization (Pocock,
1979). For the $j$th patient, let
$\boldsymbol{\phi}_{j}=(\phi_{j},1-\phi_{j})^{T}$ be a vector of conditional
randomization probabilities for treatments E and C. It is sufficient to
consider only the first component, $\phi_{j}\in[0,1]$, and define a restricted
randomization procedure as (Rosenberger and Lachin, 2015)):

$\phi_{1}=0.5\;\;\textrm{and}\;\;\phi_{j}=\Pr(\delta_{j}=1|\delta_{1},\ldots,\delta_{j-1}),\;j\geq
2.$ (3)

Below we list the randomization procedures for a 1:1 RCT implemented in the
current version of Incertus.jl. These procedures can be grouped into several
classes depending on their probabilistic structure.111For a given procedure,
there may be parameter(s) that determine a tradeoff between treatment balance
and allocation randomness.

  * •

Complete Randomization Design (CRD): every treatment assignment is made with
probability 0.5, irrespective of the previous assignments (Lachin, 1988).

  * •

Procedures that force strict final balance ($n/2$ E’s and $n/2$ C’s): Random
Allocation Rule (Rand) and Truncated Binomial Design (TBD) (Blackwell and
Hodges, 1957).

  * •

A procedure that forces strict balance after $2b$ allocation steps, where $b$
is some small positive integer: Permuted Block Design (PBD) (Matts and Lachin,
1988).

  * •

Maximum Tolerated Imbalance (MTI) procedures that control imbalance within
pre-specified limits ($\pm b$) with fewer restrictions on randomization
compared to PBD: Big Stick Design (BSD) (Soares and Wu, 1983); Biased Coin
Design With Imbalance Tolerance (BCDWIT) Chen (1999); Ehrenfest Urn Design
(EUD) (Chen, 2000); Block Urn Design (BUD) (Zhao and Weng, 2011).

  * •

Biased Coin Designs (BCDs) with the allocation function of the form
$\phi_{j+1}=F(D(j))$, where $F:\mathbb{Z}\rightarrow[0,1]$ is a non-
increasing, symmetric around zero function defined on the set of integers such
that $F(D(j))+F(-D(j))=1$. One notable procedure of this class is Efron’s BCD
with $F_{p}(x)=0.5+sign(x)(0.5-p)$, where $sign(x)$ returns the value $-1$, 0,
1 if $x$ is negative, zero, or positive, respectively, and $p\in(0.5,1]$ is
the fixed coin bias probability (Efron, 1971). Another procedure of this class
is the Adjustable Biased Coin Design (ABCD) (Baldi Antognini and Giovagnoli,
2004) with the allocation function $F_{a}(x)=\\{1+x^{a}\\}^{-1}$, where $a\geq
0$ is a user-specified constant ($a=0$ corresponds to CRD and larger values of
$a$ correspond to more balanced and less random procedures).

  * •

Generalized Biased Coin Designs (GBCDs) (Wei, 1978; Smith, 1984) with the
allocation function of the form $\phi_{j+1}=f\left(\frac{D(j)}{j}\right)$,
where $f:[-1,1]\rightarrow[0,1]$ is a continuous, differentiable at 0, non-
increasing function such that $f(x)+f(-x)=1$. We consider a family of
functions proposed by Smith (1984):
$f_{\rho}(x)=\frac{(1-x)^{\rho}}{(1-x)^{\rho}+(1+x)^{\rho}}$, where $\rho\geq
0$ is a user-specified constant ($\rho=0$ corresponds to CRD and larger values
of $a$ correspond to more balanced and less random procedures).

  * •

Bayesian Biased Coin Designs (BBCDs) (Ball et al., 1993) with the allocation
rule

$\phi_{1}=0.5;\;\;\phi_{2}=\mathbb{1}\left\\{\delta_{1}=0\right\\};\;\;\phi_{j+1}=\frac{\left\\{1+\frac{N_{C}(j)}{jN_{E}(j)}\right\\}^{1/\gamma}}{\left\\{1+\frac{N_{C}(j)}{jN_{E}(j)}\right\\}^{1/\gamma}+\left\\{1+\frac{N_{E}(j)}{jN_{C}(j)}\right\\}^{1/\gamma}},\;j\geq
2.$ (4)

In (4), $\mathbb{1}\left\\{\delta_{1}=0\right\\}=1$ (or 0), if the first
treatment assignment was C (or E), and $\gamma>0$ is a user-specified constant
(smaller values of $\gamma$ correspond to more balanced and less random
procedures).

Note that for all aforementioned 1:1 randomization procedures, the conditional
randomization probability for an incoming patient can be analytically
expressed as a function of the current treatment sample sizes, i.e., we have
$\phi_{j+1}=\Pr(\delta_{j+1}=1|N_{E}(j),N_{C}(j))$. This is referred to as
_Markov property_ (Zhao, 2018). However, there are 1:1 randomization
procedures that do not have Markov property; e.g., the maximal procedure
(Berger et al., 2003). Such procedures require the full knowledge of the
completed allocation path (not just the numbers $N_{E}(j)$ and $N_{C}(j)$) to
determine the treatment assignment for an incoming patient. The implementation
of randomization procedures without Markov property is deferred to our future
work.

####  2.1.2 Statistical characteristics for 1:1 randomization procedures

For 1:1 restricted randomization procedures, a key consideration is a tradeoff
between treatment balance and allocation randomness (Sverdlov and Ryeznik,
2023). The probability distribution of imbalance $D(j)$ for $j=1,\ldots,n$ may
be an important starting point for obtaining various operating characteristics
of a randomization design. The following four measures of imbalance for
$j=1,\ldots,n$ are implemented in Incertus.jl:

  * •

$E|D(j)|$ – expected absolute imbalance at the $j$th step.

  * •

$E|D(j)|^{2}=var\\{D(j)\\}$ – variance of imbalance at the $j$th step.

  * •

$E\left\\{\max_{1\leq i\leq j}|D(j)|\right\\}$ – expected maximum imbalance
over first $j$ allocation steps.

  * •

A cumulative average loss (Atkinson, 2002) over $j$ steps:

$L(j)=\frac{1}{j}\sum_{i=1}^{j}\frac{E|D(i)|^{2}}{i}.$ (5)

As regards lack of randomness, the following four measures for $j=1,\ldots,n$
are implemented in Incertus.jl:

  * •

$EPCG_{conv}(j)$ – expected proportion of correct guesses over first $j$
allocation steps under the _convergence_ guessing strategy (Blackwell and
Hodges, 1957). Under this model, the investigator’s guess on an upcoming
treatment assignment, $\delta_{i+1}$, depends on the sign of the current
imbalance $D(i)$: guess the $(i+1)$st assignment as E, if $D(i)<0$; as C, if
$D(i)>0$; or make a random guess with probability 0.5, if $D(i)=0$.

  * •

$EPCG_{max}(j)$ – expected proportion of correct guesses over first $j$
allocation steps under the _maximum probability_ guessing strategy: guess the
$i$th treatment assignment as E, if $\phi_{i}>0.5$; as C, if $\phi_{i}<0.5$;
or guess at random, if $\phi_{i}=0.5$.

  * •

$EPDA(j)=\frac{1}{j}\sum_{i=1}^{j}\Pr(\phi_{i}=0\;\textrm{or}\;1)$ – expected
proportion of deterministic assignments over first $j$ allocation steps.

  * •

$FI(j)=\frac{1}{j}\sum_{i=1}^{j}E|\phi_{i}-0.5|$ – forcing index (Heritier et
al., 2005) that takes values on a scale 0–1. We have $FI(j)\equiv 0$ $\forall
j$ for CRD (most random design) and $FI(j)=1$ for the permuted blocks of 2,
assuming $j$ is even (most balanced design).

To quantify the balance–randomness tradeoff at the $j$th allocation step, we
consider the Euclidean distance to the origin $(0,0)$ in the Loss–Forcing
Index space (Sverdlov and Ryeznik, 2023):

$G(j)=\sqrt{\\{L(j)\\}^{2}+\\{FI(j)\\}^{2}},$ (6)

with lower values of $G(j)$ indicating better balance–randomness tradeoff. A
plot of $G(j)$ for $j=1,\ldots,n$ can provide valuable insights into the
dynamic behavior of balance–randomness tradeoff of a chosen randomization
procedure, and it can be useful for comparing different randomization
procedures.

###  2.2 Multi-arm randomization

A general restricted randomization procedure for a trial with $K\geq 2$
treatment arms and a fixed target allocation ratio
$\rho_{1}{:}\rho_{2}{:}\ldots{:}\rho_{K}$, where $\rho_{i}\in(0,1)$ and
$\sum_{i=1}^{K}\rho_{i}=1$ is defined by Equation (1). Two important data
structures in this context are the $n\times K$ matrix of treatment assignments
$\boldsymbol{\Delta}_{n}=(\boldsymbol{\delta}_{1},\boldsymbol{\delta}_{2},\ldots,\boldsymbol{\delta}_{n})^{T}$
and the $n\times K$ matrix of conditional randomization probabilities
$\boldsymbol{\Phi}_{n}=(\boldsymbol{\phi}_{1},\boldsymbol{\phi}_{2},\ldots,\boldsymbol{\phi}_{n})^{T}$
. In the current version of Incertus.jl, we mainly consider randomization
procedures with Markov property, i.e., for which the conditional randomization
probabilities for an incoming patient depend on the past allocations only
through the current treatment numbers, in other words:

$\phi_{j+1,k}=\Pr(\boldsymbol{\delta}_{j+1}=\boldsymbol{e}_{k}|\boldsymbol{N}(j)),\;\;j\geq
1,\;k=1,\ldots,K,$ (7)

where $\boldsymbol{N}(j)=(N_{1}(j),\ldots,N_{K}(j))^{T}$ and
$N_{k}(j)=\sum_{i=1}^{j}\mathbb{1}\\{\boldsymbol{\delta}_{i}=\boldsymbol{e}_{k}\\}$
is the number of patients assigned to the $k$th treatment after $j$ allocation
steps, with $\sum_{k=1}^{K}N_{k}(j)=j$.

While most of the multi-arm restricted randomization procedures currently
implemented in Incertus.jl have the allocation rules of the form (7), some do
not (e.g., the drop-the-loser urn design (Ivanova, 2003) for which the
conditional randomization probabilities are determined by the current urn
composition rather than by the treatment sample sizes $\boldsymbol{N}(j)$). In
general, procedures without Markov property require some customized
programming, and we defer the broad implementation of such procedures to
future work.

####  2.2.1 Randomization procedures for multi-arm RCTs

Below we list the randomization procedures for a $K\geq 2$-arm trial
implemented in the current version of Incertus.jl. Note that for a special
case of $K=2$ and $\rho_{1}=\rho_{2}=0.5$, a multi-arm randomization procedure
may take the form of some 1:1 randomization procedure described in Section
2.1.

  * •

Complete Randomization Design (CRD) applies multinomial probabilities equal to
the target allocation proportions to determine the treatment assignment at
each allocation steps: $\phi_{jk}=\rho_{k}$ for $j=1,\ldots,n$ and
$k=1,\ldots,K$.

  * •

Procedures that force strict balance for the given sample size: Suppose the
target sample sizes (positive integers) for the $K$ treatment arms are
$n_{1}:n_{2}\ldots:n_{K}$ with $\sum_{k=1}^{K}n_{k}=n$ (here
$n_{k}=n\rho_{k}$, rounded to integer if necessary). Two designs that can be
used to achieve this allocation exactly are the Random Allocation Rule (Rand)
and the Truncated Multinomial Design (TMD) (Rukhin, 2007).

  * •

Permuted Block Design (PBD): Consider the target allocation ratio
$w_{1}:w_{2}:\ldots:w_{K}$, where $w_{k}$’s are positive integers with the
common divisor of 1. Let $W=\sum_{k=1}^{K}w_{k}$ the size of the “minimal
balanced” set. The PBD performs treatment assignments in blocks of size $bW$,
where $b$ is some small positive integer. Within each block, exactly $bw_{k}$
patients are randomized to the $k$th treatment, $k=1,\ldots,K$.

  * •

Randomization procedures that can achieve (approximately) any fixed integer-
valued allocation ratio $w_{1}:w_{2}:\ldots:w_{K}$, where $w_{k}$’s are
positive integers with the common divisor of 1, with fewer restrictions than
PBD: Block Urn Design (BUD) (Zhao and Weng, 2011) and Drop-the-Loser Urn
Design (DLUD) (Ivanova, 2003; Ryeznik and Sverdlov, 2018).

  * •

Randomization procedures that can achieve (approximately) any fixed, possibly
irrational-valued allocation proportions
$\rho_{1}{:}\rho_{2}{:}\ldots{:}\rho_{K}$, $\rho_{k}\in(0,1)$ and
$\sum_{k=1}^{K}\rho_{k}=1$ (in particular, they can be applied with the
integer-valued ratios $w_{1}:w_{2}:\ldots:w_{K}$). Three designs that are fit
for this purpose are Mass Weighted Urn Design (MWUD) Zhao (2015), Doubly
Adaptive Biased Coin Design (DBCD) Hu and Zhang (2004), and Maximum Entropy
Constrained Balance Randomization (MaxEnt) (Klotz, 1978; Ryeznik and Sverdlov,
2018).

####  2.2.2 Statistical characteristics for multi-arm randomization procedures

For a $K\geq 2$-arm trial, by balanced allocation, we mean the achievement of
the target allocation ratio for the given sample size. A natural measure of
lack of balance after $j$ allocation steps is the Euclidean distance between
the observed allocation $\boldsymbol{N}(j)=(N_{1}(j),\ldots,N_{K}(j))^{T}$ and
the targeted allocation $j\boldsymbol{\rho}=(j\rho_{1},\ldots,j\rho_{K})^{T}$:

$d(j)=\sqrt{\sum_{i=1}^{K}(N_{i}(j)-j\rho_{i})^{2}}.$ (8)

Since $N_{i}(j)$’s are random, $d(j)$ in (8) is a random variable, with
smaller values indicating more balanced allocation. For $K=2$ and
$\rho_{1}=\rho_{2}=0.5$, (8) is equal to
$\frac{1}{\sqrt{2}}|N_{1}(j)-N_{2}(j)|=\frac{1}{\sqrt{2}}|D(j)|$. The
following measures of lack of balance for $j=1,\ldots,n$ are implemented in
Incertus.jl:

  * •

$E[d(j)]$ – expected distance between the observed and the targeted treatment
sample sizes at the $j$th step.

  * •

$E[d^{2}(j)]=\sum_{i=1}^{K}var[N_{i}(j)]$ – total variance of the treatment
sample sizes at the $j$th step.

  * •

$E\left\\{\max_{1\leq i\leq j}d(i)\right\\}$ – expected maximum distance
between the observed and the targeted treatment sample sizes over first $j$
allocation steps.

  * •

A cumulative loss over $j$ allocation steps:

$L(j)=\frac{1}{j}\sum_{i=1}^{j}\frac{E[d^{2}(i)]}{i}.$ (9)

To quantify the lack of randomness, we consider the following measures:

  * •

$EPCG_{min}(j)$ – expected proportion of correct guesses over first $j$
allocation steps under the _minimum imbalance_ guessing strategy (Zhao et al.,
2017). Under this model, let $\Delta_{k}(i)=N_{k}(i)-i\rho_{k}$ be the
difference between the achieved sample size and the “ideal” sample size for
the $k$th treatment ($k=1,\ldots,K$) after $i$ allocation steps. Then the
investigator guesses the $(i+1)$st assignment to be treatment $s$ that is
currently most underrepresented, i.e., for which $\Delta_{s}(i)=\min_{1\leq
k\leq K}\Delta_{k}(i)$ (if several treatments meet this condition, a random
guess is made on this set of treatments).

  * •

$EPCG_{max}(j)$ –expected proportion of correct guesses over first $j$
allocation steps under the _maximum probability_ guessing strategy (Zhao et
al., 2017): guess the $i$th treatment assignment to be treatment $s$ that has
the maximum conditional randomization probability, i.e.,
$\phi_{is}=\max_{1\leq k\leq K}\phi_{ik}$ (if several treatments meet this
condition, a random guess is made on this set of treatments).

  * •

$EPDA(j)$ – expected proportion of deterministic assignments over first $j$
allocation steps. Here a “deterministic assignment” at the $i$th allocation
step means that some treatment is assigned with probability 1 at the $i$th
step. Therefore,

$EPDA(j)=\frac{1}{j}E\left[\sum_{i=1}^{j}\mathbb{1}\left\\{\exists
k:\phi_{ik}=1\;\text{and}\;\forall k^{\prime}\neq
k:\phi_{ik^{\prime}}=0\right\\}\right].$ (10)

  * •

Forcing Index (FI) defined as the cumulative average (over $j$ allocation
steps) of the expected Euclidean distance from the vector of treatment
randomization probabilities and the target allocation proportions:

$FI(j)=\frac{1}{j}\sum_{i=1}^{j}E\sqrt{\sum_{k=1}^{K}(\phi_{ik}-\rho_{k})^{2}}.$
(11)

Small values of $FI(j)$ correspond to more random procedures. For CRD,
$FI(j)\equiv 0$ $\forall j$.

To quantify balance–randomness tradeoff for a given randomization design,
suppose we have $U_{B}(j)$ and $U_{R}(j)$ – the measures of lack of balance
and lack of randomness after $j$ allocation steps, normalized to the same
scale of 0–1, such that smaller values of $U_{B}(j)$ indicate better balance,
and smaller values of $U_{R}(j)$ indicate better randomness (Sverdlov and
Ryeznik, 2019). Then, assuming balance and randomness are equally important,
we can define an “overall” performance measure as the distance to the origin
(0,0):

$G(j)=\sqrt{\\{U_{B}(j)\\}^{2}+\\{U_{R}(j)\\}^{2}},$ (12)

with lower values of $G(j)$ indicating a better balance–randomness tradeoff.

Finally, to assess the ARP property of a randomization procedure, we obtain
the unconditional treatment randomization probabilities at every step. These
probabilities can be estimated through Monte Carlo simulations:

$\pi_{jk}=E(\phi_{jk})\approx\frac{1}{N_{sim}}\sum_{\ell=1}^{N_{sim}}\phi_{jk}^{\ell},$
(13)

where $\phi_{jk}^{\ell}$ is the conditional randomization probability for the
$k$th treatment at the $j$th allocation step from the $\ell$th simulation run
($j=1,\ldots,n$, $k=1,\ldots,K$, $\ell=1,\ldots,N_{sim}$). A plot of
$\pi_{jk}$ vs. allocation step can be very useful – for a procedure that lacks
ARP property one may expect deviations from the target allocation proportions
$\rho_{k}$’s and fluctuations in $\pi_{jk}$ as $j$ varies from 1 to $n$.

##  3 Incertus.jl Package

###  3.1 Overview

Incertus.jl is available at the GitHub repository:
https://github.com/yevgenryeznik/Incertus.jl. To install the package, one
should run the command in Julia REPL:

] add "https://github.com/yevgenryeznik/Incertus.jl"

The package documentation is available at
https://yevgenryeznik.github.io/Incertus.jl/build/.

Tables 1–2 summarize the functions and methods available in Incertus.jl.

Table 1: Functions for initializing different randomization procedures in
Incertus.jl package Functions for initializing 1:1 randomization procedures
---
CRD() | initializes a complete randomization design
RAND(n) | initializes a random allocation rule; sample size equals to $n$
TBD(n) | initializes a truncated binomial design; sample size equals to $n$
PBD(b) | initializes a permuted block design with a block size equal to $2b$
BSD(b) | initializes a big stick design with MTI parameter equal to $b$
BCDWIT(p,b) | initializes biased coin design with imbalance tolerance; biased coin probability = $p$; MTI parameter = $b$
EUD(b) | initializes Ehrenfest urn design with MTI parameter equal to $b$
BUD($\lambda$) | initializes block urn design with a parameter equal to $\lambda$
EBCD(p) | initializes Efron’s biased coin design; biased coin probability = $p$
ABCD(p) | initializes adjustable biased coin design with a parameter equal to $a$
GBCD($\gamma$) | initializes generalized biased coin design with a parameter equal to $\gamma$
BBCD($\gamma$) | initializes the Bayesian biased coin design with a parameter equal to $\gamma$
Functions for initializing multi-arm randomization procedures
(an input parameter vector $\boldsymbol{w}$ specifies target allocation)
CRD(w) | initializes a complete randomization design
RAND(w,n) | initializes a random allocation rule; sample size equals to $n$
TBD(w,n) | initializes a truncated multinomial design; sample size equals to $n$
PBD(w,b) | initializes a permuted block design; block size equals to $b(w_{1}+\ldots+w_{K})$
BUD(w,$\lambda$) | initializes block urn design with a parameter equal to $\lambda$
MWUD(w,$\alpha$) | initializes mass weighted urn design with a parameter equal to $\alpha$
DLUD(w,a) | initializes drop-the-loser urn design with a parameter equal to $a$
DBCD(w,$\gamma$) | initializes doubly adaptive biased coin design with a parameter equal to $\gamma$
MaxEnt(w,$\eta$) | initializes MaxEnt design with a parameter equal to $\eta$
Table 2: Methods implemented in Incertus.jl package Simulation of
randomization procedures
---
simulate(rnd,nsbj,nsim,seed) | simulates a randomization. Here rnd represents either a single randomization
| procedure or a set of procedures to be simulated; nsbj is a sample size;
|  nsim is the number of simulations; and seed is used for initializing a
random
| number generator (a default is set to 314159)
Calculation of operating characteristics
(an input parameter sr is the output of the simulate function, representing
the simulation results)
Balance
calc_final_imb(sr) | calculates the distribution of the final imbalance
calc_expected_abs_imb(sr) | calculates expected absolute imbalance vs. allocation step
calc_variance_of_imb(sr) | calculates variance of imbalance vs. allocation step
calc_expected_max_abs_imb(sr) | calculates expected maximum absolute imbalance vs. allocation step
calc_cummean_loss(sr) | calculates cumulative average loss vs. allocation step
Randomness
calc_cummean_epcg(sr,gs) | calculates cumulative averages of expected proportions of correct guesses vs.
| allocation step; gs parameter accepts two values: “C” (converging guessing
strategy)
| or “MP” (maximum probability guessing strategy)
calc_cummean_pda(sr) | calculates cumulative averages of the proportions of deterministic assignments vs.
| allocation step
calc_fi(sr) | calculates forcing index vs. allocation step
Balance-randomness tradeoff and the allocation ratio preserving (ARP) property
calc_brt(sr) | calculates balance-randomness tradeoff vs. allocation step
eval_arp(sr) | evaluates ARP property, i.e., calculates unconditional allocation probabilities
Visualization of operating characteristics
(op is a data frame with corresponding operating characteristics; brt is an
output of calc_brt function;
arp is an output of eval_arp function; final_imb is an output of
calc_final_imb function;
kwargs is a set of key words that can be accepted by the visualization
functions
violin(final_imb; kwargs...) | visualize a distribution of the final imbalance as a violin plot
heatmap(brt; kwargs...) | visualize balance-randomness trade-off as a heatmap plot
plot(op; kwargs...) | visualize simulated operating characteristics
plot(arp; kwargs...) | visualize simulated unconditional allocation probabilities

###  3.2 Examples: Two-arm equal randomization

####  3.2.1 Generating a randomization sequence for a 1:1 RCT

The function simulate is a core function for sequentially generating treatment
randomization probabilities and treatment assignments for the chosen
randomization method(s). Suppose we want to generate a randomization sequence
for a 1:1 RCT with $n=8$. We consider four randomization procedures:

  * •

PBD with blocks of size $2b=4$ – PBD(2)

  * •

Random allocation rule – RAND(8)

  * •

BSD with $b=3$ – BSD(3)

  * •

Efron’s BCD with $p=2/3$ – EBCD(2/3)

The following code is used:

using Incertus

# sample size
nsbj = 8;

# number of simulations is set to 1
nsim = 1;

# randomization procedures to be simulated
rnd = [PBD(2), RAND(nsbj), BSD(3), EBCD(2//3)];

# labels
label_rnd = [label(item) for item in rnd]

# simulation run
sr = simulate(rnd, nsbj, nsim);

The function simulate takes as the input the set of four specified
randomization procedures (rnd = [PBD(2), RAND(nsbj), BSD(3), EBCD(2//3)]), the
sample size (nsbj = 8), and the number of simulations (nsim = 1, which
indicates that we are interested in simulating a single sequence for each
randomization method). Also, it is possible to provide a seed for a random
number generator as a fourth input argument, which is set to 314159 by
default. The output object sr contains four data structures corresponding to
four implemented randomization procedures. For example, sr[1] has the output
for PBD(2), sr[2] has the output for RAND, etc. Then the $8\times 2$ matrix of
treatment assignments and the $8\times 2$ matrix of conditional randomization
probabilities for PBD(2) can be obtained as sr[1].trt and sr[1].prb,
respectively. Similar commands can be used to extract the corresponding data
structures for the other three randomization procedures.

We can also create data frames combining treatment E assignments (trt_df) and
conditional randomization probabilities for treatment E (prb_df) of the four
considered randomization procedures using the following code:

using DataFrames

trt = hcat([sr[i].trt[:, 1, 1] for i in eachindex(sr)]...);
trt_df = DataFrame(trt, label_rnd)

prb = hcat([sr[i].prb[:, 1, 1] for i in eachindex(sr)]...);
prb_df = DataFrame(prb, label_rnd)

The resulting data frames are displayed in Figure 1. One can see that the four
randomization procedures have different vectors of conditional randomization
probabilities (displayed as columns of the right table) and result in
different treatment allocation sequences (displayed as columns of the left
table). See Supplemental Appendix, “Simulation example 1” for details.

Figure 1: Data frames of treatment assignments (left) and conditional
randomization probabilities (right) of four randomization designs – PBD(2),
RAND, BSD(3), and EBCD(2/3) – for a 1:1 RCT with $n=8$.

####  3.2.2 Assessing statistical properties of 1:1 randomization designs
through Monte Carlo simulation

Monte Carlo simulations provide a useful tool for evaluating various
randomization designs and selecting one that is most fit for purpose for the
given trial parameters. Suppose we want to compare seven randomization designs
for a 1:1 RCT with $n=40$ patients with respect to the criteria of balance and
randomness, as described in Section 2.1.2. The designs are:

  * •

Complete randomization design – CRD

  * •

PBD with blocks of size $2b=2$ – PBD(1)

  * •

Random allocation rule – RAND

  * •

Truncated binomial design – TBD

  * •

Efron’s BCD with $p=2/3$ – EBCD(2/3)

  * •

Adjustable BCD with $a=2$ – ABCD(2)

A Monte Carlo simulation (10,000 runs) is performed using the code:

using Incertus

# sample size
nsbj = 40;

# number of simulations is set to 1
nsim = 10000;

# randomization procedures to be simulated
rnd = [CRD(), PBD(1), RAND(nsbj), TBD(nsbj), BSD(3), EBCD(2//3), ABCD(2);

# simulation run
sr = simulate(rnd, nsbj, nsim);

Based on the output object, sr, the selected design operating characteristics
can be calculated and visualized using Julia graphical tools using the
following code:

# calculating cumulative average loss vs. allocation step
cummean_loss = calc_cummean_loss(sr);
# creating a plot of cumulative average losses over the first allocation steps
plot(cummean_loss, ylabel = "cumulative average loss", size = (800, 600))

# calculating forcing index vs. allocation step
fi = calc_fi(sr);
# creating a plot of forcing indices vs. allocation step
plot(fi, ylabel = "forcing index", size = (800, 600))

# Evaluating balance-randomness tradeoff vs. allocation step
brt = calc_brt(sr);
# creating a heatmap plot of balance--randomness tradeoff vs. allocation step
heatmap(brt, size = (800, 600))

# calculating final imbalance
final_imb = calc_final_imb(sr);
# creating a violin plot of final imbalances
violin(final_imb, size = (800, 600))

Figure 2(a–d) displays, respectively, the plots of cumulative average loss,
forcing index, balance-randomness tradeoff heatmap plot, and the violin plot
of final imbalance. We can see that among the seven considered designs in this
scenario, BSD(3) provides, overall, the best tradeoff between balance and
randomness. The full report of this simulation is available in the
Supplemental Appendix, “Simulation example 2”.

Figure 2: Selected operating characteristics of seven randomization designs –
CRD, PBD(1), RAND, TBD, BSD(3), EBCD(2/3), and ABCD(2) – for a 1:1 RCT with
$n=40$.

####  3.2.3 Calibrating 1:1 randomization design through Monte Carlo
simulation

Most restricted randomization designs involve a parameter that can be fine-
tuned to determine a design with a desirable balance-randomness tradeoff. For
example, consider the Bayesian Biased Coin Design (BBCD) with five different
values of the parameter $\gamma=0.01;0.05;0.1;0.2;$ and 1, for a 1:1 RCT with
$n=40$. We ran 10,000 Monte Carlo simulations to obtain the operating
characteristics of these designs. The Julia code and the full report of
simulation results are available in the Supplemental Appendix, “Simulation
example 3". Figure 3 below shows a heatmap plot of the balance–randomness
tradeoff metric. One can see that $\gamma=0.05$ is an “optimal” choice – it
results in increasingly smaller values of the metric over allocation steps,
the smallest among the five considered designs.

Figure 3: Heatmap plot of balance-randomness tradeoff for the Bayesian biased
coin design (BBCD) with five different values of the parameter
$\gamma=0.01;0.05;0.1;0.2;$ and 1, for a 1:1 RCT with $n=40$.

###  3.3 Examples: Multi-arm equal randomization

####  3.3.1 Generating a randomization sequence for a 4:3:2:1 RCT

Here, we demonstrate how to generate a randomization sequence and
randomization probabilities for a multi-arm RCT, targeting 4:3:2:1 allocation,
with $n=10$. We consider two randomization designs: random allocation rule
(RAND) and truncated multinomial design (TMD). The following code is used:

using Incertus

# sample size
nsbj = 10;

# number of simulations is set to 1
nsim = 1;

# target allocation
w = [4, 3, 2, 1];

# randomization procedures to be simulated
rnd = [RAND(w, nsbj), TMD(w, nsbj)];

# labels
label_rnd = [label(item) for item in rnd]

# simulation run
sr = simulate(rnd, nsbj, nsim);

The output object sr contains two data structures: sr[1] for RAND and sr[2]
for TMD. For RAND, the $10\times 4$ matrix of treatment assignments is
retrieved as sr[1].trt, and the $10\times 4$ matrix of conditional
randomization probabilities is retrieved as sr[1].prb. Likewise, for TMD, the
corresponding matrices are retrieved as sr[2].trt and sr[2].prb. See
Supplemental Appendix, “Simulation example 5” for details.

####  3.3.2 Assessing statistical properties of 4:3:2:1 randomization designs
through Monte Carlo simulation

Consider a four-arm RCT with $n=40$ patients, targeting 4:3:2:1 allocation. We
are interested in comparing seven randomization designs with respect to the
criteria of balance and randomness, as described in Section 2.2.2. The designs
are:

  * •

Complete randomization design – CRD

  * •

PBD with a block size $bW=1\cdot(4+3+2+1)=10$ – PBD(1)

  * •

Block urn design with $\lambda=2$ – BUD(2)

  * •

Random allocation rule – RAND

  * •

Truncated multinomial design – TMD

  * •

Drop-the-loser urn design with $a=2$ – DLUD(2)

  * •

Mass weighted urn design with $\alpha=2$ – MWUD(2)

A Monte Carlo simulation (10,000 runs) is performed using the code:

using Incertus

# sample size
nsbj = 40;

# number of simulations
nsim = 10000;

# target allocation
w = [4, 3, 2, 1];

# randomization procedures to be simulated
rnd = [CRD(w), PBD(w, 1), BUD(w, 2), RAND(w, nsbj), TMD(w, nsbj), DLUD(w, 2),
MWUD(w, 2)];

# simulation run
sr = simulate(rnd, nsbj, nsim);

The design operating characteristics are calculated based on the output object
sr, and they are visualized using Julia graphical tools; see Supplemental
Appendix, “Simulation example 6”.

####  3.3.3 Calibrating a 4:3:2:1 randomization design through Monte Carlo
simulation

Suppose an experimenter decides to use a doubly adaptive biased coin design
(DBCD) for a 4:3:2:1 RCT with $n=40$ patients. Our goal is to find a value of
the parameter $\gamma>0$ that would result in a procedure with a sensible
balance-randomness tradeoff. For illustrative purposes, let us consider DBCD
with $\gamma=0.01;1;2;5;$ and 10. We ran 10,000 Monte Carlo simulations to
obtain the operating characteristics of these designs. The Julia code and the
full report of simulation results are available in the Supplemental Appendix,
“Simulation example 7”. Based on these simulations, Figure 4 shows the violin
plot of the final imbalance (left plot) and the heatmap plot of the
balance–randomness tradeoff metric (right plot). Figure 5 shows the plot of
simulated unconditional randomization probabilities for the four treatment
arms. From Figure 5, the DBCD does not have the ARP property because there are
some fluctuations (deviations from the target allocation proportions) in the
unconditional probabilities; these fluctuations are negligible for
$\gamma=0.1$, and they are more pronounced for larger values of $\gamma$.

Figure 4: (a) Violin plot of the final imbalance; and (b) Heatmap plot of
balance–randomness tradeoff – for DBCD with five different values of the
parameter $\gamma=0.01;1;2;5;$ and 10 – for a 4:3:2:1 RCT with $n=40$. Figure
5: Simulated unconditional randomization probabilities for the four treatment
arms – for DBCD with five different values of the parameter
$\gamma=0.01;1;2;5;$ and 10 – for a 4:3:2:1 RCT with $n=40$.

##  4 Discussion and future work

Incertus.jl is a Julia software package for generating randomization sequences
for clinical trials with two or more treatment arms, with equal or unequal
allocation ratios. The package is highly efficient computationally and it can
be invoked in R. It has a utility for performing simulation studies of
selected randomization procedures and assessing their statistical properties,
which could help the investigators select a randomization design that fits for
purpose of their study. While our examples in the current paper concerned
relatively simple cases, the software provides “building blocks” for
developing randomization designs for more complex studies – e.g., trials with
stratified randomization (Sverdlov et al., 2024), platform trials (The
Adaptive Platform Trials Coalition, 2019), etc.

It is worth highlighting some opportunities for extension and future work.
First, adding more randomization procedures to the package would be useful.
The procedures with Markov property (whose conditional randomization
probabilities have analytical expressions as functions of current treatment
sample sizes) can be easily programmed following the existing template in
Incertus.jl. For the procedures that do not have Markov property, such as the
maximal procedure of Berger et al. (2003), some more customized approach to
programming will be warranted.

Second, it will be important to enrich the package with some functions that
would enable characterizing reference sets of various randomization
procedures, i.e., the sets of all possible randomization sequences and their
probabilities. Such an option is currently available in the randomizeR package
(Uschner et al., 2018) for 1:1 RCTs, and it would be also useful to implement
for more general settings, such as multi-arm RCTs with equal or unequal
allocation ratios, which are within the scope of Incertus.jl. Having a tool
that enables efficient sampling of randomization sequences from the reference
set of a given randomization procedure would allow investigators to construct
randomization-based tests (Rosenberger et al., 2019; Proschan and Dodd, 2019)
that complement standard likelihood-based tests.

Finally, would be beneficial to develop a user-friendly interface for
Incertus.jl. Since the package was developed in Julia, a natural way is to use
Julia functionality for building such an interface. At the same time, many
recent statistical software applications are developed as R Shiny tools.
Finding a way to cast Incertus.jl as an R Shiny tool would be an interesting
future research project.

## References

  * Wittes [2005] Janet Wittes.  Randomized treatment assignment.  In _Encyclopedia of Biostatistics_. John Wiley & Sons, Ltd, 2005.  doi:10.1002/0470011815.b2a01050.
  * Rosenberger and Lachin [2015] William F. Rosenberger and John M. Lachin.  _Randomization in Clinical Trials: Theory and Practice_.  Wiley, New York, 2nd edition, 2015.
  * Green [2002] Sylvan B. Green.  Design of randomized trials.  _Epidemiologic Reviews_ , 24(1):4–11, 2002.  doi:10.1093/epirev/24.1.4.
  * Friedman et al. [2015] Lawrence M. Friedman, Curt D. Furberg, David L. DeMets, David M. Reboussin, and Christopher B. Granger.  _Fundamentals of Clinical Trials_.  Springer International Publishing Switzerland, 5th edition, 2015.
  * Bretz et al. [2008] Frank Bretz, Jason Hsu, José Pinheiro, and Yi Liu.  Dose finding – a challenge in statistics.  _Biometrical Journal_ , 50(4):480–504, 2008.  doi:10.1002/bimj.200810438.
  * Schmidli et al. [2006] Heinz Schmidli, Frank Bretz, Amy Racine, and Willi Maurer.  Confirmatory seamless phase II/III clinical trials with hypotheses selection at interim: applications and practical considerations.  _Biometrical Journal_ , 48(4):635–643, 2006.  doi:10.1002/bimj.200510231.
  * The Adaptive Platform Trials Coalition [2019] The Adaptive Platform Trials Coalition.  Adaptive platform trials: definition, design, conduct and reporting considerations.  _Nature Reviews Drug Discovery_ , 18(10):797–807, 2019.  doi:s41573-019-0034-3.
  * Hilgers et al. [2017] Ralf-Dieter Hilgers, Diane Uschner, William F. Rosenberger, and Nicole Heussen.  ERDO – a framework to select an appropriate randomization procedure for clinical trials.  _BMC Medical Research Methodology_ , 17:159:1–12, 2017.  doi:10.1186/s12874-017-0428-z.
  * Berger et al. [2021] Vance W. Berger, Louis Joseph Bour, Kerstine Carter, Jonathan J. Chipman, Colin C. Everett, Nicole Heussen, Catherine Hewitt, Ralf-Dieter Hilgers, Yuqun Abigail Luo, Jone Renteria, Yevgen Ryeznik, Oleksandr Sverdlov, and Diane Uschner.  A roadmap to using randomization in clinical trials.  _BMC Medical Research Methodology_ , 21:168:1–24, 2021.  doi:10.1186/s12874-021-01303-z.
  * Sverdlov et al. [2024] Oleksandr Sverdlov, Yevgen Ryeznik, Volodymyr Anisimov, Olga M. Kuznetsova, Ruth Knight, Kerstine Carter, Sonja Drescher, and Wenle Zhao.  Selecting a randomization method for a multi-center clinical trial with stochastic recruitment considerations.  _BMC Medical Research Methodology_ , 24:52:1–23, 2024.  doi:10.1186/s12874-023-02131-z.
  * Uschner et al. [2018] Diane Uschner, David Schindler, Ralf-Dieter Hilgers, and Nicole Heussen.  randomizeR: An R package for the assessment and implementation of randomization in clinical trials.  _Journal of Statistical Software_ , 85(8):1–22, 2018.  doi:10.18637/jss.v085.i08.
  * Ryeznik and Sverdlov [2018] Yevgen Ryeznik and Oleksandr Sverdlov.  A comparative study of restricted randomization procedures for multiarm trials with equal or unequal treatment allocation ratios.  _Statistics in Medicine_ , 37:3056–3077, 2018.  doi:10.1002/sim.7817.
  * Ryeznik et al. [2018a] Yevgen Ryeznik, Oleksandr Sverdlov, and Andrew C. Hooker.  Adaptive optimal designs for dose-finding studies with time to-event outcomes.  _The AAPS Journal_ , 20:24:1–12, 2018a.  doi:10.1208/s12248-017-0166-5.
  * Ryeznik et al. [2018b] Yevgen Ryeznik, Oleksandr Sverdlov, and Andrew C. Hooker.  Implementing optimal designs for dose–response studies through adaptive randomization for a small population group.  _The AAPS Journal_ , 20:85:1–15, 2018b.  doi:10.1208/s12248-018-0242-5.
  * Fleiss [1986] Joseph L. Fleiss.  _The Design and Analysis of Clinical Experiments_.  Hoboken, NJ: John Wiley and Sons, 1986.  doi:10.1002/9781118032923.
  * Kuznetsova and Tymofyeyev [2012] Olga M. Kuznetsova and Yevgen Tymofyeyev.  Preserving the allocation ratio at every allocation with biased coin randomization and minimization in studies with unequal allocation.  _Statistics in Medicine_ , 31:701–723, 2012.  doi:10.1002/sim.4447.
  * Lachin [1988] John M. Lachin.  Properties of simple randomization in clinical trials.  _Controlled Clinical Trials_ , 9:312–326, 1988.  doi:10.1016/0197-2456(88)90046-3.
  * Pocock [1979] Stuart J. Pocock.  Allocation of patients to treatment in clinical trials.  _Biometrics_ , 35(1):183–197, 1979.
  * Blackwell and Hodges [1957] David Blackwell and J. L. Hodges.  Design for the control of selection bias.  _Annals of Mathematical Statistics_ , 28(2):449–460, 1957.  doi:10.1214/aoms/1177706973.
  * Matts and Lachin [1988] John P. Matts and John M. Lachin.  Properties of permuted-block randomization in clinical trials.  _Controlled Clinical Trials_ , 9:327–344, 1988.  doi:10.1016/0197-2456(88)90047-5.
  * Soares and Wu [1983] Jose F. Soares and C.F. Jeff Wu.  Some restricted randomization rules in sequential designs.  _Communications in Statistics—Theory and Methods_ , 12(17):2017–2034, 1983.  doi:10.1080/03610928308828586.
  * Chen [1999] Yung-Pin Chen.  Biased coin design with imbalance tolerance.  _Communications in Statistics—Stochastic Models_ , 15(5):953–975, 1999.  doi:10.1080/15326349908807570.
  * Chen [2000] Yung-Pin Chen.  Which design is better? Ehrenfest urn versus biased coin.  _Advances in Applied Probability_ , 32:738–749, 2000.  doi:10.1239/aap/1013540242.
  * Zhao and Weng [2011] Wenle Zhao and Yanqiu Weng.  Block urn design—a new randomization algorithm for sequential trials with two or more treatments and balanced or unbalanced allocation.  _Contemporary Clinical Trials_ , 32:953–961, 2011.  doi:10.1016/j.cct.2011.08.004.
  * Efron [1971] Bradley Efron.  Forcing a sequential experiment to be balanced.  _Biometrika_ , 58(3):403–417, 1971.  doi:10.1093/biomet/58.3.403.
  * Baldi Antognini and Giovagnoli [2004] Alessandro Baldi Antognini and Alessandra Giovagnoli.  A new ‘biased coin design’ for the sequential allocation of two treatments.  _Applied Statistics_ , 53(4):651–664, 2004.  doi:10.1111/j.1467- 9876.2004.00436.x.
  * Wei [1978] Lee-Jen Wei.  The adaptive biased coin design for sequential experiments.  _The Annals of Statistics_ , 6(1):92–100, 1978.  doi:10.1214/aos/1176344068.
  * Smith [1984] Richard L. Smith.  Sequential treatment allocation using biased coin designs.  _Journal of the Royal Statistical Society. Series B_ , 46(3):519–543, 1984.  doi:10.1111/j.2517-6161.1984.tb01323.x.
  * Ball et al. [1993] Frank G. Ball, Adrian F.M. Smith, and Isabella Verdinelli.  Biased coin designs with a bayesian bias.  _Journal of Statistical Planning and Inference_ , 34(3):403–421, 1993.  doi:10.1016/0378-3758(93)90148-Y.
  * Zhao [2018] Wenle Zhao.  Evolution of restricted randomization with maximum tolerated imbalance.  In Vance W. Berger, editor, _Randomization, Masking, and Allocation Concealment_ , pages 61–81. CRC Press, Boca Raton FL, 2018.  doi:10.1201/9781315305110.
  * Berger et al. [2003] Vance W. Berger, Anastasia Ivanova, and Maria Deloria Knoll.  Minimizing predictability while retaining balance through the use of less restrictive randomization procedures.  _Statistics in Medicine_ , 22:3017–3028, 2003.  doi:10.1002/sim.1538.
  * Sverdlov and Ryeznik [2023] Oleksandr Sverdlov and Yevgen Ryeznik.  Balancing the objectives of statistical efficiency and allocation randomness in randomized controlled trials.  _Statistics in Biopharmaceutical Research_ , 2023.  doi:10.1080/19466315.2023.2261671.
  * Atkinson [2002] Anthony C. Atkinson.  The comparison of designs for sequential clinical trials with covariate information.  _Journal of the Royal Statistical Society. Series A_ , 165(2):349–373, 2002.  doi:10.1111/1467-985X.00564.
  * Heritier et al. [2005] Stephane Heritier, Val Gebski, and Avinesh Pillai.  Dynamic balancing randomization in controlled clinical trials.  _Statistics in Medicine_ , 24:3729–3741, 2005.  doi:10.1002/sim.2421.
  * Ivanova [2003] Anastasia Ivanova.  A play-the-winner-type urn design with reduced variability.  _Metrika_ , 58(1):1–13, 2003.  doi:10.1007/s001840200220.
  * Rukhin [2007] Andrew L. Rukhin.  Normal order statistics and sums of geometric random variables in treatment allocation problems.  _Statistics and Probability Letters_ , 77:1312–1321, 2007.  doi:10.1016/j.spl.2007.03.017.
  * Zhao [2015] Wenle Zhao.  Mass weighted urn design—a new randomization algorithm for unequal allocations.  _Contemporary Clinical Trials_ , 43:209–216, 2015.  doi:10.1016/j.cct.2015.06.008.
  * Hu and Zhang [2004] Feifang Hu and Li-Xin Zhang.  Asymptotic properties of doubly adaptive biased coin designs for multitreatment clinical trials.  _Annals of Statistics_ , 32(1):268–301, 2004.  doi:10.1214/aos/1079120137.
  * Klotz [1978] Jerome Klotz.  Maximum entropy constrained balance randomization in clinical trials.  _Biometrics_ , 34(2):283-287, 1978.  doi:10.2307/2530018.
  * Zhao et al. [2017] Wenle Zhao, Colin C. Everett, Yanqiu Weng, and Vance W. Berger.  Letter to the Editor – Guessing strategies for treatment prediction under restricted randomization with unequal allocation.  _Contemporary Clinical Trials_ , 59:118–120, 2017.  doi:10.1016/j.cct.2017.05.009.
  * Sverdlov and Ryeznik [2019] Oleksandr Sverdlov and Yevgen Ryeznik.  Implementing unequal randomization in clinical trials with heterogeneous treatment costs.  _Statistics in Medicine_ , 38(16):2905–2927, 2019.  doi:10.1002/sim.8160.
  * Rosenberger et al. [2019] William F. Rosenberger, Diane Uschner, and Yanying Wang.  Randomization: The forgotten component of the randomized clinical trial.  _Statistics in Medicine_ , 38(1):1–12, 2019.  doi:10.1002/sim.7901.
  * Proschan and Dodd [2019] Michael A. Proschan and Lori E. Dodd.  Re-randomization tests in clinical trials.  _Statistics in Medicine_ , 38(12):2292–2302, 2019.  doi:10.1002/sim.8093.