Julia でふつうの統計解析

ｔ検定 t-test その２

#
# GSwR2 第５章 その２ ｔ検定
#
# 農場のオゾン濃度による作物の生育速度の比較
#
# 2019/03/09 Daisuke TOMINGAGA

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

# データ読み込み、データ形式の確認
dat = CSV.read("datasets/ozone.csv", allowmissing=:none)
# カラム名に . が入ってるとあとでうまく扱えないので､名前を変える
names!(dat, [:Ozone, :Location, :ID])

# 要約統計量、プロット
# データは要するに二群、各群でのオゾン濃度の平均値と分散
by(dat, [:Location], :Ozone => mean)
by(dat, [:Location], :Ozone => var)

p = plot(dat, y = :Ozone, Geom.boxplot)
draw(PDF("ozone_1.pdf", 10cm, 8cm), p)

# グループごとにまとめた棒グラフ
p = plot(dat, x=:Ozone, color=:Location, Geom.histogram())
draw(PDF("ozone_2.pdf", 10cm, 8cm), p)

p = plot(dat, ygroup=:Location, x=:Ozone,
     Geom.subplot_grid(Geom.histogram(bincount = 10)))
draw(PDF("ozone_3.pdf", 10cm, 8cm), p)

# でｔ検定
# ベクトルを二本作らねばならない
datW = @linq dat |> where(:Location .== "West")
datE = @linq dat |> where(:Location .== "East")
EqualVarianceTTest(datW[:, :Ozone], datE[:, :Ozone])
UnequalVarianceTTest(datW[:, :Ozone], datE[:, :Ozone])
pvalue(EqualVarianceTTest(datW[:, :Ozone], datE[:, :Ozone]))
pvalue(UnequalVarianceTTest(datW[:, :Ozone], datE[:, :Ozone]))

p = plot(dat, ygroup=:Location, x=:Ozone,
         Guide.title("農場の場所によるオゾン濃度の分布"),
         Guide.xlabel("オゾン濃度"),
         Guide.ylabel("測定地点数"),
         Geom.subplot_grid(Geom.histogram(bincount = 10)))
draw(PDF("ozone_4.pdf", 10cm, 8cm), p)