×

## R Code for Cox & Stuart Test for Trend Analysis Article originally posted on Data Science Central. Visit Data Science Central

Below is an R code for Cox & Stuart Test for Trend Analysis. Simply, copy and paste the code into R workspace and use it. Unlike cox.stuart.test in R package named “randtests”, this version of the test does not return a p-value greater than one. This phenomenon occurs when the test statistic, T is half of the number of untied pairs, N.
Here is a simple example that reveals the situtaion:

> x
 1 4 6 7 9 7 1 6

> cox.stuart.test(x)

Cox Stuart test

data: x
statistic = 2, n = 4, p-value = 1.375
alternative hypothesis: non randomness

> o.Cox.Stuart.test(x)

Cox & Stuart Test for Trend Analysis

data: x
Test Statistic = 2, Number of Untied Pairs = 4, p-value = 1
alternative hypothesis: any type of trend, either decreasing or increasing

R Code for Cox & Stuart Test:

o.Cox.Stuart.test <- function(x, alternative=c(“two.sided” ,”left.sided”, “right.sided”)){

dname <- deparse(substitute(x))

alternative <- match.arg(alternative)

stopifnot(is.numeric(x))

n0 <- length(x)

if (n0 < 2){
stop(“sample size must be greater than 1”)
}

n0 <- round(length(x)) %% 2

if (n0 == 1) {
remove <- (length(x)+1)/2
x <- x[ -remove ]
}

half <- length(x)/2
x1 <- x[1:half]
x2 <- x[(half+1):(length(x))]
n <- sum((x2-x1)!=0)
t <- sum(x1<x2)

if (alternative == “left.sided”) {
p.value <- pbinom(t, n, 0.5)
alternative <- “decreasing trend”
}

if (alternative == “right.sided”){
p.value <- 1-pbinom(t-1, n, 0.5)
alternative <- “increasing trend”
}

if (alternative == “two.sided”){
alternative <- “any type of trend, either decreasing or increasing”
if (1-pbinom(t-1, n, 0.5) == pbinom(t, n, 0.5)) {
pdist <- dbinom(0:n, n, 0.5)
p.value <- sum(pdist[pdist <= t+1])
}
else {
p.value <- 2*min( 1-pbinom(t-1, n, 0.5), pbinom(t, n, 0.5))
}
}

rval <- list(statistic=c(“Test Statistic”=t),
alternative=alternative,
p.value=p.value,
method=”Cox & Stuart Test for Trend Analysis”,
parameter=c(“Number of Untied Pairs”=n),
data.name=dname)

class(rval) <- “htest”
return(rval)
}