#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。
ファイル名: 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
