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
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! > > ScottYou 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 ------------------------------
>>>>> "Scott" == Scott Briggs <scott.br at gmail.com> >>>>> on Sun, 8 May 2005 21:08:46 +0200 writes:Scott> Hi, I'm trying to find an implementation of monotonic regression in R Scott> and I haven't been able to find anything that's really related to Scott> this. isoMDS in the MASS package uses monotonic regression, however, Scott> I was wondering if there is any standalone function for monotonic Scott> regression? yes, isoreg() which actually was built by ``exporting'' the relevant bits & pieaces from MASS' isoMDS. Scott> Basically what I'm trying to do is implement monotonic regression Scott> where I can see not just the results of each iteration but also be Scott> able to tweak the input in order to test for or "kick" the regression Scott> out of a local minimum so that I can make sure I have the global Scott> minimum. (Ted Harding has already answered this part.)
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 > > >