[R] Monotonic regression in R??
Prof Brian D Ripley
ripley at stats.ox.ac.uk
Sun Mar 4 09:05:48 CET 2001
On Sat, 3 Mar 2001, Dirk Eddelbuettel wrote:
>
> "Ed" == Ed <M.> writes:
> Ed> Is there R code somewhere for monotonic regression? I have some
> Ed> functions (computer response time curves) that are guaranteed by theory
> Ed> to be monotonic, and some experimental data where they are not. What
> Ed> I'd like to do is plot a smooth monotonic curve through the
> Ed> experimental points. There is no closed form expression for the curve,
> Ed> although there is a recursive difference equation formulation in some
> Ed> simple cases, so ordinary nonlinear regression is not the answer. I
> Ed> have Hardle's book (the title of which escapes me at the moment) which
>
> If you are refering to
>
> @Book{haerdle90b,
> author = {Wolfgang H{\"a}rdle},
> title = {Smoothing Techniques: With Implementation in S},
> year = 1990,
> publisher = {Springer Verlag},
> address = {New York},
> series = {Springer Series in Statistics},
> }
(My copy has a 1991 date inside.)
> then there should be S code for it on Statlib. That should at least be a
> start.
(Yes, but the master is http://www.stats.ox.ac.uk/pub/S/haerdle.sh.Z. No
isotonic regression there that I can see, though.)
> I think the topic of his code came up hear before with the consensus that
> some of the more recent smoothing packages are actually better.
for smoothing and density estimation, much so. There are now better
statistical methods and much more careful implementations. I recommend
packages KernSmooth and sm, in particular.
The topic of isotonic regression comes up here and on S-news from time to
time. There are a few algorithms, not just one. One is implemented
as part of isoMDS in package MASS, but no direct interface is available.
Here is a pure S version of the pool-adjacent-violators algorithm from
http://www.biostat.wustl.edu/hyperlists/s-news/200002/msg00284.html
pava <- function (x, wt=rep(1,length(x)))
# Compute the isotonic regression of numeric vector 'x', with
# weights 'wt', with respect to simple order. The pool-adjacent-
# violators algorithm is used. Returns a vector of the same length
# as 'x' containing the regression.
# 02 Sep 1994 / R.F. Raubertas
{
n <- length(x)
if (n <= 1) return (x)
if (any(is.na(x)) || any(is.na(wt))) {
stop ("Missing values in 'x' or 'wt' not allowed")
}
lvlsets <- (1:n)
repeat {
viol <- (as.vector(diff(x)) < 0) # Find adjacent violators
if (!(any(viol))) break
i <- min( (1:(n-1))[viol]) # Pool first pair of violators
lvl1 <- lvlsets[i]
lvl2 <- lvlsets[i+1]
ilvl <- (lvlsets == lvl1 | lvlsets == lvl2)
x[ilvl] <- sum(x[ilvl]*wt[ilvl]) / sum(wt[ilvl])
lvlsets[ilvl] <- lvl1
}
x
}
Another version is given at
http://www.biostat.wustl.edu/hyperlists/s-news/200002/msg00282.html
That uses the S-PLUS function nnls.fit, but optim can be used: for the
conversion see the MASS script ch09.R.
Both of these are very much slower that the C code isoMDS uses, and I think
too that the algorithms are slower than the one that uses.
The search facility on S-news at
http://www.biostat.wustl.edu/s-news/s-news-archive/search.html
is very useful: might be worth a link on www.r-project.org.
--
Brian D. Ripley, ripley at stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272860 (secr)
Oxford OX1 3TG, UK Fax: +44 1865 272595
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list