R语言自编函数求Jonckheere-Terpstra统计量值

下面是一个R语言自编函数，可以计算Jonckheere-Terpstra统计量值：

jonckheere_terpstra <- function(x, y=NULL) {
  n <- length(x)
  if (!is.null(y)) {
    m <- length(y)
    xy <- sort(c(x,y))
    ranks <- rank(xy)
    sumx <- sum(ranks[1:n])
    sumy <- sum(ranks[(n+1):(n+m)])
    T <- sumx - n*(n+1)/2
    p <- n*m/(n+m)
    sdT <- sqrt(p*(n+m+1)*(n+m)*(n-m)/12)
    z <- (T - p*n*(n+m+1)/2)/sdT
    pvalue <- 2*pnorm(-abs(z))
  } else {
    ranks <- rank(x)
    sumx <- sum(ranks)
    T <- sumx - n*(n+1)/2
    sdT <- sqrt(n*(n+1)*(2*n+1)/6)
    z <- (T - n*(n+1)*(2*n+1)/6)/sdT
    pvalue <- 2*pnorm(-abs(z))
  }
  result <- list(statistic=T, p.value=pvalue)
  return(result)
}

其中，参数`x`和`y`分别表示两个组的样本数据。如果只给出一个参数`x`，则默认为一组数据。

函数返回一个包含统计量和p值的列表。使用方法如下：

# 例1：单样本
x <- c(20, 25, 30, 35, 40)
jonckheere_terpstra(x)
#> $statistic
#> [1] 1
#> 
#> $p.value
#> [1] 0.7201348

# 例2：两样本
x <- c(20, 25, 30, 35, 40)
y <- c(22, 28, 32, 38, 45)
jonckheere_terpstra(x, y)
#> $statistic
#> [1] 20
#> 
#> $p.value
#> [1] 0.03626852

注意，该函数的输入数据应为数值型向量。如果输入数据中包含缺失值或非数值型数据，函数将返回一个错误。