Julia でふつうの統計解析

共分散分析 ANCOVA (ANalysis of COVAriance)

前のページに戻る 
#
# GSwR2 第６章 共分散分析 (ANCOVA)
#
# 笠貝の成長速度の季節と密度による変化
#
# 2021/06/02 Daisuke TOMINGAGA

using CSV, Gadfly, DataFrames, DataFramesMeta, Statistics, 
      HypothesisTests, Cairo, Fontconfig, GLM, Distributions,
      Distances, LinearAlgebra

set_default_plot_size(8cm, 6cm)

# データ読み込み、データ形式の確認、要約統計量、プロット
dat = CSV.File("datasets/limpet.csv") |> DataFrame
describe(dat)
p = plot(dat, x = :DENSITY, y = :EGGS, color = :SEASON)
draw(PDF("fig/limpet_1.pdf"), p)

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

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

# 診断プロット２：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))
p = 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"))
draw(PDF("fig/limpet_3.pdf"), p)

# 診断プロット３：Scale-Location
set_default_plot_size(10cm, 10cm) # ちょっと縦長にしないとｙラベルが入らない
sqrtStdRes = sqrt.(abs.(stdRes))
p = plot(x = pred, y = sqrtStdRes, Geom.point,
         Guide.title("Scale-Location"),
         Guide.xlabel("Fitted values"),
         Guide.ylabel("(|Standardized residuals|)1/2"))
draw(PDF("fig/limpet_4.pdf"), p)
set_default_plot_size(8cm, 6cm)

# 診断プロット４：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
p = plot(x = leverage, y = stdRes, Geom.point,
         Guide.xlabel("leverage"), Guide.ylabel("stdRes"))
draw(PDF("fig/limpet_5.pdf"), p)

# モデルの確認
model
coef(model)
stderror(model)

# データとモデルのプロット
# まず横軸の生成、元々のデータでは4点だったところを各季節 10 点にする
xmin = describe(dat)[1, :min]
xmax = describe(dat)[1, :max]
x    = DataFrame([range(xmin, stop = xmax, length = 10)], :auto)
rename!(x, [:DENSITY])
# で季節２水準と生育密度 10 点の全ての組み合わせのデータを作る
v = vec(collect(Base.product(x[!, :DENSITY], levels(dat[!, :SEASON]))))
newx = DataFrame(collect.(collect(zip(v...))), :auto)
rename!(newx, [:DENSITY, :SEASON])
# でモデルの値を計算してデータフレームにする
pred = predict(model, newx)
insertcols!(newx, 1, :PREDICTION => pred)
pred = newx

# プロットを作る
# まず spring と summer、SEASON を日本語にする
for i in 1:length(dat[!, :SEASON])
   dat[i, :SEASON] = (dat[i, :SEASON] == "spring" ? "春" : "夏")
end
for i in 1:length(pred[!, :SEASON])
   pred[i, :SEASON] = (pred[i, :SEASON] == "spring" ? "春" : "夏")
end
rename!(dat,  [:DENSITY, :季節, :EGGS])
rename!(pred, [:PRED, :DENSITY, :季節])
# プロット
p = plot(layer(dat,  x = :DENSITY, y = :EGGS, color = :季節, Geom.point),
         layer(pred, x = :DENSITY, y = :PRED, color = :季節, Geom.line),
         Guide.title("カサガイ産卵数に対する生育密度の影響"),
         Guide.xlabel("生育密度"),
         Guide.ylabel("産卵量"))
draw(PDF("fig/limpet_6.pdf", 10cm, 10cm), p)
2021, © Daisuke TOMINAGA.