目的

Jonckheere 検定を行う。

使用法

Jonckheere(dat, correct=FALSE, alternative = c("increasing", "decreasing"))

引数

x           測定データ
g           群を表す数値。群を表す変数の取る値は「対照群＜処置群1＜処置群2 ...」のように決めておく。
            逆でも良いが，以下の説明を全部逆に解釈することになり間違いが生じやすいので注意のこと。
correct     連続性の補正をする場合には TRUE を指定する（デフォールトは連続性の補正をしない）
alternative 対立仮説の指示（デフォールトは increasing）
　　increasing   H1: 対照群≦処置群1≦処置群2 ...
　　decreasing   H1: 対照群≧処置群1≧置群2 ...

注：increasing と decreasing を間違えて選択しないように注意すること。
注：両側検定を行うときは，increasing と decreasing を指定して二回検定を行い，いずれかの有意確率が α/2 以下ならば有意と判定する。

ソース

インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/Jonckheere.R", encoding="euc-jp")

# ヨンキー検定
Jonckheere <- function(      x,                                              # データベクトル
                        g,                                              # 群変数ベクトル
                        correct=FALSE,                                  # 連続性の補正
                        alternative = c("increasing", "decreasing"))    # 帰無仮説の種類（必ず片側検定）
{
        method <- "ヨンキー検定"
        data.name <- paste(deparse(substitute(x)), "by", deparse(substitute(g)))
        sgn <-  ifelse(match.arg(alternative) == "increasing", "<", ">")
        alternative <- sprintf("control %s= treat-1 %s= ... %s= treat-n", sgn, sgn, sgn)
        OK <- complete.cases(x, g)                                   # 欠損値を持つケースを除く
        x <- x[OK]
        g <- g[OK]
        ni <- table(g)                                                       # 各群のケース数
        gv <- as.numeric(names(ni))                                  # 群を表す数値
        a <- length(ni)                                                      # 群の個数
        n <- sum(ni)                                                 # 全ケース数
        sumSij <- sumTij <- 0
        for (i in 1:(a-1)) {
                di <- x[g==gv[i]]
                for (j in (i+1):a) {
                        dj <- x[g==gv[j]]
                        sumTij <- sumTij+sum(outer(di, dj, sgn))
                        sumSij <- sumSij+sum(outer(di, dj, "=="))
                }
        }
        tau <- table(x)
        V <- (n*(n-1)*(2*n+5)-sum(ni*(ni-1)*(2*ni+5))-sum(tau*(tau-1)*(2*tau+5)))/72 +
                sum(ni*(ni-1)*(ni-2))*sum(tau*(tau-1)*(tau-2))/(36*n*(n-1)*(n-2)) +
                sum(ni*(ni-1))*sum(tau*(tau-1))/(8*n*(n-1))             # 分散
        J0 <- sumTij+sumSij/2                                                # 検定統計量
        E <- (n^2-sum(ni^2))/4                                               # 期待値
        Z <- (abs(J0-E)-0.5*correct)/sqrt(V)                         # 標準化得点
        P <- pnorm(Z, lower.tail=FALSE)                                      # P 値
        return(structure(list(statistic=c(J=J0, "E(J)"=E, "V(J)"=V, "Z-value"=Z), p.value=P,
                method=method, alternative=alternative, data.name=data.name), class="htest"))
}


使用例

> x <- c(
+ 	153, 153, 152, 156, 158, 151, 151, 150, 148, 157,  # 第 1 群のデータ
+ 	158, 152, 152, 152, 151, 151, 157, 147, 155, 146,  # 第 2 群のデータ
+ 	153, 146, 138, 152, 140, 146, 156, 142, 147, 153,  # 第 3 群のデータ
+ 	137, 139, 141, 141, 143, 133, 147, 144, 151, 156   # 第 4 群のデータ
+ 	)
> g <- rep(1:4, each=10)

# 連続性の補正をしない場合

> Jonckheere(x, g, alternative="decreasing")

	ヨンキー検定

data:  x by g 
J = 446.5000, E(J) = 300.0000, V(J) = 1705.0776,
Z-value = 3.5479, p-value = 0.0001942
alternative hypothesis: control >= treat-1 >= ... >= treat-n 

 
# 連続性の補正をする場合

> Jonckheere(x, g, correct=TRUE, alternative="decreasing")

	ヨンキー検定

data:  x by g 
J = 446.5000, E(J) = 300.0000, V(J) = 1705.0776,
Z-value = 3.5357, p-value = 0.0002033
alternative hypothesis: control >= treat-1 >= ... >= treat-n