Isotonic regression

Problem statement and algorithms

Let ${\displaystyle (x_1,y_1),\dots ,(x_n,y_n)}$ be a given set of observations, where the ${\displaystyle y_i\in \mathbb{R}}$ and the ${\displaystyle x_i}$ fall in some partially ordered set. For generality, each observation ${\displaystyle (x_i,y_i)}$ may be given a weight ${\displaystyle w_i\ge 0}$, although commonly ${\displaystyle w_i=1}$ for all ${\displaystyle i}$.

Isotonic regression in the ${\displaystyle L_2}$ metric seeks a weighted least-squares fit ${\displaystyle {\hat{y}}_i\approx y_i}$ for all ${\displaystyle i}$, subject to the constraint that ${\displaystyle {\hat{y}}_i\le {\hat{y}}_j}$ whenever ${\displaystyle x_i\le x_j}$. This gives the following quadratic program (QP) in the variables ${\displaystyle {\hat{y}}_1,\dots ,{\hat{y}}_n}$:

${\displaystyle min\sum _{i=1}^n{w}_i({\hat{y}}_i-y_i{)}^2}$ subject to ${\displaystyle {\hat{y}}_i\le {\hat{y}}_j\text{ for all }(i,j)\in E}$

where ${\displaystyle E=\{(i,j):x_i\le x_j\}}$ specifies the partial ordering of the observed inputs ${\displaystyle x_i}$ (and may be regarded as the set of edges of some directed acyclic graph (dag) with vertices ${\displaystyle 1,2,\dots n}$). Problems of this form may be solved by generic quadratic programming techniques.

In the usual setting where the ${\displaystyle x_i}$ values fall in a totally ordered set such as ${\displaystyle \mathbb{R}}$, we may assume WLOG that the observations have been sorted so that ${\displaystyle x_1\le x_2\le \cdots \le x_n}$, and take ${\displaystyle E=\{(i,i+1):1\le i<n\}}$. In this case, a simple and efficient iterative algorithm for solving the quadratic program (of a "least squares" optimisation in the ${\displaystyle L_2}$ metric) is the pool adjacent violators algorithm. Conversely, Best and Chakravarti studied the problem as an active set identification problem, and proposed a primal algorithm. These two algorithms can be seen as each other's dual, and both have a computational complexity of ${\displaystyle O(n)}$ on already sorted data.

To complete the isotonic regression task in the ${\displaystyle L_2}$ metric, we may then choose any non-decreasing function ${\displaystyle f(x)}$ such that ${\displaystyle f(x_i)={\hat{y}}_i}$ for all i. Any such function obviously solves

${\displaystyle \underset{f}{min}\sum _{i=1}^n{w}_i(f(x_i)-y_i{)}^2}$ subject to ${\displaystyle f}$ being nondecreasing

and can be used to predict the ${\displaystyle y}$ values for new values of ${\displaystyle x}$. A common choice when ${\displaystyle x_i\in \mathbb{R}}$ would be to interpolate linearly between the points ${\displaystyle (x_i,{\hat{y}}_i)}$, as illustrated in the figure, yielding a continuous piecewise linear function:

${\displaystyle f(x)=\{\begin{array}{ll}{\hat{y}}_1& \text{if }x\le x_1\\ {\hat{y}}_i+\frac{x-x_i}{x_{i+1}-x_i}({\hat{y}}_{i+1}-{\hat{y}}_i)& \text{if }{x}_i\le x\le x_{i+1}\\ {\hat{y}}_n& \text{if }x\ge x_n\end{array}}$

It is possible -- but more time-consuming -- to compute an isotonic regression in a metric other than ${\displaystyle L_2}$. An ${\displaystyle L_1}$ fit has a least absolute deviation and may thus be more robust to large outliers than an ${\displaystyle L_2}$ fit with its least squares criterion.

Another generalisation is to unimodal regressions, i.e. to fits with a non-decreasing segment followed by a non-increasing segment (or vice versa).

Centered isotonic regression

As this article's first figure shows, in the presence of monotonicity violations the resulting interpolated curve will have flat (constant) intervals. In dose-response applications it is usually known that ${\displaystyle f(x)}$ is not only monotone but also smooth. The flat intervals are incompatible with ${\displaystyle f(x)}$'s assumed shape, and can be shown to be biased. A simple improvement for such applications, named centered isotonic regression (CIR), was developed by Oron and Flournoy and shown to substantially reduce estimation error for both dose-response and dose-finding applications. Both CIR and the standard isotonic regression for the univariate, simply ordered case, are implemented in the R package "cir". This package also provides analytical confidence-interval estimates.