From: Raubertas, Richard <richard_raubertas_at_merck.com>

Date: Tue 10 May 2005 - 03:51:03 EST

R-help@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 Received on Tue May 10 04:00:26 2005

Date: Tue 10 May 2005 - 03:51:03 EST

Rich Raubertas

Merck & Co.

> -----Original Message-----

*> From: r-help-bounces@stat.math.ethz.ch
**> [mailto:r-help-bounces@stat.math.ethz.ch] On Behalf Of
**> Ted.Harding@nessie.mcc.ac.uk
**> Sent: Sunday, May 08, 2005 5:14 PM
**> To: Scott Briggs
**> Cc: r-help@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@nessie.mcc.ac.uk>
**> Fax-to-email: +44 (0)870 094 0861
**> Date: 08-May-05 Time: 20:56:28
**> ------------------------------ XFMail ------------------------------
**>
**> ______________________________________________
**> R-help@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
**>
**>
*

>

R-help@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 Received on Tue May 10 04:00:26 2005

*
This archive was generated by hypermail 2.1.8
: Fri 03 Mar 2006 - 03:31:40 EST
*