目的

マン・ホイットニーの U 検定を行う
（R には，等価な検定を行う関数 wilcox.test が用意されている）

使用法

二つのデータベクトルを与えるとき
　U.test(x, y, correct=FALSE)
　U.test(x, y)

2 行 k 列の分割表を与えるとき
　U.test(x, correct = FALSE)
　U.test(x)

引数

x        第一群の観測値ベクトルまたは，分割表データ（y=NULL)
y        第二群の観測値ベクトル
correct  連続性の修正を行うかどうかを示す
         correct=FALSE のときには，連続性の修正を行わない
         correct=TRUE のときには，連続性の修正を行う
         この引数が省略されたときは，連続性の修正を行う

ソース

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

# マン・ホイットニーの U 検定
U.test <- function(  x,                              # 第一群の観測値ベクトルまたは，分割表データ（y=NULL)
                        y = NULL,                       # 第二群の観測値ベクトル
                        correct = TRUE)                 # 連続性の補正を行うかどうか
{
        method <- "マン・ホイットニーの U 検定"
        if (is.null(y)) {                               # 2 × C 行列の分割表として与えられたとき
                if (nrow(x) != 2) stop("2 x C matrix is expected.")
                data.name <- paste(deparse(substitute(x)), "as 2 by C matrix")
                nc <- ncol(x)                                # カテゴリー数
                y <- x[2,]                           # 第二群の度数分布
                x <- x[1,]                           # 第一群の度数分布
                tie <- x+y                           # 合計した度数分布（同順位）
                n1 <- sum(x)                         # 第一群のサンプルサイズ
                n2 <- sum(y)                         # 第二群のサンプルサイズ
                n <- n1+n2                           # 合計したサンプルサイズ
                rj <- c(0, cumsum(tie)[-nc])+(tie+1)/2       # カテゴリーに属するものの順位
                U1 <- n1*n2+n1*(n1+1)/2-sum(x*rj)    # 検定統計量
        }
        else {                                          # 2 つのデータベクトルとして与えられたとき
                data.name <- paste(deparse(substitute(x)), "and", deparse(substitute(y)))
                x <- x[!is.na(x)]                    # 欠損値を持つケースを除く
                y <- y[!is.na(y)]                    # 欠損値を持つケースを除く
                n1 <- length(x)                              # 第一群のサンプルサイズ
                n2 <- length(y)                              # 第二群のサンプルサイズ
                n <- n1+n2                           # 合計したサンプルサイズ
                xy <- c(x, y)                                # 両群のデータを結合したもの
                r <- rank(xy)                                # 順位のベクトル
                U1 <- n1*n2+n1*(n1+1)/2-sum(r[1:n1]) # 検定統計量
                tie <- table(r)                              # 同順位の集計
        }
        U <- min(U1, n1*n2-U1)                               # U 統計量
        V <- n1*n2*(n^3-n-sum(tie^3-tie))/12/(n^2-n) # 同順位を考慮した分散
        E <- n1*n2/2                                 # 期待値
        Z <- (abs(U-E)-ifelse(correct, 0.5, 0))/sqrt(V)      # Z 値
        P <- pnorm(Z, lower.tail=FALSE)*2            # 両側 P 値
        return(structure(list(statistic=c(U=U, "E(U)"=E, "V(U)"=V, "Z-value"=Z), p.value=P,
                method=method, data.name=data.name), class="htest"))
}


使用例

二つのデータベクトルを与えるとき

> A <- rep(1:4, c(9, 12, 6, 3)) # 右と同じ　c(rep(1, 9), rep(2, 12), rep(3, 6), rep(4, 3))
> B <- rep(1:4, c(4, 9, 11, 5)) # 右と同じ　c(rep(1, 4), rep(2, 9), rep(3, 11), rep(4, 5))
> U.test(A, B, correct=FALSE)

	マン・ホイットニーの U 検定

data:  A and B 
U = 310.5000, E(U) = 435.0000, V(U) = 3993.5593, Z-value = 1.9701, p-value = 0.04883

以下と比較のこと
> wilcox.test(A, B, correct=FALSE)

	Wilcoxon rank sum test

data:  A and B
W = 310.5, p-value = 0.04883
alternative hypothesis: true mu is not equal to 0

Warning message:
Cannot compute exact p-value with ties in: wilcox.test.default(A, B, correct = FALSE)

分割表を与えるとき

        C1  C2  C3　合計
第一群  10   5   1　　16
第二群　 4   7   3    14

> x <- matrix(c(10, 5, 1, 4, 7, 3), byrow=TRUE, nrow=2)
> U.test(x, correct=FALSE)

	マン・ホイットニーの U 検定

data:  x as 2 by C matrix 
U = 70.0000, E(U) = 112.0000, V(U) = 481.9862, Z-value = 1.9131, p-value = 0.05574

 解説ページ