Juliaで共分散分析

[Rをはじめよう 生命科学のためのRStudio入門](https://www.yodosha.co.jp/yodobook/book/9784758120951/index.html)のデータセット（[ここ](http://r4all.org/books/datasets/)にある）を使って､共分散分析（ANalysis of COVAriance, ANCOVA）をやるJuliaコード。データを読み込み､図で結論を予想し、検定し､図を仕上げる。コードや出力の説明などは[技術書典６](https://techbookfest.org/event/tbf06)(2019/04/14)で頒布予定。

```limpet.jl
using CSV, Gadfly, DataFrames, DataFramesMeta, Statistics, HypothesisTests,
      Cairo, Fontconfig, GLM, Distributions, Distances, LinearAlgebra

# データ読み込み、データ形式の確認、要約統計量、プロット
dat = CSV.read("datasets/limpet.csv", allowmissing=:none)
describe(dat)
describe(dat[:SEASON])
draw(PDF("fig/limpet_1.pdf", 10cm, 8cm),
     plot(dat, x = :DENSITY, y = :EGGS, color = :SEASON))

# 線形モデル当てはめ
model = lm(@formula(EGGS ~ DENSITY * SEASON), dat)

# 診断プロット１：Residuals vs Fitted
pred = predict(model)
res = dat[:EGGS] - pred
draw(PDF("fig/limpet_2.pdf", 10cm, 8cm),
     plot(x = pred, y = res, Geom.point,
          Guide.title("Residuals vs Fitted"),
          Guide.xlabel("Fitted values"), Guide.ylabel("Residuals")))

# 診断プロット２：Normal Q-Q plot
m = mean(res)
N = length(res)
stdRes = (res - m * ones(N)) / (std(res)) # 標準化誤差 = (誤差-平均)/標準偏差
qx = quantile.(Normal(), range(0.5, stop=(N-0.5), length=(N)) / (N + 1))
draw(PDF("fig/limpet_3.pdf", 10cm, 8cm),
     plot(layer(x = qx, y = sort(stdRes), Geom.point),
          layer(x = [-3,3], y = [-3,3], Geom.line, style(line_style = [:dot])),
          Guide.title("Normal Q-Q plot"),
          Guide.xlabel("Theoretical Quantiles"),
          Guide.ylabel("Standardized residuals")))

# 診断プロット３：Scale-Location
sqrtStdRes = sqrt.(abs.(stdRes))
draw(PDF("fig/limpet_4.pdf", 10cm, 10cm),
     plot(x = pred, y = sqrtStdRes, Geom.point,
         Guide.title("Scale-Location"),
         Guide.xlabel("Fitted values"),
         Guide.ylabel("(|Standardized residuals|)<sup>1/2</sup>")))

# 診断プロット４：Residuals vs Leverage
x = dat[:DENSITY]
Sxx = sum(x .* x) - sum(x)^2/N
m = mean(x)
leverage = ones(length(res))
for i in 1:N leverage[i] = (x[i] - m)^2 / Sxx + 1/N end
draw(PDF("fig/limpet_5.pdf", 10cm, 10cm),
     plot(x = leverage, y = stdRes, Geom.point,
          Guide.xlabel("leverage"), Guide.ylabel("stdRes")))

model            # もう一度係数を見る
#coef(model)     # 係数だけ見たいならこれでもよい
#stderror(model) # 標準誤差ならこれで見られる

# で、データとモデルと信頼区間のプロット
# まず DENSITY の予測範囲を決め、予測点を生成し、元データと同じ変数名にし、
# DENSITY の値と SEASON の値の全ての組み合わせを生成する -> newx
xmin = describe(dat)[1, :min]
xmax = describe(dat)[1, :max]
x    = DataFrame([range(xmin, stop = xmax, length = 10)])
names!(x, [:DENSITY])
v = vec(collect(Base.product(x[:DENSITY], levels(dat[:SEASON]))))
newx = DataFrame(collect.(collect(zip(v...))))
names!(newx, [:DENSITY, :SEASON])

# でその予測点での目的変数の期待値95％信頼区間の上下の値を計算して一つの
# データフレームにまとめる
pred = predict(model, newx, interval = :confidence)
names!(pred, [:P, :L, :U])
insertcols!(pred, 1, DENSITY = newx[:DENSITY])
insertcols!(pred, 1, SEASON  = newx[:SEASON])

# で仕上げにデータと予測直線をプロット
# 先に置いた Geom.point が Geom.line の上に描画される
# ymin と ymax は Geom.ribbon が使う
# プロットをできるだけ日本語化するために、SEASON変数の中身を日本語にする
dat2 = deepcopy(dat)
pred2 = deepcopy(pred)
for i in 1:length(dat2[:SEASON])
   dat2[i, :SEASON] = ( dat2[i, :SEASON] == "spring" ? "春" : "夏")
end
# 上の for ループと同じことを pred2 に対して違うやり方でやる
pred2[pred2[:SEASON] .== "spring", :SEASON] = "春"
pred2[pred2[:SEASON] .== "summer", :SEASON] = "夏"
names!(dat2,  [:DENSITY, :季節, :EGGS])
names!(pred2, [:季節, :DENSITY, :P, :L, :U])
draw(PDF("fig/limpet_6.pdf", 10cm, 10cm),
     plot(layer(dat2,  x = :DENSITY, y = :EGGS, color = :季節, Geom.point),
          layer(pred2, x = :DENSITY, y = :P, Geom.line,
                       ymin = :L, ymax = :U, Geom.ribbon,
                       color = :季節),
         Guide.title("カサガイ産卵数に対する生育密度の影響"),
         Guide.xlabel("生育密度"),
         Guide.ylabel("産卵量")))
```