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.