Julia でふつうの統計解析

ｔ検定 t-test その１

#
# GSwR2 第４章 ｔ検定
# compensations.csv を読み込んで、プロットして、Grazed/Ungrazed でｔ検定する
#
# 2021/06/02 Daisuke TOMINAGA

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

set_default_plot_size(8cm, 8cm)

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

# 散布図を箱ヒゲ図に重ねる
p = 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"))))
draw(PDF("fig/compensation_2.pdf"), 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)

# 全体の散布図
p = plot(dat, x="Root", y="Fruit", Geom.point)
draw(PDF("fig/compensation_31.pdf"), p)

# Grazing で色分けした散布図 https://nzw0301.github.io/2015/12/gadfly
p = plot(dat, x="Root", y="Fruit", color="Grazing", Geom.point)
draw(PDF("fig/compensation_32.pdf"), p)

# layer を使って重ねるやり方
p = plot(layer(datG, x = "Root", y = "Fruit", Geom.point,
               Theme(default_color=colorant"red")),
         layer(datU, x = "Root", y = "Fruit", Geom.point))
draw(PDF("fig/compensation_33.pdf"), p)

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

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

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

# 仕上げの散布図
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("収穫量"))
draw(PDF("fig/compensation_7.pdf"), p)