 \documentclass[nojss]{jss}

\usepackage{graphicx} % subfigures
\usepackage{subcaption} % Figures
\usepackage{booktabs} % nice tables
\usepackage{multirow}
%\usepackage{csquotes} % quotes


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% declarations for jss.cls %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% almost as usual
\author{Diane Uschner\\RWTH Aachen University \And 
        David Schindler\\RWTH Aachen University \AND 
				Nicole Heussen \\RWTH Aachen University \And
				Ralf-Dieter Hilgers\\RWTH Aachen University
				}

\title{\pkg{randomizeR}: An R Package for the Assessment and Implementation of Randomization in Clinical Trials}
%\title{randomizeR: An R Package for the Assessment of Randomization Procedures and the Generation of Randomization Sequences}

%\VignetteEngine{knitr::knitr}
%\VignetteIndexEntry{Assessment and Implementation of Randomization in Clinical Trials}

%% for pretty printing and a nice hypersummary also set:
\Plainauthor{Diane Uschner, David Schindler, Nicole Heussen, Ralf-Dieter Hilgers} %% comma-separated
\Plaintitle{randomizeR: An R Package for the Assessment and Implementation of Randomization in Clinical Trials} %% without formatting
\Shorttitle{\pkg{randomizeR}: Randomization for Clinical Trials} %% a short title (if necessary)

%% an abstract and keywords
\Abstract{
This introduction to the \proglang{R} package \pkg{randomizeR} is a slightly
modified version of \cite{Uschner2016}.

Randomization in clinical trials is the key design technique to ensure the comparability
of treatment groups. 
%
Although there exist a large number of software products which assist the researcher
to implement randomization, no tool which covers a wide range of 
procedures and allows the comparative evaluation of the procedures under practical
restrictions has been proposed in the literature. 

The R package \pkg{randomizeR} addresses this need.
%
%
The paper includes a detailed description of the \pkg{randomizeR} package that serves
as a tutorial for the generation of randomization sequences and the assessment of randomization procedures. 
%Two parts: 1) Assessment of a randomization procedure 2) Implementation of the randomization procedure, i.e. Generation of randomization sequences
}

\Keywords{restricted randomization procedure, selection bias, chronological bias, \proglang{R}}
\Plainkeywords{restricted randomization procedure, selection bias, chronological bias, R} %% without formatting
%% at least one keyword must be supplied

%% publication information
%% NOTE: Typically, this can be left commented and will be filled out by the technical editor
%% \Volume{50}
%% \Issue{9}
%% \Month{June}
%% \Year{2012}
%% \Submitdate{2012-06-04}
%% \Acceptdate{2012-06-04}

%% The address of (at least) one author should be given
%% in the following format:
\Address{
  Diane Uschner\\
  Department of Medical Statistics\\
  RWTH Aachen University\\
	Pauwelsstra\ss e 30\\
  52068 Aachen, Germany\\
  E-mail: \email{duschner@ukaachen.de}\\
  URL: \url{http://www.ideal.rwth-aachen.de/}
}
%% It is also possible to add a telephone and fax number
%% before the e-mail in the following format:
%% Telephone: +43/512/507-7103
%% Fax: +43/512/507-2851

%% for those who use Sweave please include the following line (with % symbols):
%% need no \usepackage{Sweave.sty}

%% end of declarations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Mathematica Notation
\usepackage{amsthm, amsmath, amssymb}
\newtheorem{mydef}{Definition}
\newtheorem{Def}{Definition}

% Hadamard Matrix has more than 10 columns
\setcounter{MaxMatrixCols}{12}


\begin{document}

<<setup, include = FALSE, cache = FALSE>>=
library(knitr)
# set global chunk options
knitr::render_sweave() 
#opts_chunk$set(prompt=TRUE, comment=NA) #$ remove highlight and background to get standard knitr look.
options(prompt = "R> ", continue = "+  ", width = 70, useFancyQuotes = FALSE)
#options(replace.assign=TRUE, width=90, prompt="R> ")
set.seed(1986)
@

\section{Introduction}

Randomization is a design technique to ensure the comparability
of treatment groups in clinical trials by introducing a deliberate element of 
chance. \cite{Armitage1982} states the three main goals that are supposed to be achieved by randomization: First, 
it tends to balance known and unknown covariates and, thus, to produce 
structural equality of the treatment groups. Second, by ensuring effective 
blinding of treatment allocations from investigators and patients, randomization helps to 
avoid bias caused by the selection of patients. Finally, randomization contributes to the
internal validity of a trial that provides the basis for statistical inference.
The importance of randomization for clinical trials was first noted in the 
1940s by Sir A. Bradford Hill \citep[see][]{Chalmers1999} who realized that successful blinding
of treatment allocations was impossible without randomization.
%
Since that time, regulators have advocated the use of randomization in their 
guidelines \citep[see for example][]{ICHE9} and several different randomization 
procedures have been proposed in the literature.
It has been noticed that different randomization procedures behave differently, e.g.,
%It has been noticed that different randomization procedures have different properties, e.g.
concerning their susceptibility to bias and their potential to control power and type-I-error probability.
An overview containing the latest developments can be found in \cite{Rosenberger2016}. 

Despite the importance of randomization in the context of clinical trials, the 
novelties of the last decades have hardly found their way into clinical practice.
\cite{Berger2015} remark that clinicians stick with randomization procedures that 
are easily available, although their poor properties have been widely discussed
and better alternatives have been proposed. Notably, the
most suitable randomization procedure may depend on the clinical context. For
example, when block randomization is used, \cite{Tamm2014} advocate small block sizes to control chronological bias, 
while \cite{Kennes2012} show that larger blocks are better to control selection bias. 

The  choice of an adequate randomization procedure is crucial to achieve favorable properties 
such as control of power and type-I-error rate, especially for small sample sizes 
where asymptotic assumptions do not hold.
%
The R package \pkg{randomizeR} assists clinical trialists in making a scientifically 
sound choice based on these criteria. It combines the assessment of randomization procedures with the 
generation of randomization lists for clinical trials. 
%
Existing software tools such as the \pkg{blockrand} package by \cite{Snow2013}
implement only a very limited number of randomization procedures. The special class of 
response-adaptive randomization procedures is not included in \pkg{randomizeR},
but has currently been implemented for trials with time-to-event outcome in the 
MTLAB package \pkg{RARtool} by \cite{Ryeznik2015}.
To our knowledge, no software package for the assessment of the practical 
properties of randomization procedures exists, particularly not for the assessment 
of their exact properties in small clinical trials. 

This article gives a detailed description of \pkg{randomizeR}. 
%
It is organized as follows. Section \ref{sec:assessment} gives the background for assessing
the properties of randomization procedures, such as susceptibility to bias. The randomization procedures
are presented in Section \ref{ssec:rand-proc}. Section \ref{sec:package} 
shows how \pkg{randomizeR} implements the methods from the background 
chapter. Particularly, Section \ref{ssec:gen-one-seq} provides a tutorial
for the generation of randomization lists with \pkg{randomizeR}, and 
\ref{sec:assessment-randomizeR} illustrates the assessment of 
randomization procedures with respect to the different types of bias and other criteria.

\section{Assessing the properties of randomization procedures} \label{sec:assessment}

\subsection{Notation and model}
%
% Design Stage, N patients.
We consider the case of a two armed clinical trial with parallel group design.
%
Let $A$ and $B$ be treatments that influence a continuous outcome $Y$.
%
A \textit{randomization procedure} $\mathcal{M}$ is a probability distribution on the
space $\Omega = \{0,1\}^N$. A \textit{randomization sequence} is an element 
$t \in \Omega$, where $t_i=1$ if patient $i$ is allocated to treatment $A$ and
$t_i=0$ otherwise. Let $T= (T_1, \dots, T_N)$ denote a random variable with 
probability distribution $\mathcal{M}$ taking values in $\Omega$.
%
Throughout the paper, let the outcome $y_i$ of patient $i=1,\dots,N$ be
the realization of a normally distributed random variable
\begin{align}\label{eq:model}
Y_i \sim \mathcal{N}(\mu_A \cdot T_i + \mu_B \cdot (1 - T_i), \sigma^2)
\end{align}
%with \textit{treatment effect} $d \in \mathbb{R}$, overall mean $\mu >0$ and 
with group expectations $\mu_A, \mu_B$ and equal but unknown variance $\sigma^2>0$. 
%
The outcome $y_i$ is called \textit{response}.
Higher values of the response are regarded as better.
%
%
We test the null hypothesis that the expectation of the experimental
treatment $A$ does not differ from the expectation of the control treatment $B$
against the two-sided alternative 
%that the expectations of the treatment effects differ
\begin{gather}
H_0: \mu_A=\mu_B  \quad \text{vs.} \quad H_1: \mu_A\neq \mu_B.
\end{gather}

When the response $Y_i$ of patient $i$ is 
influenced by an unobserved quantity $b_i$, we call $b_i$ \textit{bias} of the $i$-th 
patient. Bias includes selection bias or chronological bias, or both, as proposed in 
the sections \ref{ssec:cb} and \ref{ssec:sb}. 
%
% Design stage, want to choose RP such that misspecification does not matter.
We investigate how the randomization 
procedures manage the misspecification of the model
%
\begin{align} \label{eq:biasmodel}
Y_i \sim \mathcal{N}(\mu_A \cdot T_i + \mu_B \cdot (1 - T_i) + b_i, \sigma^2)
\end{align}
with $b_i$  the fixed bias of patient $i$.
%
%We test the hypothesis of no treatment effect under model misspecification
%and evaluate the effect on the test-decision for different randomization procedures.
%

\subsection{Criteria for the assessment}\label{ssec:criteria}

The potential of a randomization procedure to control the impact of bias on the
study results along with other exemplary criteria for the choice of a 
randomization procedure are summarized in this section. The assessment of the impact of bias is important if the
form of the bias is unknown or the bias is unobserved.


%The following sections summarize three exemplary criteria the choice of a 
%randomization procedure should take into account.

\subsubsection{Susceptibility to chronological bias}
\label{ssec:cb}

%The sequential enrollment of patients in randomized controlled clinical trials 
%can cause the characteristics of the patients or the study environment to change  
%in the course of the trial. 
Changes in study environment, e.g., increased diagnostic potential, may impact
the response to treatment over time. 
%
Unobserved time trends lead to a bias of the estimator of the treatment effect, 
for which \cite{Matts1978} used the term  \textit{chronological bias}.
%
%The bias due to this change in covariates is called chronological bias. 
Chronological bias is a special case of accidental bias as introduced by 
\cite{Efron1971}. Efron investigated the effects of covariates that have been
(unintentionally) ignored in the model. Although the underlying trend could in 
theory be included in the model, the bias is often unobserved, or the exact form
is unknown. 
%
%Furthermore, especially in 
%small population groups, it is undesirable to introduce too many assumptions in
%the model. 
Furthermore, especially in  small population groups, it is challenging to use 
models with many explanatory variables.  
%\cite{Tamm2014} therefore propose to choose a randomization procedure
%that balances the trend in the groups.  
%
\cite{Rosenkranz2011} measured the impact of chronological bias by the distortion of the type-I-error
rate of the $t$~test when the responses are influenced by  a trend $\tau(i, \vartheta)$:
%under the null hypothesis $\mu=\mu_A=\mu_B$:
%
\begin{align}\label{eq:timetrend}
Y_i \sim \mathcal{N} (\mu_A \cdot T_i + \mu_B \cdot (1 - T_i) + \tau(i,\vartheta), \sigma^2 ).
\end{align}
\cite{Tamm2014} proposed three shapes of trend that are summarized in Table \ref{tab:shapes-of-trend}. 
Here $\log$ denotes the natural logarithm and $1_{\{i\geq n_0\}}$ is the indicator function
yielding the value one if $i\geq n_0$ and zero if $i < n_0$.
%
\begin{table}%
\centering
\begin{tabular}{ll}
\toprule
Trend & Shape\\
\midrule
\addlinespace
Linear & $\tau(i,\vartheta)= i\cdot\vartheta$ \\
Logarithmic & $\tau(i,\vartheta)= \log(\frac{i}{N})\cdot\vartheta$ \\
Step & $\tau(i,\vartheta)= 1_{\{i\geq n_0\}}\cdot \vartheta, n_0 \in \{1,\dots,N\}$\\
\addlinespace
\bottomrule
\end{tabular}
\caption{The three types of time trend included in \pkg{randomizeR}.}
\label{tab:shapes-of-trend}
\end{table}


\subsubsection{Susceptibility to selection bias}
\label{ssec:sb}

Sir A. Bradford Hill was the first to adopt randomization 
in clinical trials \citep[see][]{Chalmers1999}. His aim was to ensure effective blinding and to avoid bias 
due to the conscious or unconscious selection of patients to treatment groups, the so-called \textit{selection bias}. 
%
We consider two measures for 
selection bias proposed in the literature, the expected number of correct 
guesses and the influence of selection bias on the test decision.

The expected number of correct guesses was introduced by \cite{Blackwell1957}.
They assume that, based on past treatment assignments, the investigator 
consciously or unconsciously guesses the next treatment based on the past 
assignments. Suppose the investigator guesses that patient $i>1$ will receive 
treatment $g(t,i) \in \{A,B\}$ based on the previous assignments 
$(t_1, \dots, t_{i-1})$.
%
The correct guesses of a randomization sequence is the number of assignments
the investigator guesses correctly:
\begin{align}\label{eq:corGuess-seq}
CG(t) = |\{i \in \{1,\dots,N\}: g(t, i)= t_{i} \}|.
\end{align}

Two guessing strategies were investigated by \cite{Blackwell1957}. 
Under the \textit{convergence strategy (CS)}, the investigator assumes that the 
next patient is assigned to the group that has so far been assigned less.
%
Under the \textit{divergence strategy (DS)} the experimenter assumes that the next
patient is assigned to the treatment that has so far been observed more often.
%
At the beginning of the trial and when there is a tie, the investigator guesses the next allocation at random.

The expected number of correct guesses $E(CG(t))$ reflects the predictability of a 
sequence $t\in\Omega$.
%
The overall predictability of a randomization procedure $\mathcal{M}$ is given by the average 
proportion of correct guesses:
\[
avpCG_{\mathcal{M}} := \frac{1}{N} \cdot E_\mathcal{M} ( E(CG(t))).
\]

\cite{Proschan1994} proposed to measure the influence of selection bias on the 
test decision when the responses are biased as a result of selecting the patients
following the convergence strategy:
%
\begin{align} \label{eq:selbiasmodel}
Y_i \sim \mathcal{N}(\mu_A \cdot T_i  + \mu_B \cdot (1 - T_i) - \text{sign}(D_i)\cdot \eta, \sigma^2),
\end{align} 
%
where $D_i:= D_i(T) = \sum_{j=1}^i T_j - \sum_{j=1}^i (1-T_j)$ denotes the \textit{imbalance} of a 
randomization sequence, $\text{sign}(x)$ denotes the sign function, and $\eta$
denotes the selection effect.
% 
\cite{Tamm2012} demonstrated the impact of selection bias for different values of the selection effect $\eta$.


\subsubsection{Balancing behavior}
According to \cite{ICHE9}, it is desirable for a randomization procedure to 
balance the group sizes throughout the trial as well as at the end of the trial, 
while avoiding predictability. Table \ref{tab:types-of-imbalance} summarizes
the measures for the imbalance that have been proposed in the literature
\citep[see for example][]{Atkinson2014}.
%
According to \cite{Lachin1988} imbalance may cause decreased power of the 
statistical test in case of continuous endpoint and homoscedasticity.
\begin{table}%
\centering
\begin{tabular}{lcr}
\toprule
Imbalance & Formula \\
\midrule
Final & $D_N$\\
Absolute final & $|D_N|$\\
Maximal & $\max_{i=1,\dots,N} |D_i|$\\
\addlinespace
Loss & $ \frac{D_N^2}{N}$\\
\addlinespace
\bottomrule
\end{tabular}
\caption{Imbalance measures implemented in \pkg{randomizeR}.}
\label{tab:types-of-imbalance}
\end{table}


\subsection{Randomization procedures} \label{ssec:rand-proc}
Randomization procedures can be described in terms of a restricted or unrestricted 
random walk \citep[see][]{Proschan1994}. The restrictions imposed
on the random walk lead to different randomization procedures.
% 
In this section, we give a short overview about the randomization procedures that 
are implemented in \pkg{randomizeR}. For a comprehensive overview we refer to \cite{Rosenberger2016}.
%over most of the procedures
%covered here. %Let $N \in \mathbb{N}$ be the total sample size. 

\begin{figure}%
\centering
<<random-walk,echo=F, out.width='.6\\textwidth', message=F>>=
library(randomizeR)
#R<- genSeq(crPar(10),seed=126)
R<- genSeq(crPar(10),seed=808898100)
par(cex=1.5)
plotSeq(R,plotAllSeq = T)
@
\caption{Random walk of the randomization sequences of CR.}%
\label{fig:random-walk}%
\end{figure}

Complete Randomization (CR) is equivalent to tossing a fair coin for the 
allocation of each patient. CR leads to $2^N$ equi-probable sequences where $N$
denotes the total sample size.
Figure \ref{fig:random-walk} shows a randomization sequence produced by a CR in 
heavy black, along with all possible sequences in light gray.

Using Random Allocation Rule (RAR), patients are allocated by drawing $N$ times 
without replacement from an urn consisting of $N/2$ balls for each treatment. 
RAR produces $\binom{N}{N/2}$ equi-probable sequences that all attain final balance. 

Permuted Block Randomization (PBR) with \textit{block constellation} $bc=(k_1, \dots, k_m)$ 
balances the allocations in the blocks of length $k_1, \dots, k_m$.
%is similar to RAR. 
For each block $j=1, \dots, m$, an urn is filled with $k_j/2$ 
balls for each of the two treatments, and $k_j$ balls are drawn without replacement from 
the urn. PBR leads to $\prod_{j=1}^m \binom{k_j}{k_j/2}$ equi-probable sequences that attain
balance after each block, particularly at the end.

Permuted Block Randomization with \textit{random block constellation} (RPBR) is similar
to PBR, but the block constellation $bc$ is sampled at random from the set of given block 
lengths $rb$. RPBR permits two variations \citep{Heussen2004,Rosenberger2016}: The entries 
$k_j$ of $bc$ can either be drawn uniformly with replacement from $rb$ until 
$\sum k_j \geq N$, or be conditioned to achieve final balance, namely $\sum k_j =N$.

The Truncated Binomial Design (TBD) consists of tossing a fair coin for the allocation
of patients until $N/2$ heads or tails have occurred. The remaining patients are
allocated deterministically to the opposite group. TBD results in sequences that attain final balance
but are not equiprobable. TBD admits the same number of sequences as RAR. As an
extension, TBD can be conducted in blocks similar to PBR, or similar to RPBR with
random block constellation (RTBD).

The Maximal Procedure (MP) was proposed by \cite{Berger2003}. MP gives equal 
probability to all sequences that attain final balance and do not exceed a 
pre-specified maximum tolerated imbalance $mti \in \mathbb{N}$. Always when the 
imbalance boundary is reached, i.e., $|D_i|=mti$, a deterministic 
assignment is made to the underrepresented group in order to reduce the 
imbalance. 

The Big Stick Design (BSD) introduced by \cite{Soares1983} consists of tossing a 
fair coin until the imbalance $|D_i|$ reaches a  maximum tolerated imbalance $mti \in \mathbb{N}$. 

In Efron's Biased Coin design (EBC) \citep[see][]{Efron1971}, a biased coin is tossed for the next 
allocation in order to reduce the imbalance with probability $0.5\leq p\leq1$. 
When the groups are balanced, a fair coin is tossed.

\cite{Chen2007} proposed a biased coin design with imbalance tolerance (CHEN) 
that combines BSD and EBC. When the groups are balanced, a 
fair coin is tossed for the allocation of patients. Otherwise, a biased coin 
with probability  $0.5\leq p\leq1$ is tossed until reaching an imbalance boundary $mti \in \mathbb{N}$. 
When the imbalance boundary is reached, a deterministic assignment is made in
order to reduce the imbalance.

The Accelerated Biased Coin Design (ABCD) proposed by \cite{Antognini2004} uses
the toss of a biased coin for the allocation of patients. The probability 
$p_{i} =p(D_{i-1},a)$ for allocating the $i$th patient to the experimental group 
depends on the imbalance $D_{i-1}$ and an \textit{acceleration parameter} $a>0$. 
The acceleration parameter exponentially weights the imbalance and thus determines the degree of randomness of the 
design.

The Bayesian Biased Coin Design (BBCD) proposed under the name  'dominant 
biased coin design' by \cite{Antognini2014} is similar to ABCD. Here, the 
probability $p_{i}$ of allocating the $i$th patient to the experimental group 
depends on the acceleration parameter $a>0$ and the ratio $N_A(i-1)/ N_B(i-1)$ 
where $N_A(i-1)$ and $N_B(i-1)$ are the numbers of patients in groups $A$ and $B$ 
respectively after allocating $i-1$ patients.

Wei's Urn Design (UD) consists of $N$ draws of an urn whose composition is 
updated after each draw. Before allocating the first patient, the urn contains 
an initial number ($ini\geq 0$) of balls of different colors for each treatment. 
For the allocation of a patient, a ball is drawn, the color is noted, and 
replaced along with an additional number ($add\geq 0$) of balls of the of
the opposite color. 

The Generalized Biased Coin Design (GBCD) was developed by \cite{Smith1984} to
extend various designs. A biased coin is tossed for the allocation of patients,
where probability of allocating the next patient to the experimental group
depends on $N_A(i-1)$ and $N_B(i-1)$ as well as a balancing parameter 
$\rho \geq 0$. 


The Hadamard Randomization (HADA) proposed by \cite{bailey2003} uses the rows 
of a special Hadamard matrix $H \in \{0,1\}^{11\times 12}$ for the allocation of patients. Rows from $H$
are sampled until the number of allocations reaches the planned sample size $N$.

Table \ref{tab:rand-proc} in the next chapter shows how the presented randomization 
procedures can be used in \pkg{randomizeR}.


\subsection{Implementation}\label{ssec:implementation}

The assessment of a randomization procedure is based on a set of allocation sequences. 
Some of the assessment criteria, for example type-I-error or power, depend on a 
response, while others, e.g., correct guesses, are independent of a response.
%Consequently, we distinguish between the generation of the randomization sequences
%and the calculation of the assessment criteria. 
%
There are two different options to generate a set of allocation sequences
depending on the sample size $N$, and two different methods to calculate the
response based assessment criteria:

%Both can be conducted exactly or by simulation. 

\begin{itemize}
	\item \textbf{Complete or simulated reference set:}
		%
		%
		In case of $N \leq 24$, it is possible to generate the set of all possible allocation sequences $\{0,1\}^N$,
		assess the eligibility of an allocation sequence and, independently, calculate 
		the associated probabilities for the allocation sequences by the randomization algorithm.
		This results in the \textit{complete set} of the sequences. The function 
		\code{getAllSeq} provides this functionality, see Section \ref{ssec:gen-mult-sec}.

		In case of $N>24$, the complete set cannot be calculated in reasonable time.
		Instead, a number $r$ of allocation sequences is generated undergoing the 
		randomization algorithm of the particular randomization procedure. The function
		\code{genSeq} provides this functionality, see Section \ref{ssec:gen-mult-sec}.
		%
		Formally, the relative frequencies of the allocation sequences can be used to estimate the 
		true probabilities of allocation sequences. Note therefore that, in a simulation, 
		the relative frequency of a randomization sequence is used. This results in the
		\textit{simulated set} of the sequences. The simulated set can also be applied in case of 
		smaller sample sizes, and procedures where the exact approach is not available
		in the literature, i.e., Block Randomization with random block constellation and Hadamard Randomization.

\item \textbf{Exact or simulated response based assessment criteria:}
		%For each randomization sequences, the criteria from Section \ref{ssec:criteria} 
		%can be assessed.
		The susceptibility of a randomization procedure to bias can be measured as the
		distortion of the type-I- or type-II-error probability. 
		%
		The \textit{exact method} computes the distribution of the exact rejection 
		probabilities. For each randomization sequence in the reference set, the
		rejection probability of  Student's $t$~test is calculated using the knowledge on 
		the $b_i$ from Model \ref{eq:biasmodel} \citep[see][Chapter 4]{Langer2014}. 
		%
		The \textit{simulation method} simulates a response vector for each allocation
		sequence in the reference set according to Model \ref{eq:biasmodel} and derives a test decision of Student's $t$~test.
		The type-I- or type-II-error \textit{rate} is computed as the proportion of falsely rejected test decisions. 
		%
		%%
		%The \textit{simulation method} computes the type-I-error \textit{rate} as 
		%the proportion of falsely positive test decisions. 
		%%
		%For each randomization sequence of the reference set, one response vector is
		%simulated, and the test decision of Student's $t$~test is derived.
		%
		This method was used for example by \cite{Proschan1994}. Section 
		\ref{ssec:ass-rand-proc} shows how the exact and simulated type-I- and type-II-error
		probabilities can be assessed with \pkg{randomizeR}.
\end{itemize}

Any combination of the above methods can be used. For example, a simulated 
reference set can be used to assess the distribution of the exact rejection 
probabilities. The combination of simulated error rate with simulated 
reference set is usually used in the literature. For small sample sizes, the 
combination of exact reference set with distribution of exact error 
probabilities yields most accurate results.

All sampling algorithms use R's standard random number generator, the 
Mersenne-Twister \citep{RCT2016}.


\section{The randomizeR package}\label{sec:package}

\subsection{Overview}

The \pkg{randomizeR} package covers two closely connected purposes: The generation of 
randomization sequences and  the assessment of randomization procedures 
according to the aforementioned criteria.
%
The previous chapter was dedicated to introducing the basic terms and literature
from the field. In the present chapter, we show how \pkg{randomizeR} addresses
these purposes.

The current version of \pkg{randomizeR} is based on R 3.3.0. It can be loaded in an R 
session via:
<<loading>>=
library("randomizeR")
@

<<load-ggplot, message=FALSE, echo=FALSE>>=
library(ggplot2)
theme_update(text=element_text(size=24), axis.title.y=element_text(margin=margin(0,13,0,0)),  axis.title.x=element_text(margin=margin(13,0,0,0)))
@

All the main components of \pkg{randomizeR} are implemented using the  S4 object 
oriented system. 

\subsection{Generating randomization sequences}

There are two main purposes for the generation of randomization sequences. 
%
The first purpose is the generation of a single sequence for the allocation of 
patients in a clinical trial. 
%
The second purpose is the generation of multiple sequences
in order to assess the properties of a randomization procedure. 
%
Both purposes can be regarded as functions that use the randomization procedure 
itself as basis for their behaviour. They are therefore implemented as methods 
that take an object representing the randomization procedure as input.

\subsubsection{Representing randomization procedures}

\pkg{randomizeR} implements randomization procedures as subclasses of the 
\code{randPar} class. For example, an object representing Complete 
Randomization with sample size $N=10$ is generated by 
<<>>=
N <- 10 
(params <- crPar(N))
@
The function \code{crPar} is a so-called constructor function, namely a function
that generates an object of the class \code{crPar} and prepares it for use. 
The object \code{param} then contains all the information about the randomization 
procedure. 
%
%The number of treatment arms \code{K} and the names of the 
%treatment groups names \code{groups} are set default, but can be passed to 
%\code{crPar} if other values are desired. 
%
Table \ref{tab:rand-proc} summarizes the constructor functions for the 
randomization procedures described in Section \ref{ssec:rand-proc}. In \pkg{randomizeR}
an overview over the implemented randomization procedures is shown by \code{?randPar}.

\begin{table}%
\begin{tabular}{lll}
\toprule
Randomization procedure & Constructor function & Parameters\\
\midrule
Complete Randomization (CR) & \code{crPar(N)} & \code{N} sample size\\
\addlinespace
Random Allocation Rule (RAR) & \code{rarPar(N)} & \code{N} sample size\\
\addlinespace
Permuted Block Randomization (PBR) & \code{pbrPar(bc)} & \code{bc} block constellation\\
\addlinespace
\multirow{2}{*} {Rand. Permuted Block Randomization (RPBR)} & \multirow{2}{*} {\code{rpbrPar(N, rb)}} & \code{N} sample size \\ 
& & \code{rb} random block lengths \\
\addlinespace
Truncated Binomial Design (TBD) & \code{tbdPar(bc)} & \code{bc} block constellation\\
\addlinespace
\multirow{2}{*} {Rand. Truncated Binomial Design (RTBD)} & \multirow{2}{*} {\code{rtbdPar(N, rb)}} & \code{N} sample size \\ 
& & \code{rb} random block lengths \\
\addlinespace
\multirow{2}{*} {Maximal Procedure (MP)} & \multirow{2}{*} {\code{mpPar(N, mti)}} & \code{N} sample size \\ 
& & \code{mti} max. tolerated imbalance\\
\addlinespace
\multirow{2}{*} {Big Stick Design (BSD)} & \multirow{2}{*} {\code{bsdPar(N, mti)}} & \code{N} sample size \\ 
& & \code{mti} max. tolerated imbalance\\
\addlinespace
\multirow{2}{*} {Efrons Biased Coin Design (EBC)} & \multirow{2}{*} {\code{ebcPar(N, p)}} & \code{N} sample size \\ 
& & \code{p} biased coin probability\\
\addlinespace
\multirow{3}{*} {Chen's Design (CHEN)} & \multirow{3}{*} {\code{chenPar(N, mti, p)}} & \code{N} sample size \\ 
& & \code{mti} max. tolerated imbalance\\
& & \code{p} biased coin probability\\
\addlinespace
\multirow{2}{*} {Generalized Biased Coin Design (GBCD)} & \multirow{2}{*} {\code{gbcdPar(N, rho)}} & \code{N} sample size \\ 
& & \code{rho} balance factor\\
\addlinespace
\multirow{2}{*} {Adjustable Biased Coin Design (ABCD)} & \multirow{2}{*} {\code{abcdPar(N, a)}} & \code{N} sample size \\ 
& & \code{a} balance factor\\
\multirow{2}{*} {Bayesian Biased Coin Design (BBCD)} & \multirow{2}{*} {\code{bbcdPar(N, a)}} & \code{N} sample size \\ 
& & \code{a} balance factor\\
\addlinespace
\multirow{3}{*} {Wei's Urn Design (UD)} & \multirow{3}{*} {\code{udPar(N, ini, add)}} & \code{N} sample size \\ 
& & \code{ini} initial urn composition\\
& & \code{add} adjustment in each step\\
\addlinespace
Hadamard Randomization (HADA) & \code{hadaPar(N)} & \code{N} sample size\\
\bottomrule
\end{tabular}
\caption{Randomization procedures included in \pkg{randomizeR}.}
\label{tab:rand-proc}
\end{table}


\subsubsection{Generation of a single sequence} \label{ssec:gen-one-seq}

We can use the function \code{genSeq} to generate a single randomization sequence for 
a particular clinical trial. It takes an object representing a randomization
procedure as input and generates a sequence at random using this procedure. 
For example, the following code generates a randomization sequence
using Complete Randomization with sample size $N=10$ as above:
<<echo=2:3>>=
set.seed(112)
params <- crPar(N)
(R <- genSeq(params))
@

To ensure the reproducibility of the results and to enhance the reporting of the
randomization procedure, the function \code{genSeq} saves all the information 
that was used for the generation of the randomization sequence in the object \code{R} along with the
randomization sequence itself.

In order to obtain the randomization sequence stored in the object \code{R}, we can use the 
function \code{getRandList}:
<<>>=
getRandList(R)
@
The randomization sequence and the other information stored in the object 
\code{R} can conveniently be saved to a \code{.csv} file by using the \code{saveRand} function:
<<results=F>>=
saveRand(R, file="myRandList.csv")
@

Figure \ref{fig:random-walk} in Section \ref{ssec:rand-proc} shows the random walk of 
the randomization sequence \code{R}. This figure can be 
generated using the function \code{plotSeq} in \pkg{randomizeR}:
%
<<plotSeq,fig.keep='none'>>=
plotSeq(R, plotAllSeq = T)
@
%
The function \code{genSeq} can generate randomization sequences for all 
randomization procedures from Table \ref{tab:rand-proc}. The function \code{genSeq} 
has a method for each randomization procedure. For all randomization procedures, 
its output is an object of a class extending the class \code{randSeq}. 


\subsubsection{Generation of a set of sequences} \label{ssec:gen-mult-sec}

Randomization procedures can be assessed based on the set of sequences they produce 
(see Section \ref{ssec:implementation}).
\pkg{randomizeR} provides two ways to generate multiple randomization sequences 
from a specific procedure. 

For small sample sizes $N\leq 24$, we can use the function \code{getAllSeq} to
generate the complete set of randomization sequences. It takes an 
object representing a randomization procedure as input and calculates the 
complete set of randomization sequences of that procedure. For example, in the above case 
of Complete Randomization with sample size $N=10$, we run the statement 
<<>>=
(allSeqs <- getAllSeq(params))
@
to get the complete set of $2^{10} = 1024$ sequences.
%
Note that the function \code{getAllSeq} does not support the random block designs
RPBR and RTBD or the Hadamard Randomization, because no algorithms have been 
established in the literature.

In those cases where the enumeration of the complete set of sequences is computationally
intensive or algorithmically not feasible, the function \code{genSeq} can be
used to generate a simulated reference set. For example for sample size $N=50$ 
and Complete Randomization, a simulated reference set of size $r=10,000$
is generated by
<<>>=
N <- 50
params <- crPar(N)
(randomSeqs <- genSeq(params, r = 10000))
@
The parameter \code{seed} can be passed to \code{genSeq} to ensure the reproducibility of the results.

As \code{genSeq} samples randomization sequences with 
replacement and with the true probability of occurrence of the randomization 
procedure represented by \code{params}, the result may contain 
duplicates. More probable sequences will occur more frequently.

\pkg{randomizeR} provides the function \code{getProb} for the computation of the 
true probability of occurrence randomization of sequences. \code{getProb} 
takes any object that inherits from \code{randSeq}:

<<echo=2:3>>=
myPaste <- function(R) {apply(R, 1, function(x) paste(x, collapse = ""))}
p <- getProb(allSeqs)
head(data.frame(Sequences = myPaste(getRandList(allSeqs)), 
	        Probability = round(p, 6)))
@
If applied to an object resulting from \code{getAllSeq}, the 
resulting probabilities will always sum up to one. In contrary, if applied to 
\code{genSeq}, the sum of the probabilities will not equal one, because 
typically, the set of  sampled sequences is only a subset of the complete set of 
sequences. Note therefore that, in a simulation, the sampled frequency of a
randomization sequence is used instead of its true probability of occurrence.
Due to the nature of the sampling algorithms, the sampled frequencies converge
to the true probabilities if the sample is large enough. 


\subsection{Assessing randomization procedures} \label{ssec:ass-rand-proc}

Most of the assessment criteria in Section \ref{ssec:criteria} depend on the 
assumtion of normally distributed responses according to \eqref{eq:model}. 
Assume that both treatments have equal expectation $\mu_A = \mu_B = 0$ and 
variances $\sigma_A = \sigma_B = 1$. Then we can use the function \code{normEndp} 
representing the normal endpoint to pass these assumtions to \pkg{randomizeR} by setting
<<model, results='hide'>>=
muA <- muB <- 0
sigmaA <- sigmaB <- 1
normalEndpoint <- normEndp(mu = c(muA, muB), sigma = c(sigmaA, sigmaB))
@
The class \code{endpoint} provides flexibility for the extension to other endpoints.
Currently, only normal endpoints are available.

\pkg{randomizeR} implements the criteria for the evaluation of randomization procedures
as subclasses of the class \code{issue}. For example, an object representing 
 the exact rejection probability in the presence of chronological bias due to 
linear time trend that has a strength of $\vartheta = 1$, is generated by
<<>>=
(cb <- chronBias(type = "linT", theta = 1, method = "exact"))
@
The parameter \code{method} indicates whether the exact distribution of the
type-I-error rate should be calculated or whether the test decision should be
simulated by generating responses that are influenced by the trend (\code{method = "sim"}).
%
The function \code{chronBias} is a constructor function for objects of the
class \code{chronBias}. The object \code{cb} contains all the information about the bias. 
Table \ref{tab:criteria} summarizes the constructor functions for all the criteria
of assessment presented in Section \ref{ssec:criteria}. In \pkg{randomizeR}
an overview over the implemented assessment criteria is shown by \code{?issues}.

%%If \code{method = "exact"}, the exact rejection probability of a randomization
%%sequence is calculated, yielding the true distribution of the power loss due
%%to imbalance of a randomization sequence. If \code{method = "sim"}, the 
%%test decision is simulated as described in section \ref{ssec:model}. % paragraph \ref{sssec:testdec}.
%
%
%The output is formally an S4 object representing the correct guesses. The 
%\code{type} attribute can take the values \code{"CS"} and \code{"DS"} 
%according to the guessing strategy.

\begin{table}%
\begin{tabular}{lll}
\toprule
Criterion & Constructor Function & Parameters\\
\midrule
\multirow{3}{*}{Chronological Bias} & \multirow{3}{*}{\code{chronBias(type,theta, method)}} & \code{type} of trend\\
  & & \code{theta} strength of trend\\
	& & \code{method} of assessment\\
\addlinespace
\multirow{3}{*}{Selection Bias} & \multirow{3}{*}{\code{selBias(type,eta, method)}} & \code{type} of guessing strategy\\
  & & \code{eta} selection effect \\
	& & \code{method} of assessment\\
\addlinespace
Correct Guesses & \code{corGuess(type)} & \code{type} of guessing strategy\\
\addlinespace
Imbalance & \code{imbal(type)} & \code{type} of imbalance\\
\addlinespace
\multirow{2}{*}{Combined Bias} & \multirow{2}{*}{\code{combineBias(selBias, chronBias)}} & \code{selBias} object\\
  & & \code{chronBias} object\\
\addlinespace
\multirow{2}{*}{Power loss due to imbalance} & \multirow{2}{*}{\code{setPower(d, method)}} & \code{d} detectable effect\\
	& & \code{method} of assessment\\
	\bottomrule
\end{tabular}
\caption{Criteria for the assessment.}
\label{tab:criteria}
\end{table}



\subsubsection{Assessment of a randomization procedure}\label{sec:assessment-randomizeR}
The \pkg{randomizeR} package includes the function \code{assess} to evaluate the 
behavior of a randomization procedure with respect to one or more of the criteria from Section \ref{ssec:criteria}.

For example, if we want to evaluate the behavior of the Big Stick Design with
sample size $N=12$ and imbalance tolerance $mti=2$ with respect to chronological 
bias, selection bias, and power loss given a difference in group means $d=1.796$, we may call:

<<parameters>>=
N <- 12
mti <- 2
bsdSeq <- getAllSeq(bsdPar(N, mti))

d <- 1.796 
sb <- selBias("CS", eta = d/4, method = "exact")
cb <- chronBias("linT", theta = 1/N, method = "exact") 
pw <- setPower(d, method = "exact") 
@

The function \code{assess} yields the true rejection probability for each 
criterion and each randomization sequence.
<<assessment>>=
(A <-assess(bsdSeq, sb, cb, pw, endp = normalEndpoint))
@
In the case of the Big Stick Design with sample size $N=12$ and maximum tolerated
imbalance $mti=2$, there are 972 possible sequences.
The first column of the assessment corresponds to the randomization sequence and the second to its
probability of occurrence of the sequence. The following columns correspond to
the criteria \code{sb}, \code{cb} and \code{pw} passed to \code{assess}.
The notation \code{P(rej)(type)} refers to the probability of rejection in the presence  of
the given \code{type} of bias. The column \code{power(exact)} gives the
exact power of each randomization sequence. Any number of assessment criteria 
can be passed to \code{assess}. For the criteria \code{imbal} and \code{corGuess} 
the endpoint \code{endp} is not relevant and can be omitted.

The \code{summary} of the assessment shows the important characteristics of the 
distribution of the allocation sequences such as mean, standard deviation, minimum, maximum and the quantiles 
for each criterion:
<<summary-assess>>=
summary(A)
@
For example, the five percent quantile \code{x05} of the 
probability of rejection in the presence of a linear time trend \code{linT}
is $0.042$.
%
Note that, until now, only the mean value of each criterion has been studied in 
the literature.
\subsubsection{Comparison of randomization procedures}

\pkg{randomizeR} provides the function \code{compare} for the comparison of several 
randomization procedures with respect to one of the criteria from Section \ref{ssec:criteria}. 
%
For example, assume we are in the same setting as above ($N=12$) and want to compare 
the Big Stick Design, the Maximal Procedure with $mti=2$ and the Permuted Block Randomization with block size four
with respect to their susceptibility to selection bias. We can 
partly recycle the previous code and set the parameters for the other randomization procedures
we want to compare:
<<parameters-comp>>=
mpSeq <- getAllSeq(mpPar(N, mti))
bc <- rep(4, N/4) 
pbrSeq <- getAllSeq(pbrPar(bc))
@

Then the following code compares the aforementioned procedures with respect to
selection bias:

<<comparison>>=
(C <- compare(sb, bsdSeq, mpSeq, pbrSeq, endp = normalEndpoint))
@

The distribution of a criterion for a particular randomization procedure can be
visualized and compared between different randomization procedures via boxplot (Figure \ref{fig:comp-box}) or
violin plot (Figure \ref{fig:comp-viol}):
<<plot-comparison-hide, fig.keep='none', warning=FALSE>>=
plot(C)
plot(C, y = "boxplot")
@

For the generation of boxplot and violin plot, we used the functions 
\code{geom_boxplot} and \code{geom_violin} of the R package \pkg{ggplot2} 
proposed by \cite{Wickham2009}.

\begin{figure}%

\centering
	\begin{subfigure}[b]{0.45\textwidth}
<<plot-comparison-show-violin, echo=FALSE>>=
plot(C)
@
	\caption{Violin plot}
	\label{fig:comp-viol}
\end{subfigure}
\hspace{1em}
\begin{subfigure}[b]{0.45\textwidth}
%<<plot-comparison-show-box, echo=FALSE, warning=FALSE, dev='pdf', dev.args=list(pointsize=22)>>=
<<plot-comparison-show-box, echo=FALSE, warning=FALSE>>=
plot(C, y='boxplot')
@
	\caption{Boxplot}
	\label{fig:comp-box}
\end{subfigure}
\caption{Comparative visualization of the distribution of the type-I-error probability under the influence of unobserved selection bias for BSD(2), MP(2) and PBR(4).}%
\label{}%
\end{figure}



\section{Summary and further research}

\pkg{randomizeR} is an R package that facilitates the generation
of randomization lists for a large number of randomization procedures and 
makes the assessment of randomization procedures with respect to various criteria
possible. The package currently implements 15 randomization procedures and six
criteria for the assessment of the procedures. It assists researchers at the 
design stage of a clinical trial by letting the choice of a randomization 
procedure and the implementation of the design go hand in hand. 

We are working towards extending \pkg{randomizeR} in various directions. The object 
oriented approach makes it easy to add new  randomization procedures and 
assessment criteria. A unified assessment criterion may be included to uniformly judge the 
suitability of a randomization procedure based on various criteria. The models 
may be extended to other endpoints such as time-to-event data, or to more than 
two treatment groups. Finally, randomization tests may be implemented to enable 
randomization based inference independent of parametric assumptions. 

\section*{Acknowledgement}
This research has received funding within the IDeAl FP7 project from the 
European Union's Seventh Framework Programme for research, technological 
development and demonstration under Grant Agreement no 602552. 

\bibliography{references}
\end{document}