Up-and-Down Designs

Up-and-Down Designs (UDDs) are a family of statistical experiment designs used in dose-finding experiments in science, engineering, and medical research. Dose-finding experiments have binary responses: each individual outcome can be described as one of two possible values, such as success vs. failure or toxic vs. non-toxic. Mathematically the binary responses are coded as 1 and 0. The goal of dose-finding experiments is to estimate the strength of treatment (i.e., the 'dose') that would trigger the "1" response a pre-specified proportion of the time. This dose can be envisioned as a percentile of the distribution of response thresholds. An example where dose-finding is used: an experiment to estimate the LD₅₀ of some toxic chemical with respect to mice.

Simulated experiments from 3 different Up-and-Down designs. '0' and '1' responses are marked by 'o' and 'x', respectively. Top to bottom: the original "Simple" UDD that targets the median, a Durham-Flournoy Biased-Coin UDD targeting approximately the 20.6% percentile, and a k-in-a-row / "Transformed" UDD targeting the same percentile.

Dose-finding designs are sequential and response-adaptive: the dose at a given point in the experiment depends upon previous outcomes, rather than be fixed a priori. Dose-finding designs are generally more efficient for this task than fixed designs, but their properties are harder to analyze, and some require specialized design software. UDDs use a discrete set of doses rather than vary the dose continuously. They are relatively simple to implement, and are also among the best understood dose-finding designs. Despite this simplicity, UDDs generate random walks with intricate properties.[1] The original UDD aimed to find the median threshold by increasing the dose one level after a "0" response, and decreasing it one level after a "1" response. Hence the name "Up-and-Down". Other UDDs break this symmetry in order to estimate percentiles other than the median, or are able to treat groups of subjects rather than one at a time.

UDDs were developed in the 1940s by several research groups independently.[2][3][4] The 1950s and 1960s saw rapid diversification with UDDs targeting percentiles other then the median, and expanding into numerous applied fields. The 1970s to early 1990s saw little UDD methods research, even as the design continued to be used extensively. A revival of UDD research since the 1990s has provided deeper understanding of UDDs and their properties,[5] and new and better estimation methods.[6][7]

UDDs are still used extensively in the two applications for which they were originally developed: psychophysics where they are used to estimate sensory thresholds and are often known as fixed forced-choice staircase procedures,[8] and explosive sensitivity testing, where the median-targeting UDD is often known as the Bruceton test. UDDs are also very popular in toxicity and anesthesiology research.[9] They are also considered a viable choice for Phase I clinical trials.[10]

Mathematical Description

Definition

Let ${\displaystyle n}$ be the sample size of a UDD experiment, and assume for now that subjects are treated one at a time. Then the doses these subjects receive, denoted as random variables ${\displaystyle X_{1},\ldots ,X_{n}}$, are chosen from a discrete, finite set of ${\displaystyle M}$ increasing dose levels ${\displaystyle {\mathcal {X}}=\left\{d_{1},\ldots ,d_{M}:\ d_{1}<\cdots <d_{M}\right\}.}$ Furthermore, if ${\displaystyle X_{i}=d_{m}}$, then ${\displaystyle X_{i+1}\in \{d_{m-1},d_{m},d_{m+1}\},}$ according to simple constant rules based on recent responses. In words, the next subject must be treated one level up, one level down, or at the same level as the current subject; hence the name "Up-and-Down". The responses themselves are denoted ${\displaystyle Y_{1},\ldots ,Y_{n}\in \left\{0,1\right\};}$ hereafter we call the "1" responses positive and "0" negative. The repeated application of the same rules (known as dose-transition rules) over a finite set of dose levels, turns ${\displaystyle X_{1},\ldots ,X_{n}}$ into a random walk over ${\displaystyle {\mathcal {X}}}$. Different dose-transition rules produce different UDD "flavors", such as the three shown in the figure above.

Despite the experiment using only a discrete set of dose levels, the dose-magnitude variable itself, ${\displaystyle x}$, is assumed to be continuous, and the probability of positive response is assumed to increase continuously with increasing ${\displaystyle x}$. The goal of dose-finding experiments is to estimate the dose ${\displaystyle x}$ (on a continuous scale) that would trigger positive responses at a pre-specified target rate ${\displaystyle \Gamma =P\left\{Y=1\mid X=x\right\},\ \ \Gamma \in (0,1)}$; often known as the "target dose". This problem can be also expressed as estimation of the quantile ${\displaystyle F^{-1}(\Gamma )}$ of a cumulative distribution function describing the dose-toxicity curve ${\displaystyle F(x)}$. The density function ${\displaystyle f(x)}$ associated with ${\displaystyle F(x)}$ is interpretable as the distribution of response thresholds of the population under study.

The Transition Probability Matrix

Given that a subject receives dose ${\displaystyle d_{m}}$, denote the probability that the next subject receives dose ${\displaystyle d_{m-1},d_{m}}$, or ${\displaystyle d_{m+1}}$, as ${\displaystyle p_{m,m-1},p_{mm}}$ or ${\displaystyle p_{m,m+1}}$, respectively. These transition probabilities obey the constraints ${\displaystyle p_{m,m-1}+p_{mm}+p_{m,m+1}=1}$ and the boundary conditions ${\displaystyle p_{1,0}=p_{M,M+1}=0}$.

Each specific set of UDD rules enables the symbolic calculation of these probabilities, usually as a function of ${\displaystyle F(x)}$. Assume for now that transition probabilities are fixed in time, depending only upon the current allocation and its outcome, i.e., upon ${\displaystyle \left(X_{i},Y_{i}\right)}$ and through them upon ${\displaystyle F(x)}$ (and possibly on a set of fixed parameters). The probabilities are then best represented via a tri-diagonal transition probability matrix (TPM) ${\displaystyle \mathbf {P} }$:

${\displaystyle {\bf {{P}=\left({\begin{array}{cccccc}p_{11}&p_{12}&0&\cdots &\cdots &0\\p_{21}&p_{22}&p_{23}&0&\ddots &\vdots \\0&\ddots &\ddots &\ddots &\ddots &\vdots \\\vdots &\ddots &\ddots &\ddots &\ddots &0\\\vdots &\ddots &0&p_{M-1,M-2}&p_{M-1,M-1}&p_{M-1,M}\\0&\cdots &\cdots &0&p_{M,M-1}&p_{MM}\\\end{array}}\right).}}}$

The Balance Point

Usually, UDD dose-transition rules bring the dose down (or at least bar it from escalating) after positive responses, and vice versa. Therefore, UDD random walks have a central tendency: dose assignments tend to meander back and forth around some dose ${\displaystyle x^{*}}$ that can be calculated from the transition rules, when those are expressed as a function of ${\displaystyle F(x)}$.[1] This dose has often been confused with the experiment's formal target ${\displaystyle F^{-1}(\Gamma )}$, and the two are often identical - but they do not have to be. The target is the dose that the experiment is tasked with estimating, while ${\displaystyle x^{*}}$, known as the "balance point", is approximately where the UDD's random walk revolves around.[11]

The Stationary Distribution of Dose Allocations

Since UDD random walks are regular Markov chains, they generate a stationary distribution of dose allocations, ${\displaystyle \pi }$, once the effect of the manually-chosen starting dose wears off. This means, long-term visit frequencies to the various doses will approximate a steady state described by ${\displaystyle \pi }$. According to Markov chain theory the starting-dose effect wears off rather quickly, at a geometric rate.[12] Numerical studies suggest that it would typically take between ${\displaystyle 2/M}$ and ${\displaystyle 4/M}$ subjects for the effect to wear off nearly completely.[11] ${\displaystyle \pi }$ is also the asymptotic distribution of cumulative dose allocations.

UDD's central tendency ensures that long-term, the most frequently visited dose (i.e., the mode of ${\displaystyle \pi }$) will be one of the two doses closest to the balance point ${\displaystyle x^{*}}$.[1] If ${\displaystyle x^{*}}$ is outside the range of allowed doses, then the mode will be on the boundary dose closest to it. Under the original median-finding UDD, the mode will be at the closest dose to ${\displaystyle x^{*}}$ in any case. Away from the mode, asymptotic visit frequencies decrease sharply, at a faster-than-geometric rate. Even though a UDD experiment is still a random walk, long excursions away from the region of interest are very unlikely.

Examples of UDD stationary distributions with ${\displaystyle M=10}$. Left: original ("classical") UDD, ${\displaystyle x^{*}=5.6}$. Right: Biased-Coin targeting the 30th percentile, ${\displaystyle x^{*}\approx 3.9.}$

Common Up-and-Down Designs

The Original ("Simple" or "Classical") UDD

The original "simple" or "classical" UDD moves the dose up one level upon a negative response, and vice versa. Therefore, the transition probabilities are

$${\displaystyle {\begin{array}{rl}p_{m,m+1}&=P\{Y_{i}=0|X_{i}=d_{m}\}=1-F(d_{m});\\p_{m,m-1}&=P\{Y_{i}=1|X_{i}=d_{m}\}=F(d_{m}).\end{array}}}$$

We use the original UDD as an example for calculating the balance point ${\displaystyle x^{*}}$. The design's 'up', 'down' functions are ${\displaystyle p(x)=1-F(x),q(x)=F(x).}$ We equate them to find ${\displaystyle F^{*}}$:

$${\displaystyle 1-F^{*}=F^{*}\ \longrightarrow \ F^{*}=0.5.}$$

As stated earlier, the "classical" UDD is designed to find the median threshold. This is a case where ${\displaystyle F^{*}=\Gamma .}$

The "classical" UDD can be seen as a special case of each of the more versatile designs described below.

Durham and Flournoy's Biased Coin Design

This UDD shifts the balance point, by adding the option of treating the next subject at the same dose rather than move only up or down. Whether to stay is determined by a random toss of a metaphoric "coin" with probability ${\displaystyle b=P\{{\textrm {heads}}\}.}$ This biased-coin design (BCD) has two "flavors", one for ${\displaystyle F^{*}>0.5}$ and one for ${\displaystyle F^{*}<0.5,}$ whose rules are shown below:

$${\displaystyle X_{i+1}={\begin{array}{ll}d_{m+1}&{\textrm {if}}\ \ Y_{i}=0\ \ \&\ \ {\textrm {'heads'}};\\d_{m-1}&{\textrm {if}}\ \ Y\_i=1;\\d_{m}&{\textrm {if}}\ \ Y_{i}=0\ \ \&\ \ {\textrm {'tails'}}.\\\end{array}}}$$

The `heads' probability ${\displaystyle b}$ can take any value in${\displaystyle [0,1]}$. The balance point is

$${\displaystyle {\begin{array}{rcl}b\left(1-F^{*}\right)&=&F^{*}\\F^{*}&=&{\frac {b}{1+b}}\in [0,0.5].\end{array}}}$$

The BCD balance point can made identical to a target rate ${\displaystyle F^{-1}(\Gamma )}$ by setting the `heads' probability to ${\displaystyle b=\Gamma /(1-\Gamma )}$. For example, for ${\displaystyle \Gamma =0.3}$ set ${\displaystyle b=3/7}$. Setting ${\displaystyle b=1}$ makes this design identical to the classical UDD, and inverting the rules by imposing the coin toss upon positive rather than negative outcomes, produces above-median balance points. Versions with two coins, one for each outcome, have also been published, but they do not seem to offer an advantage over the simpler single-coin BCD.

Group (Cohort) UDDs

Some dose-finding experiments, such as Phase I trials, require a waiting period of weeks before determining each individual outcome. It may preferable then, to be able treat several subjects at once or in rapid succession. With group UDDs, the transition rules apply rules to cohorts of fixed size ${\displaystyle s}$ rather than to individuals. ${\displaystyle X_{i}}$ becomes the dose given to cohort ${\displaystyle i}$, and ${\displaystyle Y_{i}}$ is the number of positive responses in the ${\displaystyle i}$-th cohort, rather than a binary outcome. Given that the ${\displaystyle i}$-th cohort is treated at ${\displaystyle X_{i}=d_{m}}$ on the interior of ${\displaystyle {\mathcal {X}}}$ the ${\displaystyle i+1}$-th cohort is assigned to

$${\displaystyle X_{i+1}={\begin{cases}d_{m+1}&{\textrm {if}}\ \ Y_{i}\leq l;\\d_{m-1}&{\textrm {if}}\ \ Y_{i}\geq u;\\d_{m}&{\textrm {if}}\ \ Y_{i}<u.\end{cases}}}$$

${\displaystyle Y_{i}}$ follow a Binomial distribution conditional on ${\displaystyle X_{i}}$, with parameters${\displaystyle s}$ and${\displaystyle F(X_{i})}$. The `up' and `down' probabilities are the Binomial distribution's tails, and the `stay' probability its center (it is zero if ${\displaystyle u=l+1}$). A specific choice of parameters can be abbreviated as GUD${\displaystyle _{(s,l,u)}.}$

Nominally, group UDDs generate ${\displaystyle s}$-order random walks, since the ${\displaystyle s}$ most recent observations are needed to determine the next allocation. However, with cohorts viewed as single mathematical entities, these designs generate a first-order random walk having a tri-diagonal TPM as above. Some group UDD subfamilies are of interest:

• Symmetric designs with ${\displaystyle l+u=s}$ (e.g., GUD${\displaystyle _{(2,0,2)}}$) obviously target the median.
• The family GUD${\displaystyle _{(s,0,1)},}$ encountered in toxicity studies, allows escalation only with zero positive responses, and de-escalate upon any positive response. The escalation probability at ${\displaystyle x}$ is ${\displaystyle \left(1-F(x)\right)^{s},}$ and since this design does not allow for remaining at the same dose, at the balance point it will be exactly ${\displaystyle 1/2}$. Therefore,

$${\displaystyle F^{*}=1-\left({\frac {1}{2}}\right)^{1/k}.}$$

With ${\displaystyle s=2,3,4}$ would be associated with ${\displaystyle F^{*}\approx 0.293,0.206}$ and ${\displaystyle 0.159}$, respectively. The mirror-image family GUD${\displaystyle _{(s,s-1,s)}}$ has its balance points at one minus these probabilities.

For general group UDDs, the balance point can be calculated only numerically, by finding the dose ${\displaystyle x^{*}}$ with toxicity rate ${\displaystyle F^{*}}$ such that

$${\displaystyle \sum _{r=u}^{s}\left({\begin{array}{c}s\\r\\\end{array}}\right)\left(F^{*}\right)^{r}(1-F^{*})^{s-r}=\sum _{t=0}^{l}\left({\begin{array}{c}s\\t\\\end{array}}\right)\left(F^{*}\right)^{t}(1-F^{*})^{s-t}.}$$

Any numerical root-finding algorithm, e.g., Newton-Raphson, can be used to solve for ${\displaystyle F^{*}}$.[13]

The ${\displaystyle k}$-in-a-Row (or "Transformed" or "Geometric") UDD

This is the most commonly used non-median UDD. It was introduced by Wetherill in 1963,[14] and proliferated by him and colleagues shortly thereafter to psychophysics,[15] where it remains one of the standard methods to find sensory thresholds.[8] Wetherill called it "Transformed" UDD; Gezmu who was the first to analyze its random-walk properties, called it "Geometric" UDD in the 1990s;[16] and in the 2000s the more straightforward name "${\displaystyle k}$-in-a-row" UDD was adopted.[11] The design's rules are deceptively simple:

$${\displaystyle X_{i+1}={\begin{cases}d_{m+1}&{\textrm {if}}\ \ Y_{i-k+1}=\cdots =Y_{i}=0,\ \ {\textrm {all}}\ {\textrm {observed}}\ {\textrm {at}}\ \ d_{m};\\d_{m-1}&{\textrm {if}}\ \ Y_{i}=1;\\d_{m}&{\textrm {otherwise}},\end{cases}}}$$

In words, every dose escalation requires ${\displaystyle k}$ non-toxicities observed on consecutive data points, all at the current dose, while de-escalation only requires a single toxicity. It closely resembles GUD${\displaystyle _{(s,0,1)}}$ described above, and indeed shares the same balance point. The difference is that ${\displaystyle k}$-in-a-row can bail out of a dose level upon the first toxicity, whereas its group UDD sibling might treat the entire cohort at once, and therefore might see more than one toxicity before descending.

The method used in sensory studies is actually the mirror-image of the one defined above, with ${\displaystyle k}$ successive responses required for a de-escalation and only one non-response for escalation, yielding ${\displaystyle F^{*}\approx 0.707,0.794,0.841,\ldots }$ for ${\displaystyle k=2,3,4,\ldots }$.[17]

${\displaystyle k}$-in-a-row generates a ${\displaystyle k}$-th order random walk because knowledge of the last ${\displaystyle k}$ responses might be needed. It can be represented as a first-order chain with ${\displaystyle Mk}$ states, or as a Markov chain with ${\displaystyle M}$ levels, each having ${\displaystyle k}$ internal states labeled ${\displaystyle 0}$ to ${\displaystyle k-1}$ The internal state serves as a counter of the number of immediately recent consecutive non-toxicities observed at the current dose. This description is closer to the physical dose-allocation process, because subjects at different internal states of the level ${\displaystyle m}$, are all assigned the same dose ${\displaystyle d_{m}}$. Either way, the TPM is ${\displaystyle Mk\times Mk}$ (or more precisely, ${\displaystyle \left[(M-1)k+1)\right]\times \left[(M-1)k+1)\right]}$, because the internal counter is meaningless at the highest dose) - and it is not tridiagonal.

Here is the expanded ${\displaystyle k}$-in-a-row TPM with ${\displaystyle k=2}$ and ${\displaystyle M=5}$, using the abbreviation ${\displaystyle F_{m}\equiv F\left(d_{m}\right).}$ Each level's internal states are adjacent to each other.

${\displaystyle {\begin{bmatrix}F_{1}&1-F_{1}&0&0&0&0&0&0&0\\F_{1}&0&1-F_{1}&0&0&0&0&0&0\\F_{2}&0&0&1-F_{2}&0&0&0&0&0\\F_{2}&0&0&0&1-F_{2}&0&0&0&0\\0&0&F_{3}&0&0&1-F_{3}&0&0&0\\0&0&F_{3}&0&0&0&1-F_{3}&0&0\\0&0&0&0&F_{4}&0&0&1-F_{4}&0\\0&0&0&0&F_{4}&0&0&0&1-F_{4}\\0&0&0&0&0&0&F_{5}&0&1-F_{5}\\\end{bmatrix}}.}$

${\displaystyle k}$-in-a-row is often considered for clinical trials targeting a low-toxicity dose. In this case, the balance point and the target are not identical; rather, ${\displaystyle k}$ is chosen to aim close to the target rate, e.g., ${\displaystyle k=2}$ for studies targeting the 30th percentile, and ${\displaystyle k=3}$ for studies targeting the 20th percentile.

Estimating the Target Dose

Example for reversal-averaging estimation of a psychophysics experiment. Reversal points are circled, and the first reversal was excluded from the average. The design is a two-stage, with the second (and main) stage ${\displaystyle k}$-in-a-row targeting the 70.7% percentile. The first stage (until the first reversal) uses the "classical" UDD, a commonly-employed scheme to speed up the arrival to the region of interest.

Unlike other design approaches, UDDs do not have a specific estimation method "bundled in" with the design as a default choice. Historically, the more common choice has been some weighted average of the doses administered, usually excluding the first few doses to mitigate the starting-point bias. This approach antedates deeper understanding of UDDs' Markov properties, but its success in numerical evaluations relies upon the eventual sampling from ${\displaystyle \pi }$, since the latter is centered roughly around ${\displaystyle x^{*}.}$[5]

The single most popular among these averaging estimators was introduced by Wetherill et al. in 1966, and only includes reversal points (points where the outcome switches from 0 to 1 or vice versa) in the average.[18] See example on the right. In recent years, the limitations of averaging estimators have come to light, in particular the many sources of bias that are very difficult to mitigate. Reversal estimators suffer from both multiple biases (although there is some inadvertent cancelling out of biases), and increased variance due to using a subsample of doses. However, the knowledge about averaging-estimator limitations has yet to disseminate outside the methodological literature and affect actual practice.[5]

By contrast, regression estimators attempt to approximate the curve ${\displaystyle y=F(x)}$ describing the dose-response relationship, in particular around the target percentile. The raw data for the regression are the doses ${\displaystyle d_{m}}$ on the horizontal axis, and the observed toxicity frequencies,

$${\displaystyle {\hat {F}}_{m}={\frac {\sum _{i=1}^{n}Y_{i}I\left[X_{i}=d_{m}\right]}{\sum _{i=1}^{n}I\left[X_{i}=d_{m}\right]}},\ m=1,\ldots ,M,}$$

on the vertical axis. The target estimate is the abscissa of the point where the fitted curve crosses ${\displaystyle y=\Gamma .}$

Probit regression has been used for many decades to estimate UDD targets, although far less commonly than the reversal-averaging estimator. In 2002, Stylianou and Flournoy introduced an interpolated version of isotonic regression (IR) to estimate UDD targets and other dose-response data.[6] More recently, a modification called "centered isotonic regression" (CIR) was developed by Oron and Flournoy, promising substantially better estimation performance than ordinary isotonic regression in most cases, and also offering the first viable interval estimator for isotonic regression in general.[7] Isotonic regression estimators appear to be the most compatible with UDDs, because both approaches are nonparametric and relatively robust.[5] The publicly available R package "cir" implements both CIR and IR for dose-finding and other applications.[19]

References

1. Durham, SD; Flournoy, N. "Up-and-down designs. I. Stationary treatment distributions.". In Flournoy, N; Rosenberger, WF (eds.). IMS Lecture Notes Monograph Series. 25: Adaptive Designs. pp. 139–157.
2. Dixon, WJ; Mood, AM (1948). "A method for obtaining and analyzing sensitivity data". Journal of the American Statistical Association. 43: 109–126. doi:10.1080/01621459.1948.10483254.
3. von Békésy, G (1947). "A new audiometer". Acta Oto-Laryngologica. 35: 411–422. doi:10.3109/00016484709123756.
4. Anderson, TW; McCarthy, PJ; Tukey, JW (1946). 'Staircase' method of sensitivity testing (Technical report). Naval Ordnance Report. 65-46.
5. Flournoy, N; Oron, AP. "Up-and-Down Designs for Dose-Finding". In Dean, A (ed.). Handbook of Design and Analysis of Experiments. CRC Press. pp. 858–894.
6. Stylianou, MP; Flournoy, N (2002). "Dose finding using the biased coin up-and-down design and isotonic regression". Biometrics. 58: 171–177. doi:10.1111/j.0006-341x.2002.00171.x.
7. Oron, AP; Flournoy, N (2017). "Centered Isotonic Regression: Point and Interval Estimation for Dose-Response Studies". Statistics in Biopharmaceutical Research. 9: 258–267. arXiv:1701.05964. doi:10.1080/19466315.2017.1286256.
8. Leek, MR (2001). "Adaptive procedures in psychophysical research". Perception and Psychophysics. 63: 1279–1292. doi:10.3758/bf03194543.
9. Pace, NL; Stylianou, MP (2007). "Advances in and Limitations of Up-and-down Methodology: A Precis of Clinical Use, Study Design, and Dose Estimation in Anesthesia Research". Anesthesiology. 107: 144–152. doi:10.1097/01.anes.0000267514.42592.2a.
10. Oron, AP; Hoff, PD (2013). "Small-Sample Behavior of Novel Phase I Cancer Trial Designs". Clinical Trials. 10: 63–80. arXiv:1202.4962. doi:10.1177/1740774512469311.
11. Oron, AP; Hoff, PD (2009). "The k-in-a-row up-and-down design, revisited". Statistics in Medicine. 28: 1805–1820. doi:10.1002/sim.3590.
12. Diaconis, P; Stroock, D (1991). "Geometric bounds for eigenvalues of Markov chain". The Annals of Applied Probability. 1: 36–61. doi:10.1214/aoap/1177005980.
13. Gezmu, M; Flournoy, N (2006). "Group up-and-down designs for dose-finding". Journal of Statistical Planning and Inference. 6: 1749–1764.
14. Wetherill, GB; Levitt, H (1963). "Sequential estimation of quantal response curves". Journal of the Royal Statistical Society, Series B. 25: 1–48. doi:10.1111/j.2517-6161.1963.tb00481.x.
15. Wetherill, GB (1965). "Sequential estimation of points on a Psychometric Function". British Journal of Mathematical and Statistical Psychology. 18: 1–10. doi:10.1111/j.2044-8317.1965.tb00689.x.
16. Gezmu, Misrak (1996). The Geometric Up-and-Down Design for Allocating Dosage Levels (PhD). American University.
17. Garcia-Perez, MA (1998). "Forced-choice staircases with fixed step sizes: asymptotic and small-sample properties". Vision Research. 38 (12): 1861–81. doi:10.1016/s0042-6989(97)00340-4. PMID 9797963.
18. Wetherill, GB; Chen, H; Vasudeva, RB (1966). "Sequential estimation of quantal response curves: a new method of estimation". Biometrika. 53: 439–454. doi:10.1093/biomet/53.3-4.439.
19. Oron, Assaf. "Package 'cir'". CRAN. R Foundation for Statistical Computing. Retrieved 26 December 2020.