Rをはじめよう 生命科学のためのRStudio入門のデータセット(ここにある)を使って、カイ二乗検定をやるJuliaコード。データを読み込み、図で結論を予想し、検定し、図を仕上げる。コードの説明などは技術書典6(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 で4通りの組み合わせに
# ついててんとう虫の個体数を合計し、合計値のカラム名は lady にする
# ここでは表(行列の形式)にはせず、4次元のベクトルのまま
cont = by(dat, [:Habitat, :morph_colour], lady = :number => sum)
# カイ二乗検定:ベクトルのまま渡すと違う計算が行われてp値も変わるので
# 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())))