Julia に翻訳--107 二要因の分散分析，ASB タイプ，SPFp・q デザイン，混合計画

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

二要因の分散分析（ASB タイプ；SPFp・q デザイン；混合計画）
http://aoki2.si.gunma-u.ac.jp/R/ASB.html

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

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

apply(X, MARGIN, FUN) との関係が明示されていれば，移植は簡単。

==========#

using Statistics, Rmath, Printf

function asb(data)
    n, Nb, Na = size(data)
    nm1 = n - 1
    grandmean = mean(data)
    ea = sum((mean(data, dims=[1,2]) .- grandmean) .^ 2) * Nb * n
    diffobj = sum((mean(data, dims=2) .- grandmean) .^ 2) * Nb - ea
    eb = sum((mean(data, dims=[1,3]) .- grandmean) .^ 2) * Na * n
    cross = sum((mean(data, dims=1) .- grandmean) .^ 2) * n - ea - eb
    err = sum(var(data, dims=1) * nm1) - diffobj
    SS = [ea, diffobj, eb, cross, err]
    dfa = Na - 1
    dfb = Nb - 1
    df = [dfa, Na * nm1, dfb, dfa * dfb,  Na * nm1 * dfb]
    MS = SS ./ df
    F = fill(NaN, 5)
    P = fill(NaN, 5)
    F[[1, 3, 4]] = MS[[1, 3, 4]] ./ MS[[2, 5, 5]]
    P[[1, 3, 4]] = pf.(F[[1, 3, 4]], df[[1, 3, 4]], df[[2, 5, 5]], false)
    name = ["Factor A", "S", "Factor B", "AxS", "SxB"]
    @printf("%8s %11s  %3s  %11s  %9s  %7s\n",
            "", "SS", "df", "MS", "F value", "P value")
    for i = 1:5
        if i in [1, 3, 4]
            @printf("%8s %11.8g  %3d  %11.8g  %9.5f  %7.5f\n",
                    name[i], SS[i], df[i], MS[i], F[i], P[i])
        else
            @printf("%8s %11.8g  %3d  %11.8g\n",
                    name[i], SS[i], df[i], MS[i])
        end
    end
    Dict(:name => name, :SS => SS, :df => df, :MS => MS, :F => F, :P => P)
end

data = [4, 4, 5,  3, 4, 4,  6, 6, 4,  3, 4, 3,  5, 6, 7,  5, 4, 4];
data = reshape(data, 3, 3, 2);
asb(data)
#                   SS   df           MS    F value  P value
# Factor A 0.055555556    1  0.055555556    0.50000  0.51852
#        S  0.44444444    4   0.11111111
# Factor B           4    2            2    2.32258  0.16020
#      AxS   11.111111    2    5.5555556    6.45161  0.02145
#      SxB   6.8888889    8   0.86111111

Julia に翻訳--106 二要因の分散分析，SAB タイプ，RBFpq デザイン，被検者内計画

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

二要因の分散分析（SAB タイプ；RBFpq デザイン；被検者内計画）
http://aoki2.si.gunma-u.ac.jp/R/SAB.html

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

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

R の apply(X, MARGIN, FUN) と Julia の FUN(X, dims) の関係
reshape は不要な場合もある

x = reshape(1:24, (2, 3, 4))            # == x = array(1:24, dim=c(2,3,4))
reshape(mean(x, dims=1), 3, 4)          # == apply(x, 2:3, mean)
reshape(mean(x, dims=2), 2, 4)          # == apply(x, c(1, 3), mean)
reshape(mean(x, dims=3), 2, 3)          # == apply(x, 1:2, mean)
reshape(mean(x, dims=[1,2]), 4, 1, 1)   # == apply(x, 3, mean)
reshape(mean(x, dims=[2,3]), 2, 1, 1)   # == apply(x, 1, mean)
reshape(mean(x, dims=[1,3]), 3, 1, 1)   # == apply(x, 2, mean)
reshape(mean(x, dims=[1,2,3]), 1, 1, 1) # == mean(x)

==========#

using Rmath, Statistics, Printf

function sab(data)
    Nb, Na, N = size(data)
    Xbar = vec(mean(data, dims=3))
    SD = vec(std(data, dims=3)) .* sqrt((N - 1) / N)
    v1 = mean(data)
    v2 = sum((mean(data, dims=[1,3]) .- v1) .^ 2) * Nb * N
    v3 = sum((mean(data, dims=[2,3]) .- v1) .^ 2) .* Na .* N
    v4 = sum((mean(data, dims=[3]) .- v1) .^ 2) * N
    v42 = v4-v2-v3
    v5 = sum(SD .^ 2) * N
    v6 = sum(data) .^ 2 / (N * Na * Nb)
    v62 = sum(sum(data, dims=[1,2]) .^ 2) / (Na * Nb) - v6
    v7 = sum(sum(data, dims=1) .^ 2) / Nb - v6 - v62 - v2
    v8 = sum(sum(data, dims=2) .^ 2) / Na - v6 - v62 - v3
    v9 = v5 - v62 - v7 - v8
    SS = [v62, v2, v7, v3, v8, v42, v9]
    dfs = N - 1
    dfa = Na - 1
    dfb = Nb - 1
    df = [dfs, dfa, dfs * dfa, dfb, dfs * dfb, dfa * dfb, dfs * dfa * dfb]
    MS = SS ./ df
    F = fill(NaN, 7)
    P = fill(NaN, 7)
    F[[2, 4, 6]] = MS[[2, 4, 6]] ./ MS[[3, 5, 7]]
    P[[2, 4, 6]] = pf.(F[[2, 4, 6]], df[[2, 4, 6]], df[[3, 5, 7]], false)
    @printf("%5s %11s  %3s  %11s  %9s  %7s\n",
            "", "SS", "df", "MS", "F value", "P value")
    name = ["S", "A", "SxA", "B", "SxB", "AxB", "SxAxB"]
    for i = 1:7
        if i % 2 == 1
            @printf("%5s %11.8g  %3d  %11.8g\n", name[i], SS[i], df[i], MS[i])
        else
            @printf("%5s %11.8g  %3d  %11.8g  %9.5f  %7.5f\n",
                    name[i], SS[i], df[i], MS[i], F[i], P[i])
        end
    end
    Dict(:name => name, :SS => SS, :df => df, :MS => MS, :F => F, :P => P)
end

data = [4, 3, 6, 3, 5, 5, 4, 4, 6, 4, 6, 4, 5, 4, 4, 3, 7, 4];
data = reshape(data, 3, 2, 3)
sab(data)
#                SS   df           MS    F value  P value
#     S  0.33333333    2   0.16666667
#     A 0.055555556    1  0.055555556    1.00000  0.42265
#   SxA  0.11111111    2  0.055555556
#     B           4    2            2    1.71429  0.28994
#   SxB   4.6666667    4    1.1666667
#   AxB   11.111111    2    5.5555556   10.00000  0.02778
# SxAxB   2.2222222    4   0.55555556

Julia の小ネタ--012 R の apply(X, MARGIN, FUN) と Julia の FUN(X, dims) の関係

R の apply(X, MARGIN, FUN) と Julia の FUN(X, dims) の関係
reshape は不要な場合もある

x = reshape(1:24, (2, 3, 4))            # == x = array(1:24, dim=c(2,3,4))
reshape(mean(x, dims=1), 3, 4)          # == apply(x, 2:3, mean)
reshape(mean(x, dims=2), 2, 4)          # == apply(x, c(1, 3), mean)
reshape(mean(x, dims=3), 2, 3)          # == apply(x, 1:2, mean)
reshape(mean(x, dims=[1,2]), 4, 1, 1)   # == apply(x, 3, mean)
reshape(mean(x, dims=[2,3]), 2, 1, 1)   # == apply(x, 1, mean)
reshape(mean(x, dims=[1,3]), 3, 1, 1)   # == apply(x, 2, mean)
reshape(mean(x, dims=[1,2,3]), 1, 1, 1) # == mean(x)

Julia に翻訳--105 二元配置分散分析

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

二元配置分散分析
http://aoki2.si.gunma-u.ac.jp/R/twoway-anova.html

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

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

元のプログラムにある二次元配列 mab は不要だ。

==========#

using Statistics, Rmath, Printf

function twowayanova(x, a, b)
    indexa, indexb, tbl = table(a, b)
    n = length(x)
    dft = n-1
    na, nb = size(tbl)
    dfa = na - 1
    dfb = nb - 1
    gm = mean(x)
    MSt = var(x)
    sst = dft * MSt
    ma = [mean(x[a .== i]) for i in indexa]
    mb = [mean(x[b .== j]) for j in indexb]
    if any(tbl .== 0)
        error("いくつかの水準の組合せにおいて，データがないセルがあります")
    elseif all(tbl .== 1)
        ssa = nb * sum((ma .- gm) .^ 2)
        ssb = na * sum((mb .- gm) .^ 2)
        sse = sst-ssa-ssb
        dfe = dft-dfa-dfb
        ss = [ssa, ssb, sse, sst]
        df = [dfa, dfb, dfe, dft]
        ms = ss ./ df
        f = ms ./ ms[3]
        p = pf.(f, df, df[3], false)
        f[3:4] .= NaN
        p[3:4] .= NaN
        result = hcat(ss, df, ms, f, p)
        @printf("%12s  %6s  %10s  %10s  %10s\n", "SS", "df", "MS", "F value", "P value")
        @printf("a %10.7g  %6d  %10.7g  %10.7f  %10.7f\n", ssa, dfa, ms[1], f[1], p[1])
        @printf("b %10.7g  %6d  %10.7g  %10.7f  %10.7f\n", ssb, dfb, ms[2], f[2], p[2])
        @printf("e %10.7g  %6d  %10.7g\n", sse, dfe, ms[3])
        @printf("T %10.7g  %6d  %10.7g\n", sst, dft, ms[4])
        Dict(:result => result)
    elseif all(tbl .>= 2) && hirei(tbl)
        sse = 0
        for i in 1:na, j in 1:nb
            selected = x[(a .== indexa[i]) .& (b .== indexb[j])]
            sse += sum((selected .- mean(selected)) .^ 2)
        end
        ssa = sum(sum(tbl, dims=2) .* (ma .- gm) .^ 2)
        ssb = sum(transpose(sum(tbl, dims=1)) .* (mb .- gm) .^ 2)
        ss = [ssa, ssb, sst-ssa-ssb-sse, sse, sst]
        dfab = dfa * dfb
        df = [dfa, dfb, dfab, dft-dfa-dfb-dfab, dft]
        ms = ss ./ df
        f = ms ./ ms[4]
        p = pf.(f, df, df[4], false)
        f[4:5] .= NaN
        p[4:5] .= NaN
        result = hcat(ss, df, ms, f, p)
        names = ["a", "b", "a*b", "e", "T"]
        output("Model-I", names, ss, df, ms, f, p)
        f[1:2] = ms[1:2] ./ ms[3]
        p[1:2] = pf.(f[1:2], df[1:2], df[3], false)
        result2 = hcat(ss, df, ms, f, p)
        println()
        output("Model-II", names, ss, df, ms, f, p)
        Dict(:model1 => result, :model2 => result2)
    else
        d1 = makedummy(a, indexa);
        d2 = makedummy(b, indexb);
        d3 = crossdummy(d1, d2);
        rabab = mreg(hcat(d3, x))
        rab = mreg(hcat(d1, d2, x))
        ra = mreg(hcat(d1, x))
        rb = mreg(hcat(d2, x))
        ss = [rab[1]-rb[1], rab[1]-ra[1], rabab[1]-rab[1], rabab[2], rabab[3]]
        df = [dfa, dfb, dfa*dfb, n-na*nb, n-1]
        ms = ss ./ df
        f = ms ./ ms[4]
        p = pf.(f, df, df[4], false)
        f[4:5] .= NaN
        p[4:5] .= NaN
        result = hcat(ss, df, ms, f, p)
        names = ["a", "b", "a*b", "e", "T"]
        output(names, ss, df, ms, f, p)
        Dict(:result => result)
    end
end

function hirei(tbl)
    for i in 2:size(tbl, 1)
        temp = tbl[i, :] ./ tbl[1, :]
        all(temp .== temp[1]) || return false
    end
    return true
end

function output(model, names, ss, df, ms, f, p)
    println("$model")
    @printf("%14s  %6s  %10s  %10s  %10s\n", "SS", "df", "MS", "F value", "P value")
    for i in 1:3
        @printf("%4s%10.7g  %6d  %10.7g  %10.7f  %10.7f\n",
                names[i], ss[i], df[i], ms[i], f[i], p[i])
    end
    for i in 4:5
        @printf("%4s%10.7g  %6d  %10.7g\n", names[i], ss[i], df[i], ms[i])
    end
end

function mreg(dat) # 回帰分析を行い，回帰変動，誤差変動，全変動を返す関数
    r = cor(dat)
    betas = r[1:end-1, 1:end-1] \ r[1:end-1, end]
    SS = cov(dat, corrected=false) .* size(dat, 1)
    prop = [SS[i, i] for i in 1:size(SS, 1)]
    Sr = sum(SS[1:end-1, end] .* (betas ./ sqrt.(prop[1:end-1] ./ prop[end])))
    St = SS[end, end]
    (Sr, St - Sr, St)  # 回帰変動，誤差変動，全変動
end

function makedummy(x, y)
    n = length(x)
    k = length(y)
    d = zeros(n, k)
    for (i, j) = enumerate(indexin(x, y))
        d[i, j] = 1
    end
    d[:, 2:end]
end

function crossdummy(a, b)
    n1, k = size(a)
    n2, l = size(b)
    n1 == n2 || error("data error")
    d = zeros(n1, k * l)
    for i = 1:l
        for j = 1:k
            d[:, j + (i - 1) * k] = a[:, j] .* b[:, i]
        end
    end
    hcat(a, b, d)
end

function table(x, y)
    indicesx = sort(unique(x))
    indicesy = sort(unique(y))
    counts = zeros(Int, length(indicesx), length(indicesy))
    for (i, j) in zip(indexin(x, indicesx), indexin(y, indicesy))
        counts[i, j] += 1
    end
    return indicesx, indicesy, counts
end

a = [1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3];
b = repeat(1:4, 3);
x = [9, 17, 12, 16, 1, 21, 16, 11, 7, 19, 6, 9];
twowayanova(x, a, b)
#           SS      df          MS     F value     P value
# a       21.5       2       10.75   0.6592845   0.5510300
# b   268.6667       3    89.55556   5.4923339   0.0371924
# e   97.83333       6    16.30556
# T        388      11    35.27273

x = [24.8, 23.9, 24.1, 28.8, 22.6, 28.0, 26.4, 27.4, 29.4, 30.0, 28.7, 
    29.2, 25.0, 26.6, 27.9, 28.5, 27.1, 25.2, 27.9, 29.2, 26.7, 25.0, 
    29.9, 29.4, 27.5, 32.5, 29.5, 26.3, 28.2, 31.8, 30.2, 31.7, 29.2, 
    30.3, 30.9, 28.4, 29.8, 26.7, 30.7, 28.5, 26.5, 26.7, 31.7, 27.2, 
    25.5, 29.9, 27.3, 28.8, 28.5, 25.7, 28.7, 31.3, 29.4, 29.8, 30.3, 
    29.6, 31.7, 33.6, 32.0, 34.3];
a = [1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 1, 1, 1, 2, 2, 2, 3, 3, 
    3, 4, 4, 4, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 1, 1, 1, 2, 2, 
    2, 3, 3, 3, 4, 4, 4, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4];
b = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 
    2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 
    4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5];
twowayanova(x, a, b)
#  Model-I
#            SS      df          MS     F value     P value
#   a  51.39733       3    17.13244   4.9370667   0.0052039
#   b  106.2873       4    26.57183   7.6572211   0.0001115
# a*b  52.85267      12    4.404389   1.2692154   0.2738926
#   e  138.8067      40    3.470167
#   T   349.344      59    5.921085
#
#  Model-II
#            SS      df          MS     F value     P value
#   a  51.39733       3    17.13244   3.8898573   0.0373908
#   b  106.2873       4    26.57183   6.0330352   0.0067194
# a*b  52.85267      12    4.404389   1.2692154   0.2738926
#   e  138.8067      40    3.470167
#   T   349.344      59    5.921085

a = [1, 1, 1, 1, 2, 2, 2, 2, 2, 2]
b = [1, 1, 2, 2, 1, 1, 2, 2, 2, 2]
x = [17, 16, 25, 22, 18, 26, 34, 30, 34, 30]
twowayanova(x, a, b)
#            SS      df          MS     F value     P value
# a    121.4405       1    121.4405  13.7479784   0.0099953
# b    177.1905       1    177.1905  20.0592992   0.0041979
# a*b  5.142857       1    5.142857   0.5822102   0.4743687
# e          53       6    8.833333
# T       415.6       9    46.17778

Julia に翻訳--104 クラスカル・ウォリス検定と多重比較

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

クラスカル・ウォリス検定（plus 多重比較）
http://aoki2.si.gunma-u.ac.jp/R/kruskal-wallis.html

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

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

==========#

using StatsBase, Rmath, Combinatorics, Printf

function kruskalwallis(data, group)
    idx, ni = table(group)
    n = sum(ni)
    r = tiedrank(data)
    Ri = [sum(r[group .== i]) for i in idx]
    chisq = 12 * sum(Ri .^ 2 ./ ni) / (n * (n + 1)) - 3 * (n + 1)
    if length(unique(data)) != n
        idx2, tie = table(data)
        chisq /= 1 - sum(tie .^ 3 - tie) / (n ^ 3 - n)
    end
    a = length(ni)
    df = a - 1
    pvalue = pchisq(chisq, df, false)
    println("chisq. = $chisq,  df = $df,  p value = $pvalue")
    Dict()
    amean = (n + 1) / 2
    Rmean = Ri ./ ni
    V = sum((r .- amean) .^ 2) / (n - 1)
    pair = combinations(1:a, 2)
    S = [(Rmean[i]- Rmean[j]) ^ 2 / (V / ni[i] + V / ni[j]) for (i, j) in pair]
    P = pchisq.(S, df, false)
    pair2 = [(i, j) for (i, j) in pair]
    @printf("%2s: %2s  %7s  %9s\n", "i", "j", "Chi sq.", "p value")
    for i in 1:a
        @printf("%2d: %2d  %7.3f  %9.7f\n", pair2[i][1], pair2[i][2], S[i], P[i])
    end
    Dict(:chisq => chisq, :df => df, :pvalue => pvalue, :S => S, :P => P)
end

function table(x)
    indices = sort(unique(x))
    counts = zeros(Int, length(indices))
    for i = 1:length(indices)
        counts[i] = count(isequal(indices[i]), x)
    end
    return indices, counts
end

data = [ 3.42, 3.84, 3.96, 3.76,
         3.17, 3.63, 3.47, 3.44, 3.39,
         3.64, 3.72, 3.91 ];

group = [1, 1, 1, 1,  2, 2, 2, 2, 2,  3, 3, 3];

kruskalwallis(data, group)
# chisq. = 5.54871794871795,  df = 2,  p value = 0.06238945711580609
#  i:  j  Chi sq.    p value
#  1:  2    4.104  0.1284601
#  1:  3    0.004  0.9981702
#  2:  3    3.703  0.1570357