|
John,
Here
is code for the pool adjacent violators algorithm implemented
entirely in S.
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
}
I would be please to learn of any
Splus/R code which implements isotonic regression without recourse to compiled
code. I have downloaded (from Statlib) a Fortran version of the
Pool-Adjacent-Violators algorithm as described in the standard textbook on the
subject:-
Barlow RE, Bartholomew DJ, Bremner
JM and Brunk HD, Statistical Inference under Order Restriction. John Wiley
& Sons 1972.
I was wondering if anyone had
coded a version of this or a similar algorithm in S+/R? The Fortran code
worked well but I wonder if it was necessary.
John Steward
-----------------------------------------------
Dr J A Steward MBBCh BA MSc PhD FFPH
Director Welsh Cancer Intelligence & Surveillance
Unit
Rombourne House
14 Cathedral Rd
Cardiff CF11 9LJ
------------------------------------------------------
tel +44 (0)29 20373500
fax +44 (0)29 20373511
email:-
john.steward@velindre-tr.wales.nhs.uk
website:-
http://www.wcisu.wales.nhs.uk
------------------------------------------------------------------------------
Notice: This e-mail message, together with any attachments, contains information of Merck & Co., Inc. (One Merck Drive, Whitehouse Station, New Jersey, USA 08889), and/or its affiliates (which may be known outside the United States as Merck Frosst, Merck Sharp & Dohme or MSD and in Japan, as Banyu) that may be confidential, proprietary copyrighted and/or legally privileged. It is intended solely for the use of the individual or entity named on this message. If you are not the intended recipient, and have received this message in error, please notify us immediately by reply e-mail and then delete it from your system.
------------------------------------------------------------------------------ |