Julia に翻訳--158 一般化 Wilcoxon 検定

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

一般化 Wilcoxon 検定
http://aoki2.si.gunma-u.ac.jp/R/Gen-Wil.html

ファイル名: genwil.jl　　関数名: genwil

翻訳するときに書いたメモ

==========#

using Rmath

function genwil(group, event, time)
    n = length(group)
    o = sortperm(group)
    time = time[o]
    group = group[o]
    event = event[o]
    index, (na, nb) = table(group)
    W = 0
    for i = 1:na
        for j = na + 1:n
            W += getU(time[i], event[i], time[j], event[j])
        end
    end
    VarW = 0
    for i = 1:n
        temp = 0
        for j = 1:n
            temp += getU(time[i], event[i], time[j], event[j])
        end
        VarW += temp ^ 2
    end
    VarW = na * nb * VarW / n / (n - 1)
    Z = W / sqrt(VarW)
    P = pnorm(abs(Z), false) * 2
    println("W = $W,  V(W) = $VarW,  Z = $Z,  p value =$P")
    Dict(:W => W, :V => VarW, :Z => Z, :pvalue => P)
end

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

function table(x) # indices が少ないとき
    indices = sort(unique(x))
    counts = zeros(Int, length(indices))
    for i in indexin(x, indices)
        counts[i] += 1
    end
    return indices, counts
end

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.3126436781608,  Z = 1.7338126143696353,  p value =0.08295133672723795

Julia に翻訳--157 Cox-Mantel 検定

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

Cox-Mantel 検定
http://aoki2.si.gunma-u.ac.jp/R/Cox-Mantel.html

ファイル名: coxmantel.jl　　関数名: coxmantel

翻訳するときに書いたメモ

2 つのベクトルを返す関数で，後のベクトルの2要素を2変数に代入するときに
index, (nia, nib) = table(group) のようにすればよい

==========#

using Rmath, Printf

function coxmantel(group, event, time)
    len = length(group)
    len == length(event) == length(time) || error("length differ")
    indexx, indexy, tg = table(time, group * 10 + event)
    k = size(tg, 1)
    index, (nia, nib) = table(group)
    na = vcat(nia, (repeat([nia], k) - cumsum(tg[:, 1] + tg[:, 2]))[1:k - 1])
    nb = vcat(nib, (repeat([nib], k) - cumsum(tg[:, 3] + tg[:, 4]))[1:k - 1])
    m = tg[:, 2] + tg[:, 4]
    s = m .> 0
    m = m[s]
    r = (na + nb)[s]
    A = nb[s] ./ r
    U = sum(tg[:, 4]) - sum(m .* A)
    V = sum(m .* (r .- m) ./ (r .- 1) .* A .* (1 .- A))
    Z = U / sqrt(V)
    P = pnorm(abs(Z), false) * 2
    println("U = $U,  V(U) = $V,  Z = $Z,  p value = $P")
    Dict(:U => U, :V => V, :Z => Z, :pvalue => P)
end

function table(x) # indices が少ないとき
    indices = sort(unique(x))
    counts = zeros(Int, length(indices))
    for i in indexin(x, indices)
        counts[i] += 1
    end
    return indices, counts
end

function table(x, y) # 二次元
    indicesx = sort(unique(x))
    indicesy = sort(unique(y))
    counts = zeros(Int, length(indicesx), length(indicesy))
    for (i, j) in zip(indexin(x, indicesx), indexin(y, indicesy))
        counts[i, j] += 1
    end
    return indicesx, indicesy, counts
end

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];
coxmantel(group, event, time)
# U = 3.02979797979798,  V(U) = 2.715452780220283,  Z = 1.8386223884688957,  p value = 0.06597074594879868

Julia に翻訳--156 Cutler-Ederer 法による生命表

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

Cutler-Ederer 法による生命表
http://aoki2.si.gunma-u.ac.jp/R/ce-surv.html

ファイル名: cesurv.jl　　
関数名:　　 cesurv(time, event, width)
　　　　　  cesurv(ni, d, u, interval)

翻訳するときに書いたメモ

なんでもかんでもドット演算子にする必要はないということ

==========#

using Rmath, Printf, Plots

function cesurv(time, event, width)
    event .+= 1
    time = floor.(Int, time ./ width) .+ 1
    k = maximum(time)
    tbl = zeros(Int, k, 2)
    for i = 1:length(time)
        tbl[time[i], event[i]] += 1
    end
    cesurv(sum(tbl), vec(tbl[:, 2]), vec(tbl[:, 1]), 0:width:(k - 1) * width)
end

function cesurv(ni, d, u, interval)
    k = length(d)
    length(u) == k && length(interval) == k || error("length differ")
    n = repeat([ni], k) - cumsum(d + u)
    n = vcat(ni, n[1:end-1])
    np = n - u/2
    q = d ./ np
    p = 1 .- q
    P = cumprod(p)
    SE = P .* sqrt.(cumsum(q ./ (n - u/2 - d)))
    @printf("%8s %4s %3s %3s %5s %11s %11s %11s %11s\n",
            "interval", "n", "d", "u", "np", "q", "p", "P", "SE")
    for i = 1:k
        @printf("%8d %4d %3d %3d %5.1f %11.8f %11.8f %11.8f %11.8f\n",
                interval[i], n[i], d[i], u[i], np[i], q[i], p[i], P[i], SE[i])
    end
    interval = vcat(interval, 2interval[end] - interval[end-1])
    P = vcat(1, P)
    pyplot()
    plt = plot(interval, P, grid=false, tick_direction=:out, xlabel="interval",
               ylabel="P", ylimits=(-0.01, 1.01), label="")
    scatter!(interval, P, label="")
    display(plt)
end

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];
agroup = group .== 1;
cesurv(time[agroup], event[agroup], width)
# interval    n   d   u    np           q           p           P          SE
#        0   16   2   5  13.5  0.14814815  0.85185185  0.85185185  0.09668593
#        4    9   1   2   8.0  0.12500000  0.87500000  0.74537037  0.13068362
#        8    6   0   2   5.0  0.00000000  1.00000000  0.74537037  0.13068362
#       12    4   1   1   3.5  0.28571429  0.71428571  0.53240741  0.20275239
#       16    2   0   1   1.5  0.00000000  1.00000000  0.53240741  0.20275239
#       20    1   0   1   0.5  0.00000000  1.00000000  0.53240741  0.20275239

agroup = group .== 2;
cesurv(time[agroup], event[agroup], width)
# interval    n   d   u    np           q           p           P          SE
#        0   14   5   2  13.0  0.38461538  0.61538462  0.61538462  0.13493200
#        4    7   1   1   6.5  0.15384615  0.84615385  0.52071006  0.14359608
#        8    5   1   2   4.0  0.25000000  0.75000000  0.39053254  0.15591118
#       12    2   1   1   1.5  0.66666667  0.33333333  0.13017751  0.15904665

ni = 37
interval = 0:10:120
d = [3, 2, 3, 4, 1, 2, 3, 1, 0, 1, 0, 0, 0]
u = [8, 3, 1, 4, 0, 0, 0, 0, 0, 0, 0, 0, 1]
cesurv(ni, d, u, interval)
# interval    n   d   u    np           q           p           P          SE
#        0   37   3   8  33.0  0.09090909  0.90909091  0.90909091  0.05004381
#       10   26   2   3  24.5  0.08163265  0.91836735  0.83487941  0.06812545
#       20   21   3   1  20.5  0.14634146  0.85365854  0.71270193  0.08734827
#       30   17   4   4  15.0  0.26666667  0.73333333  0.52264808  0.10356244
#       40    9   1   0   9.0  0.11111111  0.88888889  0.46457607  0.10710681
#       50    8   2   0   8.0  0.25000000  0.75000000  0.34843206  0.10729149
#       60    6   3   0   6.0  0.50000000  0.50000000  0.17421603  0.08908649
#       70    3   1   0   3.0  0.33333333  0.66666667  0.11614402  0.07599690
#       80    2   0   0   2.0  0.00000000  1.00000000  0.11614402  0.07599690
#       90    2   1   0   2.0  0.50000000  0.50000000  0.05807201  0.05594695
#      100    1   0   0   1.0  0.00000000  1.00000000  0.05807201  0.05594695
#      110    1   0   0   1.0  0.00000000  1.00000000  0.05807201  0.05594695
#      120    1   0   1   0.5  0.00000000  1.00000000  0.05807201  0.05594695