Julia に翻訳--124 ダネット法，平均値の多重比較

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

ダネットの方法による平均値の多重比較
http://aoki2.si.gunma-u.ac.jp/R/dunnett.html

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

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

mvtnorm の pmvt 関数は Julia にはないので，RCall で R を呼び出して計算させる。

==========#

using Statistics, Rmath, RCall, Printf

function dunnett(data, group)
    index, ni = table(group)
    a = length(ni)
    n = length(data)
    mi = [mean(data[group .== g]) for g in index]
    vi = [var(data[group .== g]) for g in index]
    Vw = sum(vi .* (ni .- 1)) / (n .- a)
    rho = getrho(ni)
    t = (abs.(mi .- mi[1]) ./ sqrt.(Vw .* (1 ./ ni .+ 1 / ni[1])))[2:a]
    p = pdunnett.(t, a, n - a, rho)
    p = vcat([float.(p[i]) for i = 1:a-1]...)
    @printf(" pair  %7s  %12s\n", "t value", "p value")
    for i = 1:a-1
        @printf("%2d:%2d  %7.5f  %12.7g\n", 1, i+1, t[i], p[i])
    end
    Dict(:t => t, :pvalue =>p)
end

func(x, y, ni1) = sqrt(x / (x + ni1) * y / (y + ni1))

function getrho(ni)
    k = length(ni)
    rho = func.(ni, transpose(ni), ni[1])
    [rho[i, i] = 0 for i =1:k]
    sum(rho[2:k, 2:k]) / (k - 2) / (k - 1)
end

function pdunnett(x, a, df, r)
    R"""
    library(mvtnorm)
    corr = diag($a - 1)
    corr[lower.tri(corr)] = $r
    1 - pmvt(lower = -$x, upper = $x, delta = numeric($a - 1), df = $df,
             corr = corr, abseps = 0.0001)[1]
    """
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 = [
  14, 15, 14, 16, 15, 17, 17,
  17, 16, 17, 16, 15, 18, 19, 15,
  18, 19, 20, 19, 17, 17,
  20, 21, 19, 20, 19, 22, 20,
  19, 20, 19, 17, 17, 17, 18 ]

group = [repeat([i], count) for (i, count) in enumerate([7, 8, 6, 7, 7])]
group = vcat(group...)

dunnett(data, group)

# pair  t value       p value
# 1: 2  1.85409     0.2156511
# 1: 3  4.18754  0.0008677534
# 1: 4  7.07368  1.02906e-07
# 1: 5  4.07273   0.001131645