Julia でふつうの統計解析

χ²検定 chi-squared test

#
# GSwR2 第５章　その１ カイ二乗検定
#
# 2019/03/09 Daisuke TOMINGAGA

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

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

# 要約統計量
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)  # 分散

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

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

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

# きれいなプロット
p = 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()))
draw(PDF("ladybird_2.pdf", 10cm, 8cm), p)