[R] Monotonic regression

[Raubertas, Richard]

The 'pava' function below looks like code that I wrote
(with all the comments removed).  I have posted it two or
three times over the years to the S/R lists.

As Martin Maechler has noted, the function 'isoreg' will also
do monotonic regression (much faster for large data sets).
However, it does not allow weights (at least as of R 2.0.1).

I don't understand the original poster's comment about
"local minimum".  Isotonic regression is a convex optimization
problem and 'pava' or 'isoreg' will always produce the 
unique solution.

Rich Raubertas
Merck & Co.

> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch 
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of 
> Ted.Harding at nessie.mcc.ac.uk
> Sent: Sunday, May 08, 2005 5:14 PM
> To: Scott Briggs
> Cc: r-help at stat.math.ethz.ch
> Subject: RE: [R] Monotonic regression
> 
> 
> On 08-May-05 Scott Briggs wrote:
> > Hi, I'm trying to find an implementation of monotonic 
> regression in R
> > and I haven't been able to find anything that's really related to
> > this.  isoMDS in the MASS package uses monotonic 
> regression, however,
> > I was wondering if there is any standalone function for monotonic
> > regression?
> > 
> > Basically what I'm trying to do is implement monotonic regression
> > where I can see not just the results of each iteration but also be
> > able to tweak the input in order to test for or "kick" the 
> regression
> > out of a local minimum so that I can make sure I have the global
> > minimum.
> > 
> > Any help would be much appreciated.  Thanks!
> > 
> > Scott
> 
> You may probably find PAVA ("Pool Adjacent Violators Algorithm")
> useful. Below is code for a simple version which I have been using
> for a few years. I forget where I found it!
> 
> An R site search comes up with code for a version with more
> complex functionality at
> 
>   http://finzi.psych.upenn.edu/R/Rhelp02a/archive/9807.html
> 
> (contributed to r-help by Jan de Leeuw on 01 Jul 2004). I have
> not tested this code.
> 
> 
> pava<-function(x,wt=rep(1,length(x)))
> {
>   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)
>     if(!(any(viol))) break
>     i<-min( (1:(n-1))[viol])
>     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
> }
> # Examples:
> # > x<-c(1,0,1,0,0,1,0,1,1,0,1,0)
> # > x
> # [1] 1 0 1 0 0 1 0 1 1 0 1 0
> # > pava(x)
> #  [1] 0.4 0.4 0.4 0.4 0.4 0.5 0.5 0.6 0.6 0.6 0.6 0.6
> # > 
> # > pava(c(0,0,2/4,1/5,2/4,1/2,4/5,5/8,7/11,10/11),
> #        c(5,4,4,5,4,2,5,8,11,11))
> # [1] 0.0000000 0.0000000 0.3333333 0.3333333 0.5000000
> # [6] 0.5000000 0.6666667 0.6666667 0.6666667 0.9090909
> # (example where data are {ri,ni} so x={ri/ni} and w={ni})
> 
> 
> Best wishes,
> Ted.
> 
> 
> --------------------------------------------------------------------
> E-Mail: (Ted Harding) <Ted.Harding at nessie.mcc.ac.uk>
> Fax-to-email: +44 (0)870 094 0861
> Date: 08-May-05                                       Time: 20:56:28
> ------------------------------ XFMail ------------------------------
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html
> 
> 
>