Juliaでカイ二乗検定

Rをはじめよう 生命科学のためのRStudio入門のデータセット（ここにある）を使って､カイ二乗検定をやるJuliaコード。データを読み込み､図で結論を予想し、検定し､図を仕上げる。コードの説明などは技術書典６(2019/04/14)で頒布予定。

ladybird.jl

using CSV, Gadfly, DataFrames, DataFramesMeta, Statistics, HypothesisTests, Cairo, Fontconfig

# データを読み込んで、表の形式を確認
dat = CSV.read("datasets/ladybirds_morph_colour.csv")

describe(dat) # 要約統計量
by(dat, :Habitat, :number => mean) # 平均値、地域ごと
by(dat, [:Habitat, :morph_colour], lady = :number => sum) # 個体数
by(dat, [:Habitat, :morph_colour], lady = :number => mean) # 平均値、地域、色別
by(dat, [:Habitat, :morph_colour], lady = :number => var) # 分散

# データの様子を図で見てみる
draw(PDF("fig/ladybird_1.pdf", 10cm, 8cm),
     plot(dat, y = :number, Geom.bar))
draw(PDF("fig/ladybird_2.pdf", 10cm, 8cm),
     plot(dat, x = :morph_colour, y = :number, Geom.bar))

# グループごとにまとめた棒グラフ
draw(PDF("fig/ladybird_3.pdf", 10cm, 8cm),
     plot(dat, xgroup = :Habitat, x = :morph_colour, y = :number,
               color = :morph_colour,
               Scale.y_continuous(minvalue = 0),
               Geom.subplot_grid(Geom.bar())))

# 2x2分割表：Habitat x morph_colour = 2x2 で４通りの組み合わせに
# ついててんとう虫の個体数を合計し、合計値のカラム名は lady にする
# ここでは表（行列の形式）にはせず、４次元のベクトルのまま
cont = by(dat, [:Habitat, :morph_colour], lady = :number => sum)

# カイ二乗検定：ベクトルのまま渡すと違う計算が行われてｐ値も変わるので
# reshape() で行列の形にして ChisqTest() に渡す
res = ChisqTest(reshape(cont[:, :lady], 2, 2))

# きれいでわかりやすいプロットを作る
# きれいでわかりやすいプロットを作る
draw(PDF("fig/ladybird_4.pdf", 10cm, 9cm),
     plot(dat, xgroup = :Habitat, x = :morph_colour, y = :number,
               color = :morph_colour,
               Scale.color_discrete_manual("black", "red"),
               Scale.y_continuous(minvalue = 0),
               Guide.ylabel("観測された個体数"),
               Guide.title("住むところ別の各色の天道虫の数"),
               Geom.subplot_grid(Geom.bar())))