Julia に翻訳--212 二群の比率の差の検定

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

二群の比率の差の検定
http://aoki2.si.gunma-u.ac.jp/lecture/Hiritu/diff-p-test-r.html

ファイル名: proptest.jl
関数名:    proptest(x::Vector{Int64}, n::Vector{Int64}; correction=true)
          proptest(x::Array{Int64,2}; correction=true)

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

K 群の比率の差の検定は，2xK 分割表での独立性の検定と同じである。ということで，chisqtest2 関数を下請けにしている。
下のプログラムは，2群の比率の差の検定だけでなく，K群の比率の差の検定にも対応している。

==========#

using Rmath, Printf

function proptest(x::Vector{Int64}, n::Vector{Int64}; correction=true)
    chisqtest2(hcat(x, n .- x); correction)
end

function proptest(x::Array{Int64,2}; correction=true)
    chisqtest2(x; correction)
end

function chisqtest2(x; correction=true)
    n, m = size(x)
    rowsum = sum(x, dims=2)
    colsum = sum(x, dims=1)
    expectation =  (rowsum * colsum) ./ sum(x)
    if sum(expectation .< 1) >= 1 ||sum(expectation .< 5) >= 0.2n*m
        println("Warning message: Chi-squared approximation may be incorrect")
    end
    difference = abs.(x - expectation)
    yates = 0
    correction && n == 2 && m == 2 && (yates = min(0.5, minimum(difference)))
    chisq = sum(((difference .- yates) .^ 2) ./ expectation)
    df =(n - 1) * (m - 1)
    pvalue = pchisq(chisq, df, false)
    @printf("chisq = %.5f,  df = %d,  p value = %.5f", chisq, df, pvalue)
end

smokers  = [83, 90, 129, 70]
patients = [86, 93, 136, 82];
proptest(smokers, patients) # chisq = 12.60041,  df = 3,  p value = 0.00559

data = [83 90 129 70
         3  3   7 12];
proptest(data)              # chisq = 12.60041,  df = 3,  p value = 0.00559

proptest([145, 157], [300, 250], correction=false) # chisq = 11.52663,  df = 1,  p value = 0.00069
proptest([145, 157], [300, 250])                   # chisq = 10.94973,  df = 1,  p value = 0.00094