一般化 Wilcoxon 検定

# 一般化 Wilcoxon 検定

打ち切り例のない場合の Wilcoxon 検定（ マン・ホイットニーの $U$ 検定と等価 ）の拡張である。<br>2 群の生存時間の長さに差があるかどうかを検定する。 

A 群と B 群のサンプルサイズを $n_a$，$n_b$，$n = n_a + n_b$ とする。
A 群，B 群の生存時間を $X_{i}$，$Y_{j}$，打ち切り例については $X_{i}'$，$Y_{j}'$ とする（ $i = 1, 2,  \dots , n_a\ ;\ j = 1, 2,  \dots , n_b$ ）。<br>
両群の生存時間を比較するとき，以下のような検定統計量 $W$ を定義する。

```math
    W = \sum_{i=1}^{n_a} \sum_{j=1}^{n_b} U_{ij}
```
        
ただし，$U_{ij}$ は以下のように定義する。

```math
    U_{ij} = \left \{
    \begin{align*}
        -1 &, 　X_i \lt Y_j\ または\ X_i \leqq Y_j' \\
        +1 &, 　X_i \gt Y_j\ または\ X_i' \geqq Y_j \\
         0 &, 　上記以外
    \end{align*}
    \right .
```

すなわち，生存時間の大小比較において，一方が打ち切り例のとき，例えば $X_{i}$ と $Y_{j}'$ の比較において，もし $X_{i} = Y_{j}'$ であった場合にも，$Y_{j}' = Y_{j} + \delta \gt Y_{j} \gt X_{i}$ であると考える点が「一般化」の意味するところである。

検定統計量 $W$ の平均値は 0 であり，分散 $Var ( W )$ は次式のようになる。

```math
        Var(W) = \frac{n_a\ n_b}{n\ (n-1)}\ \sum_{i=1}^{n} \left ( \sum_{j=1}^{n} H_{ij} \right ) ^2
```
        
ただし，$H_{ij}$ は上述の $U_{ij}$ の定義に準ずるが，2 群のデータををひとまとめにして，生存時間を $X_{i}$，$Y_{j}$，打ち切り例については $X_{i}'$，$Y_{j}'\ ( i, j = 1, 2, \dots , n)$ として，以下のように定義される。

```math   
    H_{ij} = \left \{
    \begin{align*}
        -1 &, 　X_i \lt Y_j\ または\ X_i \leqq Y_j' \\
        +1 &, 　X_i \gt Y_j\ または\ X_i' \geqq Y_j \\
         0 &, 　上記以外
    \end{align*}
    \right .
```
        
検定は，次式の $Z_o$ が近似的に標準正規分布にしたがうことを利用する。

```math
    Z_o = \frac{W}{\sqrt{Var(W)}}
```

$p$ 値は，標準正規分布において $|Z_o|$ より大きな値を取る確率を 2 倍したものになる（両側検定）。

```math
    p = 2 \times \Pr{(Z \gt |Z_o|)}
```

## プログラム例


```julia
using Distributions, DataFrames, Printf

outPval(P) = P < 0.00001 ? "< 0.00001" : @sprintf("= %.5f", P)

function getU(tx, sx, ty, sy)
    ((tx < ty && sx == 1 && sy == 1) || (tx <= ty && sx == 1 && sy == 0)) && return -1
    ((tx > ty && sx == 1 && sy == 1) || (tx >= ty && sx == 0 && sy == 1)) && return  1
    return 0
end

function genwil(group, event, time)
    df = DataFrame([group, time, event], [:group, :time, :event])
    gdf = groupby(df, :group)
    length(gdf) != 2 && return "error"
    na, nb = combine(gdf, nrow).nrow
    n = na + nb
    W = sum(getU(gdf[1].time[i], gdf[1].event[i], gdf[2].time[j], gdf[2].event[j])
        for i = 1:na, j = 1:nb)
    outer = [getU(time[i], event[i], time[j], event[j]) for i = 1:n, j = 1:n]
    VarW = na * nb * sum(sum(outer, dims=1) .^ 2) / n / (n - 1)
    Zo = W / sqrt(VarW)
    pvalue = 2ccdf(Normal(), abs(Zo))
    @printf("W = %d,  V(W) = %.5f,  Z = %.5f,  p value %s", W, VarW, Zo, outPval(pvalue))
    Dict(:W => W, :V => VarW, :Z => Zo, :pvalue => pvalue)
end;
```

## 使用例


```julia
group = [1, 1, 2, 2, 1, 1, 2, 2, 1, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 2, 1, 2, 1, 2, 1, 2, 2, 1, 1, 1];
event = [1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0];
time = [2, 20, 5, 1, 3, 17, 2, 3, 15, 14, 12, 13, 11, 11, 10, 8, 8, 3, 7, 3, 6, 2, 5, 4, 2, 3, 1, 3, 2, 1];

genwil(group, event, time); # W = 63,  V(W) = 1320.31264,  Z = 1.73381,  p value = 0.08295
```

    W = 63,  V(W) = 1320.31264,  Z = 1.73381,  p value = 0.08295


```julia
a = genwil(group, event, time);
a[:pvalue]
```

    W = 63,  V(W) = 1320.31264,  Z = 1.73381,  p value = 0.08295
    0.08295133672723796