Julia でふつうの統計解析

分散分析 ANOVA (ANalysis Of Variance)

前のページに戻る 
#
# 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)

# データ読み込み
dat = CSV.read("datasets/Daphniagrowth.csv", allowmissing=:none)
# カラム名に . が入ってるとあとでうまく扱えないので､名前を変える
names!(dat, [:parasite, :rep, :growth])
describe(dat)
describe(parasite)
std(dat[:rep])
std(dat[:growth])

# 全体、および寄生虫ごとの箱ヒゲ図
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))

# 分散分析
# ANOVA を一発でやってくれる関数はないので、必要な計算をそのままやる
# そのためにデータフレームを、解析しやすい用に unstack() で整形する
dat2 = unstack(dat, :parasite, :growth)         # 寄生虫ごとに列になるように
deletecols!(dat2, :rep)                         # 不要な列を消し列名を変える
names!(dat2, [:Metschnikowia_bicuspidata,       # 寄生虫１ メチニコビア酵母
              :Pansporella_perplexa,            # 寄生虫２ アメーバの一種
              :Pasteuria_ramosa,                # 寄生虫３ パスツリア菌
              :Control])                        # 寄生虫無し
# Control が最初の列になるように並べ替え
permutecols!(dat2, [:Control, :Metschnikowia_bicuspidata, :Pansporella_perplexa,
                    :Pasteuria_ramosa])         # アルファベット順に並べる
dat3 = dat2[:, :]                               # あとで使うので変更前に複製
nG = length(levels(dat[:, :parasite]))          # 群の数
nR = length(dat2[:, :Control])                  # 各群のサンプル数
means = by(dat, :parasite, M = :growth => mean) # 群ごとの平均
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]
  end
end
dDoF = (nG * (nR - 1))                          # 群内の自由度
denominator = (sum(DataFrames.colwise(sum, dat2))) / dDoF # 分母
m = mean(DataFrames.colwise(mean, dat3))        # 全サンプルの平均
for i in 1:nG
  means[i] = means[i] - m                       # 各群の平均から全体の平均を
  means[i] = means[i] * means[i]                # 引いて二乗する
end
nDoF = (nG - 1)                                 # 群間の自由度
numerator = sum(means) * nR / nDoF              # 分子
ccdf(FDist(nDoF, dDoF), numerator/denominator)  # ｐ値

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