FAQ/isoR_eg - CBU statistics Wiki
location: FAQ / isoR_eg

# 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.

None: FAQ/isoR_eg (last edited 2013-03-08 10:17:24 by localhost)