Hello all.
Bill Venables very kindly revised my somewhat sloppy (and American-centric -
see below) isotonic regression function. Here is the revised version for
those who might be interested.
> Tim,
>
> Are you the same person who wrote the "isoreg" function? I
> presume so even if your email has changed. I was looking over
> the function and found one bug and a few ways that might improve
> it from a user point-of-view that you might like to consider.
>
> The original function was
>
> isoreg <- function(x, y, w = rep(1, length(x)), untie = "optimal", ...) {
> # Purpose: isotonic regression of y on x. Returns fitted values.
> # Arguments: x is the regressor and y is the regressand.
> # Weights are optional. Methods of dealing with ties in x
> # are 'optimal', 'random', and 'none' (i.e., preserve ties).
> # Author: Tim R. Johnson (tjohnson@s.psych.uiuc.edu)
> # Revised 2/24/99
> #
> x.tmp <- x[w != 0]
> y.tmp <- y[w != 0]
> w.tmp <- w[w != 0]
> n <- length(x.tmp)
> x.tmp <- switch(untie,
> optimal = rank(rank(x.tmp) + sort(runif(n))[rank(rank(y.tmp) +
> runif(n))]),
> random = rank(rank(x.tmp) + runif(n)),
> none = rank(x.tmp))
> if(untie == "none")
> m <- outer(x.tmp, sort(unique(x.tmp)), function(z1, z2)
> z1 >= z2)
> else m <- outer(x.tmp, 1:n, function(z1, z2)
> z1 >= z2)
> b <- nnls.fit(cbind(-1, m), y.tmp, weights = w.tmp, ...)$
> coefficients
> x[w != 0] <- as.vector(cbind(-1, m) %*% b)
> x
> }
>
> The bug is in the last two lines where I suppose you mean "y"
> where you have "x", although it will not make any difference
> unless you have some values excluded from consideration through
> zero weights.
>
> The revised version I suggest is as follows:
>
> isoreg <- function(x, y, weights = rep(1, length(x)),
> untie = c("optimal", "random", "none"),
> subset = weights > 0, ...) {
> # Purpose: isotonic regression of y on x. Returns fitted values.
> # Arguments: x is the predictor and y is the response.
> # Weights are optional. Methods of dealing with ties in x
> # are 'optimal', 'random', and 'none' (i.e., preserve ties).
> # Author: Tim R. Johnson (tjohnson@s.psych.uiuc.edu)
> # Revised 2/24/99
> #
> # Revised 26/2/2000 by Bill Venables.
> #
> x.tmp <- x[subset]
> y.tmp <- y[subset]
> w.tmp <- weights[subset]
> n <- length(x.tmp)
> untie <- match.arg(untie)
> x.tmp <- switch(untie,
> optimal = rank(rank(x.tmp) +
> sort(runif(n))[rank(rank(y.tmp) + runif(n))]),
> random = rank(rank(x.tmp) + runif(n)),
> none = rank(x.tmp))
> X <- cbind(-1, outer(x.tmp,
> if(untie == "none") sort(unique(x.tmp)) else 1:n,
> ">="))
> b <- nnls.fit(X, y.tmp, weights = w.tmp, ...)$coefficients
> y[subset] <- as.vector(X %*% b)
> y
> }
>
> The main cosmetic and other changes are:
>
> 1. untie may now be abbreviated to "n", "r" or "o"
>
> 2. weights and subset are available for separate purposes and
> have conventional argument names (so, for example, a user
> could deal with subsets defined by the levels of a factor
> without bothering about weights).
>
> 3. Calculating the model matrix for the nnls.fit call is just a
> little more swish,
>
> 4. The bug is fixed.
>
> Do with it what you wish. (There is another small bug: your
> revision date is given in the m/d/y form, which is considered
> just plain wrong everywhere else but in the USA! :-)
Indeed! ;) With the revised revision date we have:
isoreg <- function(x, y, weights = rep(1, length(x)),
untie = c("optimal", "random", "none"),
subset = weights > 0, ...) {
# Purpose: isotonic regression of y on x. Returns fitted values.
# Arguments: x is the predictor and y is the response.
# Weights are optional. Methods of dealing with ties in x
# are 'optimal', 'random', and 'none' (i.e., preserve ties).
# Author: Tim R. Johnson (tjohnson@s.psych.uiuc.edu)
# Revised 24/2/1999
#
# Revised 26/2/2000 by Bill Venables.
#
x.tmp <- x[subset]
y.tmp <- y[subset]
w.tmp <- weights[subset]
n <- length(x.tmp)
untie <- match.arg(untie)
x.tmp <- switch(untie,
optimal = rank(rank(x.tmp) +
sort(runif(n))[rank(rank(y.tmp) + runif(n))]),
random = rank(rank(x.tmp) + runif(n)),
none = rank(x.tmp))
X <- cbind(-1, outer(x.tmp,
if(untie == "none") sort(unique(x.tmp)) else 1:n,
">="))
b <- nnls.fit(X, y.tmp, weights = w.tmp, ...)$coefficients
y[subset] <- as.vector(X %*% b)
y
}
Best regards,
Tim
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Timothy R. Johnson
Division of Quantitative Psychology
Department of Psychology
University of Illinois @ Urbana-Champaign
mailto:tjohnson@s.psych.uiuc.edu
http://www.psych.uiuc.edu/~tjohnson/tjohnson.html
-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news
|