「Julia でANOVA」 という記事がある。記事の日付は2019年3月15日とわずか2年足らず前のものである。Julia v1.1.1 がリリースされたのが 2019/05/16 であるから，少なくともそれより前のバージョンということになる。
しかしその後の Julia の発展が著しく，今では廃止された関数があったりして，そのままでは参考にならない。また，示された一元配置分散分析は各群の例数が等しく，また，各群の等分散を仮定する検定なので，等分散を仮定しない検定も示しておく。
廃止あるいは仕様変更のあった関数
dat = CSV.read("foo.csv", allowmissing=:none)
は
dat = CSV.read("lfoo.csv", DataFrame)
と書くようになった。今でも多くのページで単にファイル名を指定すればよいと書いてあるので，小生も最初は困り果てた。
また，allowmissing は不要になった（当然のものとして組み込まれた）。
by() によるグループごとの統計処理は廃止された。
by(dat, :parasite, :growth => mean)
は groupby() と combine() で以下のように書く。 combine() は非常に強力である。
gd = groupby(dat, :parasite);
combine(gd, :growth => mean)
stat = combine(gd, nrow => :sample_size, :growth_rate => mean => :group_mean, :growth_rate => var => :group_var)
名前の変更は names!() ではなく rename!()，
列の削除は deletecols!() ではなく，select!()
列の並べ替えも permutecols!() ではなく，select!() で行える。
さて，以上の点に留意して，以下のような処理手順を示す。
まずは，データ入力から前処理まで。
using CSV, DataFrames, Statistics, Distributions
dat = CSV.read("Daphniagrowth.csv", DataFrame);
# '.' を使っている列名は，最初に '_' に変更しておく
rename!(dat, "growth.rate" => "growth_rate");
# 群ごとの統計は groupby()， combine() を使う
gd = groupby(dat, :parasite);
stat = combine(gd, nrow => :sample_size, :growth_rate => mean => :group_mean, :growth_rate => var => :group_var)
println(stat)
4×4 DataFrame
Row │ parasite sample_size group_mean group_var
│ String Int64 Float64 Float64
─────┼───────────────────────────────────────────────────────────────
1 │ control 10 1.21391 0.0130376
2 │ Metschnikowia bicuspidata 10 0.801154 0.0465879
3 │ Pansporella perplexa 10 1.07636 0.0322315
4 │ Pasteuria ramosa 10 0.482203 0.0375737
後は，計算処理を黙々と行う（ドット演算子に注意）
まずは，各群の分散が等しいことを仮定する検定。
k = length(gd)
sample_size = stat[!, :sample_size];
group_mean = stat[!, :group_mean ];
group_var = stat[!, :group_var ];
n = sum(sample_size)
grand_mean = sum(sample_size .* group_mean) / n
St = var(dat[!, :growth_rate]) * (n - 1)
Sb = sum(sample_size .* (group_mean .- grand_mean) .^ 2)
Sw = St - Sb
dfb = k - 1
dfw = n - k
dft = n - 1
Vb = Sb / dfb
Vw = Sw / dfw
# 最終的に必要なのは，F 統計量と，P 値
F = Vb / Vw
p_value = ccdf(FDist(dfb, dfw), F)
# 見やすくするためにデータフレームにしておこう
ANOVA_table = DataFrame(SS=[Sb, Sw, St], df=[dfb, dfw, dft],
MS=[Vb, Vw, ""], F=[F, "", ""], p_value=[p_value, "", ""])
println(ANOVA_table)
3×5 DataFrame
Row │ SS df MS F p_value
│ Float64 Int64 Any Any Any
─────┼─────────────────────────────────────────────────
1 │ 3.13791 3 1.04597 32.3252 2.57128e-10
2 │ 1.16488 36 0.0323577
3 │ 4.30278 39
次は，各群の分散が等しいとは仮定しない，いわゆるウェルチの方法によるものである。
grand_mean = mean(dat[!, :growth_rate])
k = length(gd)
sample_size = stat[!, :sample_size];
group_mean = stat[!, :group_mean];
group_var = stat[!, :group_var];
wi = sample_size ./ group_var;
w = sum(wi)
m = sum(wi .* group_mean) / w
a = sum(((1 .- wi ./ w) .^ 2) ./ (sample_size .- 1)) / (k^2 - 1)
F = sum(wi .* (group_mean .- m) .^ 2) / (k - 1) / (1 + 2a*(k - 2))
dfb = k - 1
dfw = 1 / 3a
p_value = ccdf(FDist(dfb, dfw), F)
println("F = $F, p_value = $p_value")
