Horace, Last year I found another Jonckheere-Terpstra function at http://www.biostat.wustl.edu/archives/html/s-news/2000-10/msg00126.html Using that and the version in SAGx and wanting a few other options, I created the version below.
Hope it is useful, Chris -- Christopher Andrews, PhD SUNY Buffalo, Department of Biostatistics ##### ##### ##### j.test <- function(x, y, alternative = c("two.sided", "decreasing", "increasing"), asymp=TRUE, correct=FALSE, perm=0, na.action=c("omit", "fail"), permgraph=FALSE, permreps=FALSE, ...) { #Jonckheere-Terpstra test #x = response #y = group membership #alternative hypothesis #asymptotic = asymptotic formula for variance or don't bother #correct = continuity correction #perm = number of repetitions for a permutation test #na.action = what to do with missing data #permgraph = draw a histogram of the permutations? #permreps = return the permutations? #... for graph alternative <- match.arg(alternative) complete <- complete.cases(x,y) nn <- sum(complete) cat(nn, "complete cases.\n") if (!all(complete) && (na.action=="fail")) { cat("option na.action is 'fail' and", sum(!complete), "cases are not complete.\n") return(NULL) } response <- x[complete] groupvec <- y[complete, drop=TRUE] if(!is.ordered(groupvec)) { groupvec <- as.ordered(groupvec) cat("y was not an ordered factor. Redefined to be one.\n") } lev <- levels(groupvec) nlev <- length(lev) if(nlev <= 2) { if (nlev < 2) { cat("Not enough groups (k =", nlev, ") for Jonckheere.\n") return(NULL) } else { cat("Two groups. Could use wilcox.test.\n") } } computestat <- function(response, groupvec, n.groups) { pairwise <- function(x, y) { length(x)*(length(y) + (1+length(x))/2) - sum(rank(c(x,y))[seq(along=x)]) } H<-0 for (i in seq(n.groups-1)) for (j in seq(i+1, n.groups)) H <- H+pairwise(response[groupvec == lev[i]], response[groupvec == lev[j]]) return(H) } retval <- list(statistic=computestat(response, groupvec, nlev), alternative = paste(alternative, paste(levels(groupvec),collapse=switch(alternative, two.sided = ", ", decreasing= " > ", increasing = " < ")), sep=": "), method = "Jonckheere", data.name = paste(deparse(substitute(x)), "by", deparse(substitute(y)))) if (asymp) { ns <- tabulate(as.numeric(groupvec)) retval$EH <- (nn * nn - sum(ns * ns))/4 retval$VH <- if (!any(duplicated(response))) { (nn * nn * (2 * nn + 3) - sum(ns * ns * (2 * ns + 3)))/72 } else { ds <- as.vector(table(response)) ((nn*(nn-1)*(2*nn+5) - sum(ns*(ns-1)*(2*ns+5)) - sum(ds*(ds-1)*(2*ds+5)))/72 + sum(ns*(ns-1)*(ns-2))*sum(ds*(ds-1)*(ds-2))/(36*nn*(nn-1)*(nn-2)) + sum(ns*(ns-1))*sum(ds*(ds-1))/(8*nn*(nn - 1))) } pp <- if (!correct) { pnorm(retval$statistic, retval$EH, sqrt(retval$VH)) } else if (retval$statistic >= retval$EH + 1) { pnorm(retval$statistic - 1, retval$EH, sqrt(retval$VH)) } else if (retval$statistic <= retval$EH - 1) { pnorm(retval$statistic + 1, retval$EH, sqrt(retval$VH)) } else { .5 } retval$p.value <- switch(alternative, two.sided = 2*min(pp,1-pp), decreasing = pp, increasing = 1-pp) } if (perm>0) { reps <- numeric(perm) for (i in seq(perm)) { reps[i] <- computestat(response, sample(groupvec), nlev) } retval$EH.perm <- mean(reps) retval$VH.perm <- var(reps) ppp <- if (!correct) { pnorm(retval$statistic, retval$EH.perm, sqrt(retval$VH.perm)) } else if (retval$statistic >= retval$EH.perm + 1) { pnorm(retval$statistic - 1, retval$EH.perm, sqrt(retval$VH.perm)) } else if (retval$statistic <= retval$EH.perm - 1) { pnorm(retval$statistic + 1, retval$EH.perm, sqrt(retval$VH.perm)) } else { .5 } retval$p.value.perm <- switch(alternative, two.sided = 2*min(ppp,1-ppp), decreasing = ppp, increasing = 1-ppp) rr <- rank(c(retval$statistic, reps))[1]/(perm+2) retval$p.value.rank <- switch(alternative, two.sided = 2*min(rr,1-rr), decreasing = rr, increasing = 1-rr) if (permreps) retval$reps <- reps if (permgraph) { hist(reps, xlab="Jonckheere Statistic", prob=TRUE, main="Histogram of Permutations", sub=paste(perm, "permutations")) abline(v=retval$statistic, col="red") points((-3:3)*sqrt(retval$VH.perm)+retval$EH.perm, rep(0,7), col="blue", pch=as.character(c(3:1,0:3))) } } class(retval) <- "htest" names(retval$statistic) <- "J" return(retval) #returns list (of class htest) with the following components #Always: #statistic = value of J #alternative = same as input #method = "Jonckheere" #data.name = source of data based on call #Most of the time (i.e., when asymp=TRUE): #EH = expected value based on sample size #VH = variance (adjusted for ties if necessary) #p.value = asymptotic p value #Sometimes (i.e., when perm>0): #EH.perm = expected value based on permutations #VH.perm = variance based on permutations #p.value.perm = p value based on normal approximation to permutations #p.value.rank = p value based on rank within permutation #perms = vector of permutation. } j.test(rnorm(54), ordered(rep(letters[1:27],2)), na.action="fail") j.test(rnorm(54), ordered(rep(letters[1:27],2)), na.action="omit", asymp=FALSE, perm=1000) j.test(ttt <- rnorm(54), ordered(rep(letters[1:2],27)), alt="d") wilcox.test(ttt ~ordered(rep(letters[1:2],27)), correct=FALSE, exact=FALSE, alt="g") #Higgins p 178 Example totals <- c(10,12,17,30,9,9,11,35,7,8,12,43) tlevels <- factor(c("A", "B", "C"), levels=c("A", "B", "C"), labels=c("A", "B", "C"), ordered=TRUE) Treatment <- rep( rep(tlevels, each=4) , times = totals) rlevels <- factor(c("Severe", "Moderate", "Slight", "None"), levels=rev(c("Severe", "Moderate", "Slight", "None")), labels=rev(c("Severe", "Moderate", "Slight", "None")), ordered=TRUE) Response <- rep( rep(rlevels, times=3) , times=totals) result <- j.test(Response, Treatment, alt="t", perm=1000, permgraph=TRUE, correct=TRUE) result # note that the standard print does not report the permutation pvalues. # But they are there result$p.value.perm result$p.value.rank ______________________________________________ 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