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