Julia に翻訳--086 コクランの Q 検定

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

コクランの Q 検定
http://aoki2.si.gunma-u.ac.jp/R/Cochran-Q-test.html

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

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

==========#

using Rmath

function cochranqtest(x)
    k = size(x, 2)
    g = sum(x, dims = 1)
    l = sum(x, dims = 2)
    Q = ((k - 1) * (k * sum(g .^ 2) - sum(g)^2)) / (k * sum(l) - sum(l .^ 2))
    df = k - 1
    p = pchisq(Q, df, false)
    println("chisq. = $Q,  df = $df,  p value = $p")
end

x = [ 0 0 0
      0 0 0
      0 0 0
      0 0 1
      0 1 1
      0 1 1
      0 1 1
      1 1 1
      1 1 1
      1 1 1 ]

cochranqtest(x)
# chisq. = 6.5,  df = 2,  p value = 0.03877420783172201

Julia に翻訳--085 Linear-by-Linear 検定，Mantel の傾向検定

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

Linear-by-Linear 検定（Mantel の傾向検定）
http://aoki2.si.gunma-u.ac.jp/R/Linear_by_Linear.html

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

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

X を作るのが面倒

==========#

using Statistics, Rmath

function mantel(r, n, x=1:length(r))
    N = sum(n)
    R = sum(r)
    X = []
    for i = 1:length(x)
        append!(X, repeat([x[i]], r[i]))
    end
    for i = 1:length(x)
        append!(X, repeat([x[i]], n[i] - r[i]))
    end
    y = vcat(repeat([1], R), repeat([2], N-R));
    s = (N - 1) * cor(X, y)^2
    df = 1
    p = pchisq(s, df, false)
    println("chi2 = $s,  df = $df,  p value = $p")
    Dict(:chi2 => s, :df => df, :pvalue => p)
end

x = [10, 20, 30, 40, 50]
n = [30, 35, 47, 21, 45]
r = [2, 4, 14, 13, 39]

mantel(r, n, x)
# chi2 = 68.18749145474456,  df = 1,  p value = 1.486658203454285e-16

function wilcoxonscore(n)
    a = cumsum(n)
    a = vcat(0, a[1:end-1])
    a .+ (n .+ 1) ./ 2
end

score = wilcoxonscore(n)
# 5-element Array{Float64,1}:
#   15.5
#   48.0
#   89.0
#  123.0
#  156.0

mantel(r, n, score)
# chi2 = 67.7032747727581,  df = 1,  p value = 1.9004762464034343e-16

Julia に翻訳--084 コクラン・アーミテージ検定

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

コクラン・アーミテージ検定
http://aoki2.si.gunma-u.ac.jp/R/Cochran-Armitage.html

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

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

Julia のコーディング規則では変数名・関数名は小文字のみで，単語を _ で繋ぐこともしないというのは，可読性が落ちる。

==========#

using Rmath, Printf

function cochranarmitage(r, n, x=1:length(r))
    k = length(r)
    (length(n) == k && length(x) == k) || error("length differ")
    p = r ./ n
    R = sum(r)
    N = sum(n)
    pm = R / N
    xm = sum(n .* x) / N
    xx = x .- xm
    b = sum(n .* (p .- pm) .* xx) / sum(n .* xx .^ 2)
    a = pm - b * xm
    xt = b^2 * sum(n .* xx .^ 2) / (pm * (1 - pm))
    xh = N^2 * (sum(r .^ 2 ./ n) - R^2 / N) / R / (N - R)
    xq = xh - xt
    chisq = [xt, xq, xh]
    df = [1, k-2, k-1]
    pvalue = [pchisq(chisq0, df0, false) for (chisq0, df0) in zip(chisq, df)]
    rownames = ["トレンド", "直線からの乖離", "非一様性"]
    @printf("%26s  %6s  %10s\n", "カイ二乗値", "自由度", "P 値")
    for i = 1:3
        @printf("%14s  %10.7f  %6d  %10.7g\n", rownames[i], chisq[i], df[i], pvalue[i])
    end
    Dict(:chisq => chisq[1], :df => df[1], :pvalue => pvalue[1])
end

function wilcoxonscore(n)
    a = cumsum(n)
    a = vcat(0, a[1:end-1])
    a .+ (n .+ 1) ./ 2
end

x = [10, 20, 30, 40, 50]
n = [30, 35, 47, 21, 45]
r = [2, 4, 14, 13, 39]

cochranarmitage(r, n, x)

#                 カイ二乗値  自由度        P 値
#       トレンド  68.5727315       1  1.222834e-16
# 直線からの乖離   4.0144843       3   0.2599043
#       非一様性  72.5872158       4  6.449422e-15

cochranarmitage(r, n)

#                 カイ二乗値  自由度        P 値
#       トレンド  68.5727315       1  1.222834e-16
# 直線からの乖離   4.0144843       3   0.2599043
#       非一様性  72.5872158       4  6.449422e-15

ウィルコクソンスコアを使う場合

score = wilcoxonscore(n)

# 5-element Array{Float64,1}:
#   15.5
#   48.0
#   89.0
#  123.0
#  156.0

cochranarmitage(r, n, score)

#                 カイ二乗値  自由度        P 値
#       トレンド  68.0857792       1  1.565355e-16
# 直線からの乖離   4.5014366       3   0.2121622
#       非一様性  72.5872158       4  6.449422e-15

 

Julia に翻訳--083 マクネマー検定と拡張

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

マクネマー検定（拡張を含む）
http://aoki2.si.gunma-u.ac.jp/R/McNemar.html

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

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

二項検定の方が長い。

==========#

using Rmath

function mcnemar(tbl)
    nr, nc = size(tbl)
    n1 = sum([tbl[i, j] for i = 1:nr, j = 1:nc if i > j])
    n2 = sum([tbl[i, j] for i = 1:nr, j = 1:nc if i < j])
    p = binomtest(n1, n1 + n2, 0.5)
    println("n1 = $n1,  n2 = $n2,  p value = $p")
    Dict(:n1 => n1, :n2 => n2, :pvalue => p)
end

function binomtest(x, n, p = 0.5)
    if p == 0
        pvalue = (x == 0) + 0
    elseif p == 1
        pvalue = (x == n) + 0
    else
        relErr = 1 + 1e-07
        d = dbinom(x, n, p)
        m = n * p
        if x == m
            pvalue = 1
        elseif x < m
            i = ceil(Int, m):n
            y = sum(dbinom.(i, n, p) .<= d * relErr)
            pvalue = pbinom(x, n, p) + pbinom(n - y, n, p, false)
        else
            i = 0:floor(Int, m)
            y = sum(dbinom.(i, n, p) .<= d * relErr)
            pvalue = pbinom(y - 1, n, p) + pbinom(x - 1, n, p, false)
        end
    end
end

tbl = [ 11 3 2
         4 9 5
         1 2 13 ]

mcnemar(tbl)

# n1 = 7,  n2 = 10,  p value = 0.6290588378906251

Julia に翻訳--082 AIC による分割表の独立性の判定

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

AIC による分割表の独立性の判定
http://aoki2.si.gunma-u.ac.jp/R/AIC-independence.html

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

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

==========#

function aicindependence(x)
    ln(n) = sum([x * log(x) for x in n if x > 0])
    rt = sum(x, dims=2)
    ct = sum(x, dims=1)
    n = sum(x)
    lr, lc = size(x)
    AIC0 = 2 * (2 * ln(n) - ln(rt) - ln(ct) + lr + lc - 2)
    AIC1 = 2 * (ln(n) - ln(x) + lr * lc - 1)
    result = AIC0 < AIC1 ? "二要因は独立である" : "二要因は独立ではない"
    println("独立であるとしたときの AIC = $AIC0")
    println("独立でないとしたときの AIC = $AIC1")
    println("結論: $result")
end

aicindependence([749 445; 83 636])

# 独立であるとしたときの AIC = 5156.273864532803
# 独立でないとしたときの AIC = 4630.196189178783
# 結論: 二要因は独立ではない

Julia に翻訳--081 2×2 分割表のフィッシャーの正確確率検定

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

2×2 分割表のフィッシャーの正確確率検定
http://aoki2.si.gunma-u.ac.jp/R/my-fisher.html

ファイル名: myfisher.jl
関数名:    myfisher(x::Array{Int64, 2})
          myfisher(a::Int64, b::Int64, c::Int64, d::Int64)

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

同じ名前で 2 個の関数を定義。

==========#

using SpecialFunctions, Printf

function myfisher(x::Array{Int64, 2})
    myfisher(x[1,1], x[1,2], x[2,1], x[2,2])
end

function myfisher(a::Int64, b::Int64, c::Int64, d::Int64)
    mark(i) = i == 1 ? "@" : " "
    ct = [a + c, b + d]
    rt = [a + b, c + d]
    n = sum(ct)
    mx = min(rt[1], ct[1])
    mi = max(0, rt[1]+ct[1]-n)
    A = mi:mx
    B = rt[1] .- A
    C = ct[1] .- A
    D = ct[2] .- B
    nrow = length(A)
    Chisq = zeros(nrow);
    OddsRatio = similar(Chisq);
    Probability = similar(Chisq);
    for i = 1:nrow
        Chisq[i], OddsRatio[i], Probability[i] = stats(A[i], B[i], C[i], D[i])
    end
    CumProb1 = cumsum(Probability);
    CumProb2 = reverse(cumsum(reverse(Probability)));
    Chisq0, OddsRatio0, Probability0 = stats(a, b, c, d)
    Pearson = Chisq .>= Chisq0 - 1e-15
    Fisher = Probability .<= Probability0 + 1e-15
    OR = OddsRatio .>= OddsRatio0 - 1e-15
    pPearson = sum(Probability[Pearson])
    pFisher = sum(Probability[Fisher])
    pOR = sum(Probability[OR])
    println("   ピアソンの基準による p 値 = $pPearson")
    println("フィッシャーの基準による p 値 = $pFisher")
    println("   オッズ比を基準とした p 値 = $pOR")
    @printf("%3s %3s %3s %3s %8s %12s %10s %10s %1s%1s%1s %9s\n",
            "A", "B", "C", "D", "Chi.sq", "Probability", "OddsRatio",
            "Cum.Prob1", "P", "F", "O", "Cum.Prob2")
    for i = 1:nrow
        @printf("%3d %3d %3d %3d %8.3f %12.9f %10.3f %10.7f %1s%1s%1s %9.7f\n",
            A[i], B[i], C[i], D[i], Chisq[i], Probability[i], OddsRatio[i],
            CumProb1[i], mark(Pearson[i]), mark(Fisher[i]), mark(OR[i]), CumProb2[i])
    end
    Dict(:pPearson => pPearson, :pFisher => pFisher, :pOR => pOR)
end

function oddsratio(a, b, c, d)
    if a*b*c*d == 0
        a, b, c, d = a + 0.5, b + 0.5, c + 0.5, d + 0.5
    end
    res = a * d / (b * c)
    max(res, 1 / res)
end

lchoose(n, x) = logfactorial(n) - logfactorial(x) - logfactorial(n - x)

function stats(a, b, c, d)
    e, f, g, h = a + b, c + d, a + c, b + d
    n = e + f
    n * (a * d - b * c)^2 / (e * f * g * h), oddsratio(a, b, c, d),
        exp(lchoose(e, a) + lchoose(f, c) - lchoose(n, g))
end

myfisher([10 13; 16 61])
myfisher(10, 13, 16, 61)

#    ピアソンの基準による p 値 = 0.03529413313496242
# フィッシャーの基準による p 値 = 0.05508990606358157
#    オッズ比を基準とした p 値 = 0.05508990606358157
#   A   B   C   D   Chi.sq  Probability  OddsRatio  Cum.Prob1 PFO Cum.Prob2
#   0  23  26  51   10.495  0.000331755     24.184  0.0003318 @@@ 1.0000000
#   1  22  25  52    7.278  0.003815185     10.577  0.0041469 @@@ 0.9996682
#   2  21  24  53    4.649  0.019795773      4.755  0.0239427  @@ 0.9958531
#   3  20  23  54    2.606  0.061586849      2.840  0.0855296     0.9760573
#   4  19  22  55    1.151  0.128772503      1.900  0.2143021     0.9144704
#   5  18  21  56    0.282  0.192238950      1.350  0.4065410     0.7856979
#   6  17  20  57    0.000  0.212474629      1.006  0.6190156     0.5934590
#   7  16  19  58    0.305  0.177934419      1.336  0.7969501     0.3809844
#   8  15  18  59    1.198  0.114601829      1.748  0.9115519     0.2030499
#   9  14  17  60    2.677  0.057300915      2.269  0.9688528     0.0884481
#  10  13  16  61    4.743  0.022356750      2.933  0.9912096 @@@ 0.0311472
#  11  12  15  62    7.396  0.006818481      3.789  0.9980280 @@@ 0.0087904
#  12  11  14  63   10.636  0.001623448      4.909  0.9996515 @@@ 0.0019720
#  13  10  13  64   14.463  0.000300494      6.400  0.9999520 @@@ 0.0003485
#  14   9  12  65   18.877  0.000042928      8.426  0.9999949 @@@ 0.0000480
#  15   8  11  66   23.878  0.000004683     11.250  0.9999996 @@@ 0.0000051
#  16   7  10  67   29.465  0.000000384     15.314  1.0000000 @@@ 0.0000004
#  17   6   9  68   35.640  0.000000023     21.407  1.0000000 @@@ 0.0000000
#  18   5   8  69   42.402  0.000000001     31.050  1.0000000 @@@ 0.0000000
#  19   4   7  70   49.751  0.000000000     47.500  1.0000000 @@@ 0.0000000
#  20   3   6  71   57.686  0.000000000     78.889  1.0000000 @@@ 0.0000000
#  21   2   5  72   66.209  0.000000000    151.200  1.0000000 @@@ 0.0000000
#  22   1   4  73   75.318  0.000000000    401.500  1.0000000 @@@ 0.0000000
#  23   0   3  74   85.015  0.000000000   1000.429  1.0000000 @@@ 0.0000000
# Dict{Symbol,Float64} with 3 entries:
#   :pOR      => 0.0550899
#   :pPearson => 0.0352941
#   :pFisher  => 0.0550899

Julia に翻訳--080 対数尤度比検定，G2 検定，独立性の検定

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

対数尤度比（G2）に基づく独立性の検定
http://aoki2.si.gunma-u.ac.jp/R/g2.html

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

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

R だと
    ln <- function(n) sum(ifelse(n == 0, 0, n*log(n)))
でよいのだが，これを
    ln(n) = sum(n == 0 ? 0 : n .* log.(n))
などとしただけでは
    ln([2, 0, 4, 5])
は NaN になってしまうし，
    ln(n) = sum(n .== 0 ? 0 : n .* log.(n))
ではエラーになる。
結局
    ln(n) = sum([x * log(x) for x in n if x > 0])
ということかな。

==========#

using Rmath

function g2(mat; correct=false)
    ln(n) = sum([x * log(x) for x in n if x > 0])
    method = "対数尤度比に基づく独立性の検定（G-squared test）"
    n = sum(mat)
    rowsums = sum(mat, dims=2)
    colsums = sum(mat, dims=1)
    G2 = 2 * (ln(mat) - ln(rowsums) - ln(colsums) + ln(n))
    nrow, ncol = size(mat)
    df = (nrow - 1) * (ncol - 1)
    if correct == true
        method *= "連続性の補正"
        G2 /= 1 + (n * sum(1 ./ rowsums) - 1) * (n * sum(1 ./ colsums) - 1) /
                (6 * n * nrow * ncol)
    end
    P = pchisq(G2, df, false)
    println("$method\nG2 = $G2  自由度 = $df  p 値 = $P")
    Dict(:G2 => G2, :df => df, :pvalue => P)
end

g2([5 5 1; 6 4 4; 2 2 1])
# 対数尤度比に基づく独立性の検定（G-squared test）
# G2 = 1.8024287201126583  自由度 = 4  p 値 = 0.772037973903317

g2([4 5 2 0; 0 7 6 1; 1 0 3 1])
# 対数尤度比に基づく独立性の検定（G-squared test）
# G2 = 15.364591286599591  自由度 = 6  p 値 = 0.01760288650305146

g2([4 5 2 0; 0 7 6 1; 1 0 3 1], correct=true)
# 対数尤度比に基づく独立性の検定（G-squared test）連続性の補正
# G2 = 13.776490649307341  自由度 = 6  p 値 = 0.03223501559558125

Julia に翻訳--079 カイ二乗分布検定，独立性の検定，残差分析

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

カイ二乗分布を使用する独立性の検定と残差分析
http://aoki2.si.gunma-u.ac.jp/R/my-chisq-test.html

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

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

@printf の format として変数を使えないのが難

==========#

using Rmath, Printf

function mychisqtest(x)
    nr, nc = size(x)
    rowtotal = sum(x, dims=2)
    columntotal = sum(x, dims=1)
    n = sum(x)
    expectation = (rowtotal * columntotal) / n
    if any(expectation .< 1)
        println("警告: 期待値が 1 未満です")
    elseif sum(expectation .<= 5) / (nr * nc) > 0.2
        println("警告: 全セルの 20% 以上のセルの期待値が 5 以下です")
    end
    e = (x .- expectation) ./ sqrt.(expectation)
    d = e ./ sqrt.((1 .- rowtotal ./ n) * (1 .- columntotal ./ n))
    chi2 = sum(e .^ 2)
    df = (nr - 1) * (nc - 1)
    P = pchisq(chi2, df, false)
    P2 = [pnorm(abs(d[i, j]), false)*2 for i = 1:nr, j = 1:nc]
    println("カイ二乗値 = $chi2  自由度 = $df  p 値 = $P")
    printmatrix("標準化残差", e)
    printmatrix("調整された残差", d)
    printmatrix("P 値", P2)
    println()
    Dict(:chisq => chi2, :df => df, :pvalue => P,
         :standardizedresiduals => e, :adjustedresiduals => d, :P2 => P2)
end

function printmatrix(name, a)
    println(name)
    nr, nc = size(a)
    for i = 1:nr
        for j = 1:nc
            @printf(" %9.5f", a[i, j])
        end
        println()
    end
end

x = [4 5 2 0
     0 7 6 1
     1 0 3 1]

mychisqtest(x)

# 警告: 期待値が 1 未満です
# カイ二乗値 = 11.344332939787485 自由度 = 6 p 値 = 0.07829981227251993
# 標準化残差
#    1.60019   0.28604  -1.01246  -0.85635
#   -1.52753   0.59161   0.38252   0.06901
#    0.18257  -1.41421   0.86164   1.15470
# 調整された残差
#    2.20265   0.46402  -1.59862  -1.11382
#   -2.29129   1.04583   0.65817   0.09781
#    0.21909  -2.00000   1.18604   1.30931
# P 値
#    0.02762   0.64264   0.10991   0.26536
#    0.02195   0.29564   0.51043   0.92209
#    0.82658   0.04550   0.23561   0.19043