Example code in R for performing isotonic regressions and Jonckheere-Terpstra test

( taken from here.)

Example data and reading into R

grp <- c(rep(0,6),rep(1,6),rep(2,6))
dat <- c(40,35,38,43,44,41,38,40,47,44,40,42,48,40,45,43,46,44)

0 40
0 35
0 38
0 43
0 44
0 41
1 38
1 40
1 47
1 44
1 40
1 42
2 48
2 40
2 45
2 43
2 46
2 44

Alternatively you can just use a single series of data

grp <- c(1:10)
dat <- c(5,9,1,2,5,6,7,8,3,8)

NUMBER IS DAT (RESPONSE) AND INFORMATION ID THE GRP (OR GROUP) IN THE BELOW SCATTERPLOT

dat <- number
grp <- as.ordered(information)

R statements
 stripchart(grp ~ dat, method = "stack")
 stripchart(information ~ number, method = "stack")

The isotonic regression can then be performed in R by running the R syntax below

dat <- number
grp <- as.ordered(information)
mu0 <- mean(dat)
sig0 <- sd(dat)
n <- length(dat)
ss.null <- sum((dat - mu0)^2)

xbar <- sapply(split(dat, grp), mean)
nbar <- sapply(split(dat, grp), length)
k <- length(xbar)

pava <- function(x, w) {
    if (any(w <= 0))
        stop("weights must be positive")
    if (length(x) != length(w))
        stop("arguments not same length")
    n <- length(x)
    design <- diag(n)
    repeat {
        out <- lm(x ~ design + 0, weights = w)
        mu <- coefficients(out)
        dmu <- diff(mu)
        if (all(dmu >= 0)) break
        j <- min(seq(along = dmu)[dmu < 0])
        design[ , j] <- design[ , j] + design[ , j + 1]
        design <- design[ , - (j + 1), drop = FALSE]
    }
    return(as.numeric(design %*% mu))
}
 
test.stat <- function(x, w) {
    mu <- pava(x, w)
    mu0 <- sum(w * x) / sum(w)
    ss.alt <- sum(w * (x - mu)^2)
    ss.null <- sum(w * (x - mu0)^2)
    return(ss.null - ss.alt)
}

print(tstat <- test.stat(xbar, nbar))
print(mu) 

nsim <- 1e3 - 1
tsim <- double(nsim)
for (i in 1:nsim) {
    xbarsim <- rnorm(k, mu0, sig0 / sqrt(nbar))
    tsim[i] <- test.stat(xbarsim, nbar)
}
phat <- mean(tsim >= tstat)
(nsim * phat + 1) / (nsim + 1)
nsim / (nsim + 1) * sqrt(phat * (1 - phat) / nsim)
cat("Calculation took", proc.time()[1], "seconds\n")

and the Jonckheere-Terpstra test can be performed using the R code below

dat <- number
grp <- as.ordered(information)

jkstat <- function(x, g) {
    foo <- (sign(outer(x, x, "-")) + 1) / 2
    sum(foo[g[row(foo)] > g[col(foo)]])
}
 
print(jstat <- jkstat(dat, grp))

The results agree for this data set with ISOREG which was brought int to do isotonic regression after the R code in the other files was written. There are breaks at N=5 and N=9 giving

print(mu)
[1] 4.25 4.25 4.25 4.25 5.00 6.00 6.00 6.00 6.00 8.00

The main lesson here is that the normal theory test (isotonic regression) does more or less the same as the nonparametric Jonckheere-Terpstra test. This is no surprise since the data are fairly normal looking.