Julia に翻訳--128 ボンフェローニ， ホルム，シェイファー，ホランド・コペンハーバー，平均値の多重比較

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

ボンフェローニの方法による平均値の多重比較
ホルムの方法による平均値の多重比較
シェイファーの方法による平均値の多重比較
ホランド・コペンハーバーの方法による平均値の多重比較
http://aoki2.si.gunma-u.ac.jp/R/Bonferroni.html

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

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

==========#

using Statistics, Rmath, Combinatorics, Printf

function bonferroni(data, group; alpha = 0.05, method = "Bonferroni") # "Holm", "Shaffer", "Holland-Copenhaver"
    index, n = table(group)
    a = length(n)
    k = binomial(a, 2)
    m = [mean(data[group .== g]) for g in index]
    v = [var(data[group .== g]) for g in index]
    phi = length(data) - a
    V = sum((n .- 1) .* v) / phi
    pair = combinations(1:a, 2)
    t  = [abs(m[i] - m[j]) / sqrt(V / n[i] + V / n[j]) for (i, j) in pair]
    P = pt.(t, phi, false) .* 2
    pairs  = [i for i in pair]
    result2 = hcat(pairs, t, P);
    if method == "Bonferroni"
        Pthreshold = repeat([alpha / k], k)
        judge = (P .<= Pthreshold) .+ 0
    else
        result2 = result2[sortperm(result2[:, 3]), :]
        if method != "Holm" && a > 9
            error("Too many groups Use Holm")
            method = "Holm"
        end
        if method == "Holm"
            Pthreshold = alpha ./ (k:-1:1)
        else
            tbl = vcat([3, 1, 1, 6, repeat([3], 3), 2, 1, 10, repeat([6], 4),
                        4, 4:-1:1, 15, repeat([10], 5), repeat([7], 3),
                        6, 4, 4:-1:1, 21, repeat([15], 6), repeat([11], 4),
                        10, 9, 7, 7:-1:1, 28, repeat([21], 7), repeat([16], 5),
                        15, 13, 13:-1:1, 36, repeat([28], 8), repeat([22], 6),
                        21, repeat([18], 3), 16, 16, 15, 13, 13:-1:1]...)
            h = tbl[(a * (a - 1) * (a - 2) ÷ 6):((a ^ 3 - a) ÷ 6 - 1)]
            Pthreshold = method == "Shaffer" ? alpha ./ h : 1 .- (1 - alpha) .^ (1 ./ h)
        end
        judge = (cumprod(result2[:, 3] .<= Pthreshold) .!= 0) .+ 0
        if method == "Holm" || method == "Shaffer" || method == "Holland-Copenhaver"
            minpos = argmin(judge)
            if judge[minpos] == 0
                judge[minpos:k] .= 2
            end
        end
    end
    judge2 = [["not significant", "significant", "suspended"][j+1] for j in judge]
    result2 = hcat(result2, Pthreshold, judge2)
    @printf(" pair  %7s  %12s  %12s  %7s\n", "t value", "p value", "p threshold", "judge")
    for i = 1:k
         @printf("%2d:%2d  %7.5f  %12.7g  %12.7g  %7s\n",
                 result2[i, 1][1], result2[i, 1][2], result2[i, 2],
                 result2[i, 3], Pthreshold[i], judge2[i])
    end
    Dict(:alpha => alpha, :noftests => k, :result2 => result2)
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

data = [
  10.7, 9.7, 8.5, 9.4, 8.8, 8.4, 10.6,
  8.1, 8.3, 8.7, 6.9, 5.7, 9.5, 6.7,
  7.9, 7.5, 7.4, 9.2, 5.7, 8.3, 9.7,
  6.2, 7.1, 5.5, 4.7, 6.3, 6.9, 7.5
]
group = vcat([repeat([i], 7) for i = 1:4]...)

bonferroni(data, group) # method = "Bonferroni"
# pair  t value       p value   p threshold    judge
# 1: 2  2.83517   0.009148816   0.008333333  not significant
# 1: 3  2.41686    0.02362029   0.008333333  not significant
# 1: 4  5.08935  3.315126e-05   0.008333333  significant
# 2: 3  0.41830     0.6794449   0.008333333  not significant
# 2: 4  2.25419    0.03358738   0.008333333  not significant
# 3: 4  2.67249    0.01331927   0.008333333  not significant

bonferroni(data, group, method = "Holm")
# pair  t value       p value   p threshold    judge
# 1: 4  5.08935  3.315126e-05   0.008333333  significant
# 1: 2  2.83517   0.009148816          0.01  significant
# 3: 4  2.67249    0.01331927        0.0125  suspended
# 1: 3  2.41686    0.02362029    0.01666667  suspended
# 2: 4  2.25419    0.03358738         0.025  suspended
# 2: 3  0.41830     0.6794449          0.05  suspended

bonferroni(data, group, method = "Shaffer")
# pair  t value       p value   p threshold    judge
# 1: 4  5.08935  3.315126e-05   0.008333333  significant
# 1: 2  2.83517   0.009148816    0.01666667  significant
# 3: 4  2.67249    0.01331927    0.01666667  significant
# 1: 3  2.41686    0.02362029    0.01666667  suspended
# 2: 4  2.25419    0.03358738         0.025  suspended
# 2: 3  0.41830     0.6794449          0.05  suspended

bonferroni(data, group, method = "Holland-Copenhaver")
# pair  t value       p value   p threshold    judge
# 1: 4  5.08935  3.315126e-05   0.008512445  significant
# 1: 2  2.83517   0.009148816    0.01695243  significant
# 3: 4  2.67249    0.01331927    0.01695243  significant
# 1: 3  2.41686    0.02362029    0.01695243  suspended
# 2: 4  2.25419    0.03358738    0.02532057  suspended
# 2: 3  0.41830     0.6794449          0.05  suspended