Julia でふつうの統計解析

χ²検定 chi-squared test

#
# GSwR2 第５章　その１ カイ二乗検定
#
# 2021/06/02 Daisuke TOMINGAGA

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

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

# 要約統計量
describe(dat)
combine(groupby(dat, :Habitat), :number => mean) # 各地域での個体数の平均値

# 列に m という名前を付けてみる
combine(groupby(dat, :Habitat), :number => mean => :m)

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

plot(dat, x = :morph_colour, y = :number)
plot(dat, x = :morph_colour, y = :number, Geom.bar)
plot(dat, x = :morph_colour, y = :number, Geom.boxplot)
plot(dat, x = :Habitat,  y = :number, color = :morph_colour, Geom.bar)
plot(dat, x = :Habitat,  y = :number, color = :morph_colour, 
     Geom.boxplot)

# グループごとにまとめた棒グラフを作って、PDF ファイルに保存する
p = plot(dat, xgroup = :Habitat, x = :morph_colour, y = :number,
         color = :morph_colour, Geom.subplot_grid(Geom.bar()))
draw(PDF("fig/ladybird_7.pdf", 8cm, 6cm), p)


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

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

# きれいなプロット
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("fig/ladybird_8.pdf", 12cm, 8cm), p)

# フィッシャー検定
FisherExactTest(30,70,115,85)