Julia でデータフレームのデータで一元配置分散分析，Oneway ANOVA

Julia でふつうの統計解析　分散分析 ANOVA (ANalysis Of Variance)
https://staff.aist.go.jp/tominaga-daisuke/julia/daphnia.html

で，かなり面倒なことをやっている。

しかも，各群の例数が異なることに対応していない。

また，そんなに古い記事でもないのに，Obsolete になってしまった関数があったり。

別法は，最後の方に記述。
#
# GSwR2 第５章 その４ ANOVA
#
# ミジンコの寄生虫と生育速度のデータ
# スクリプト全体ではで測ると３２秒くらい
#
# 2019/03/10 Daisuke TOMINGAGA

using CSV, Gadfly, DataFrames, DataFramesMeta, Statistics, HypothesisTests,
      Cairo, Fontconfig, GLM, Distributions, Distances, LinearAlgebra

set_default_plot_size(10cm, 8cm)

# データ読み込み
# NO! dat = CSV.read("datasets/Daphniagrowth.csv", allowmissing=:none)
dat = CSV.read("datasets/Daphniagrowth.csv", DataFrame)

# カラム名に . が入ってるとあとでうまく扱えないので､名前を変える
# NO! names!(dat, [:parasite, :rep, :growth])
rename!(dat, [:parasite, :rep, :growth])

describe(dat)
# NO! describe(parasite)
describe(dat[!, :parasite]) # describe(dat.parasite) is more better
# NO! std(dat[:rep])
std(dat[!, :rep]) # std(dat.rep) is more better
# NO! std(dat[:growth])
std(dat[!, :growth]) # std(dat.growth) is more better

# 全体、および寄生虫ごとの箱ヒゲ図
p = plot(dat, y = :growth, Geom.boxplot)
draw(PDF("daphnia_1.pdf", 8cm, 18cm), p)

p = plot(dat, x= :parasite, y = :growth, Geom.boxplot)
draw(PDF("daphnia_2.pdf", 8cm, 18cm), p)

# これでも同じようにプロットできる
# plot(dat, xgroup = :parasite, y = :growth, Geom.subplot_grid(Geom.boxplot))
# draw(PDF("daphnia_2-2.pdf", 8cm, 18cm), p)
# y 軸ラベルの位置が変

# 分散分析
# ANOVA を一発でやってくれる関数はないので、必要な計算をそのままやる
# そのためにデータフレームを、解析しやすい用に unstack() で整形する
# WARNIG！！ 群ごとの例数が同じだから以下でよいが，そうでない場合のことも考慮して，別法を採るべし。また， unstack() よりは groupby() のほうがよい。

dat2 = unstack(dat, :parasite, :growth);         # 寄生虫ごとに列になるように
first(dat2)
# Obsolete！ deletecols!(dat2, :rep)             # 不要な列を消し列名を変える
select!(dat2, Not(:rep))

# Obsolete！ names!(dat2, [:Metschnikowia_bicuspidata, # 寄生虫１ メチニコビア酵母
#              :Pansporella_perplexa,            # 寄生虫２ アメーバの一種
#              :Pasteuria_ramosa,                # 寄生虫３ パスツリア菌
#              :Control])                        # 寄生虫無し
# Control が最初の列になるように並べ替え
# Obsolete! permutecols!(dat2, [:Control, :Metschnikowia_bicuspidata, :Pansporella_perplexa,
#                    :Pasteuria_ramosa])        # アルファベット順に並べる
rename!(dat2, "control" => :Control,
      "Metschnikowia bicuspidata" => :Metschnikowi_bicuspidata,
      "Pansporella perplexa" => :Pansporella_perplexa,
      "Pasteuria ramosa" => :Pasteuria_ramosa)

# dat3 = dat2[:, :]                               # あとで使うので変更前に複製
dat3 = copy(dat2)

nG = length(levels(dat[:, :parasite]))          # 群の数
nR = length(dat2[:, :Control])                  # 各群のサンプル数
# NO! means = by(dat, :parasite, M = :growth => mean) # 群ごとの平均
means = combine(groupby(dat, :parasite), :growth => mean => :M) # 群ごとの平均
# means = means[:M]                               # いるのは平均値の列だけ
means = means.M                               # いるのは平均値の列だけ

for i in 1:nR                                   # 各群のサンプルからその群の
  for j in 1:nG                                 # 平均を引いて二乗する
    # dat2[i, j] = dat2[i, j] - means[j]
    # dat2[i, j] = dat2[i, j] * dat2[i, j]
    dat2[i, j] = (dat2[i, j] - means[j])^2
  end
end

dDoF = (nG * (nR - 1))                          # 群内の自由度
# Obsolete! denominator = (sum(DataFrames.colwise(sum, dat2))) / dDoF # 分母
denominator = sum(mapcols(sum, dat2)[1,:]) / dDoF # 分母
# Obsolete! m = mean(DataFrames.colwise(mean, dat3))        # 全サンプルの平均
# m = mean(mapcols(mean, dat3)[1,:])        # 全サンプルの平均
m = mean(dat.growth)
# for i in 1:nG
#   means[i] = means[i] - m                       # 各群の平均から全体の平均を
#   means[i] = means[i] * means[i]                # 引いて二乗する
# end
means = [(x - m)^2 for x in means]

nDoF = nG - 1                                 # 群間の自由度
numerator = sum(means) * nR / nDoF              # 分子
F = numerator/denominator
ccdf(FDist(nDoF, dDoF), F)  # ｐ値

# きれいなプロットを作る
# これを時計回りに 90° 回転させれば見やすい
# 回転後の横軸の数値は白い箱で隠してテキストボックスを置くなどして描き直す
# means = by(dat, :parasite, M = :growth => mean) # 群ごとの平均を求め直す
means = combine(groupby(dat, :parasite), :growth => mean => :M) # 群ごとの平均

p = plot(layer(means, x = :parasite, y = :M, Geom.point,
               shape = [Shape.diamond],
               Theme(default_color = colorant"red", point_size = 2mm)),
         layer(dat,   x = :parasite, y = :growth, Geom.point),
         layer(dat,   x = :parasite, y = :growth, Geom.boxplot),
         Theme(key_position = :none),
         Guide.ylabel("growth rate"),
         Guide.xticks(orientation=:vertical))
draw(PDF("daphnia_3.pdf", 8cm, 18cm), p)

# 最初から横向きにできれば見やすいが、今の Geom.boxplot は横向きに描いてくれない
set_default_plot_size(18cm, 8cm)
p = plot(layer(means, x = :M, y = :parasite, Geom.point,
               shape = [Shape.diamond],
               Theme(default_color = colorant"red", point_size = 2mm)),
         layer(dat, x = :growth, y = :parasite, Geom.point),
         layer(dat, x = :growth, y = :parasite, Geom.boxplot),
         Theme(key_position = :none),
         Guide.xlabel("growth rate"))
draw(PDF("daphnia_4.pdf", 18cm, 10cm), p)

#############

using CSV, DataFrames, Statistics, Rmath
dat = CSV.read("datasets/Daphniagrowth.csv", DataFrame);
rename!(dat, "growth.rate" => :growth);

n = size(dat, 1)
St = var(dat.growth) * (length(dat.growth) - 1)
grand_mean = mean(dat.growth)

gdat = groupby(dat, :parasite);
k = size(gdat, 1)

Sb = sum(length(gdat[i].growth) * (mean(gdat[i].growth) - grand_mean)^2 for i = 1:k)
Sw = St - Sb
dfb = k - 1
dfw = n - k
MSb = Sb / dfb
MSw = Sw / dfw
F = MSb / MSw
p = pf(F, dfb, dfw, false)
println("F = $F, df1 = $dfb, df2 = $dfw, p value = $p")