Juliaでt検定

Rをはじめよう 生命科学のためのRStudio入門のデータセット（ここにある）を使って､ｔ検定をやるJuliaコード。データを読み込み､図で結論を予想し、検定し､図を仕上げる。コードや出力の説明などは技術書典６(2019/04/14)で頒布予定。

compensation.jl

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

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

# データ点も入れた箱ヒゲ図
# 後に置いた layer から絵が描かれるようなので、Geom.boxplot を先に置くと
# Geom.point が箱に隠れて見えなくなる
draw(PDF("fig/compensation_2.pdf", 10cm, 8cm),
     hstack(plot(dat, layer(x = "Grazing", y = "Fruit", Geom.point),
                      layer(x = "Grazing", y = "Fruit", Geom.boxplot,
                           Theme(default_color=colorant"red"))),
           plot(dat, layer(x = "Grazing", y = "Root", Geom.point),
                     layer(x = "Grazing", y = "Root", Geom.boxplot,
                           Theme(default_color=colorant"red")))))

# ｔ検定のために、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)

# 全体の散布図
draw(PDF("fig/compensation_3.pdf", 10cm, 8cm),
     plot(dat, x="Root", y="Fruit", color="Grazing", Geom.point))

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


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

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

# 散布図その２ https://qiita.com/Nabetani/items/ab8601f02176d63bc9e2
# 、プロットの色を選んで、ラベルも日本語にする
# 対話的に操作している時はラベルに日本語が使えるが、スクリプトを直接読み込ませて
# 実行すると、日本語の部分でエラーになる
draw(PDF("fig/compensation_5.pdf", 10cm, 8cm),
     plot(dat, x = :Root, y = :Fruit, color = :Grazing,
          Guide.title("リンゴの幹の太さと収穫量"),
          Guide.xlabel("台木直径"), Guide.ylabel("収穫量")))