[R] Monotonic regression in R??

[Prof Brian D Ripley]

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._