[R] Monotonic regression
Raubertas, Richard
richard_raubertas at merck.com
Mon May 9 19:51:03 CEST 2005
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
>
>
>
More information about the R-help
mailing list