biviso: Bivariate isotonic regression.

In Iso: Functions to Perform Isotonic Regression

Bivariate isotonic regression.

Description

Bivariate isotonic regression with respect to simple (increasing) linear ordering on both variables.

Usage

biviso(y, w = NULL, eps = NULL, eps2 = 1e-9, ncycle = 50000,
       fatal = TRUE, warn = TRUE)

Arguments

y	The matrix of observations to be isotonized. It must of course have at least two rows and at least two columns.
w	A matrix of weights, greater than or equal to zero, of the same dimension as y. If left NULL then w is created as a matrix all of whose entries are equal to 1.
eps	Convergence criterion. The algorithm is deemed to have converged if each entry of the output matrix, after the completion of the current iteration, does not differ by more than eps from the corresponding entry of the matrix after the completion of the previous iteration. If this argument is not supplied it defaults to sqrt(.Machine$double.eps).
eps2	Criterion used to determine whether isotonicity is “violated”, whence whether (further) application of the “pool adjacent violators” procedure is required.
ncycle	The maximum number of cycles of the iteration procedure. Must be at least 2 (otherwise an error is given). If the procedure has not converged after ncycle iterations then an error is given. (See below.)
fatal	Logical scalar. Should the function stop if the subroutine returns an error code other than 0 or 4? If fatal is FALSE then output is returned by the function even if there was a “serious” fault. One can set fatal=FALSE to inspect the values of the objective matrix at various interim stages prior to convergence. See Examples.
warn	Logical scalar. Should a warning be produced if the subroutine returns a value of ifault equal to 4 (or to any other non-zero value when fatal has been set to FALSE)?

Details

See the paper by Bril et al., (References) and the references cited therein for details.

Value

A matrix of the same dimensions as y containing the corresponding isotonic values. It has an attribute icycle equal to the number of cycles required to achieve convergence of the algorithm.

Error Messages

The subroutine comprising Algorithm AS 206 produces an error code ifault with values from 0 to 6 The meaning of these codes is as follows:

• 0: No error.

• 1: Convergence was not attained in ncycle cycles.

• 2: At least one entry of w was negative.

• 3: Either nrow(y) or ncol(y) was less than 2.

• 4: A near-zero weight less than delta=0.00001 was replaced by delta.

• 5: Convergence was not attained and a non-zero weight was replaced by delta.

• 6: All entries of w were less than delta.

If ifault==4 a warning is given. All of the other non-zero values of ifault result in an error being given.

WARNING

This function appears not to achieve exact isotonicity, at least not quite. For instance one can do:

    set.seed(42)
    u  <- matrix(runif(400),20,20)
    iu <- biviso(u)
    any(apply(iu,2,is.unsorted))

and get TRUE. It turns out that columns 13, 14, and 16 of iu have exceptions to isotonicity. E.g. six of the values of diff(iu[,13]) are less than zero. However only one of these is less than sqrt(.Machine$double.eps), and then only “marginally” smaller.

So some of these negative values are “numerically different” from zero, but not by much. The largest in magnitude in this example, from column 16, is -2.217624e-08 — which is probably not of “practical importance”.

Note also that this example occurs in a very artificial context in which there is no actual isotonic structure underlying the data.

Author(s)

Rolf Turner [e-mail]

References

Bril, Gordon; Dykstra, Richard; Pillers Carolyn, and Robertson, Tim ; Isotonic regression in two independent variables; Algorithm AS 206; JRSSC (Applied Statistics), vol. 33, no. 3, pp. 352-357, 1984.

See Also

pava() pava.sa() ufit()

Examples

x <- 1:20
y <- 1:10
xy <- outer(x,y,function(a,b){a+b+0.5*a*b}) + rnorm(200)
ixy <- biviso(xy)

set.seed(42)
u <- matrix(runif(400),20,20)
v <- biviso(u)
progress <- list()
for(n in 1:9) progress[[n]] <- biviso(u,ncycle=50*n,fatal=FALSE,warn=FALSE)

Iso documentation built on Oct. 2, 2023, 9:06 a.m.