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