Julia でふつうの統計解析

ｔ検定 t-test その１

#
# GSwR2 第４章 ｔ検定
# compensations.csv を読み込んで、プロットして、Grazed/Ungrazed でｔ検定する
#
# スクリプト全体では約70秒
#
# 2019/03/10 Daisuke TOMINAGA

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

set_default_plot_size(10cm, 8cm)

# データの読み込みと形式の確認、要約統計量の表示、箱ヒゲ図のプロット
dat = CSV.read("datasets/compensation.csv", allowmissing=:none)
describe(dat)
p = hstack(plot(dat, x = :Grazing, y = :Fruit, Geom.boxplot),
           plot(dat, x = :Grazing, y = :Root,  Geom.boxplot))
draw(PDF("compensation_1.pdf", 10cm, 8cm), p)

# ｔ検定のために、Grazing の値によって dat を二つに分ける
datG = @linq dat |> where(:Grazing .== "Grazed") |>
                    select(:Root, :Fruit, :Grazing)
datU = @linq dat |> where(:Grazing .== "Ungrazed") |>
                    select(:Root, :Fruit, :Grazing)
describe(datG)
describe(datU)

# 全体の散布図 https://nzw0301.github.io/2015/12/gadfly
p = plot(dat, x="Root", y="Fruit", color="Grazing", Geom.point)
draw(PDF("compensation_2.pdf", 10cm, 8cm), p)

# 要約統計量
# Grazing の内容ごとに群に分け、各変数で平均と分散を計算する
# ３群ｘ２変数ｘ平均と分散が表示される
by(dat, :Grazing, :Root => mean, :Root => var, :Fruit => mean, :Fruit => var)

# 箱ヒゲ図 http://gadflyjl.org/v0.6.3/lib/geoms/geom_boxplot.html
# 後に置いた layer から絵が描かれるようなので、Geom.boxplot を先に置くと、
# Geom.point が箱に隠れて見えなくなる
p = plot(dat, layer(x = "Grazing", y = "Fruit", Geom.point),
              layer(x = "Grazing", y = "Fruit", Geom.boxplot,
                    Theme(default_color=colorant"red")), Guide.ylabel("収穫量"))
#                   Theme(default_color=colorant"red")), Guide.ylabel("Fruit"))
draw(PDF("compensation_3.pdf", 10cm, 8cm), p)
p = plot(dat, layer(x = "Grazing", y = "Root", Geom.point),
              layer(x = "Grazing", y = "Root", Geom.boxplot,
                  Theme(default_color=colorant"red")), Guide.ylabel("台木直径"))
#                 Theme(default_color=colorant"red")), Guide.ylabel("Root"))
draw(PDF("compensation_4.pdf", 10cm, 8cm), p)

# パラメトリック検定をやっていいかどうか、ヒストグラムで見てみる
p = plot(dat, ygroup = :Grazing, x = :Fruit,
     Geom.subplot_grid(Geom.histogram(bincount = 15)))
draw(PDF("compensation_4.pdf", 10cm, 8cm), p)

# Grazed と Ungrazed の平均値の差を検定、等分散と非等分散、収穫量と台木直径
EqualVarianceTTest(datG[:, :Fruit], datU[:, :Fruit])
EqualVarianceTTest(datG[:, :Root], datU[:, :Root])
UnequalVarianceTTest(datG[:, :Fruit], datU[:, :Fruit])
UnequalVarianceTTest(datG[:, :Root], datU[:, :Root])
#pvalue(EqualVarianceTTest(datG[:, :Fruit], datU[:, :Fruit]))
#pvalue(EqualVarianceTTest(datG[:, :Root], datU[:, :Root]))
#pvalue(UnequalVarianceTTest(datG[:, :Fruit], datU[:, :Fruit]))
#pvalue(UnequalVarianceTTest(datG[:, :Root], datU[:, :Root]))

# 散布図その２ https://qiita.com/Nabetani/items/ab8601f02176d63bc9e2
# 対話的に操作している時はラベルに日本語が使えるが、スクリプトを直接読み込ませて
# 実行すると、日本語の部分でエラーになる
p = plot(layer(datG, x = "Root", y = "Fruit",
               Geom.point, Theme(default_color=colorant"red")),
         layer(datU, x = "Root", y = "Fruit", Geom.point),
         Guide.title("リンゴの幹の太さと収穫量"),
         Guide.xlabel("台木直径"), Guide.ylabel("収穫量"))
#        Guide.title("Apple tree size vs Harvest"),
#        Guide.xlabel("Root"), Guide.ylabel("Fruit"))
draw(PDF("compensation_4.pdf", 10cm, 8cm), p)