JuliaでGLM (一般化線形モデル)

Rをはじめよう 生命科学のためのRStudio入門のデータセット（ここにある）を使って､一般化線形モデル（Generalized Linear Model, GLM）をデータに当てはめるJuliaコード。データを読み込み､図で結論を予想し、モデルの健全性を診断し、検定し､図を仕上げる。コードや出力の説明などは技術書典６(2019/04/14)で頒布予定。

soay.jl

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

# データ読み込み、データ形式の確認、要約統計量
dat = CSV.read("datasets/SoaySheepFitness.csv", allowmissing=:none)
names!(dat, [:fitness, :body_size])
describe(dat)

# プロット、二つの実数変数だけのデータなので、散布図を見る
draw(PDF("fig/soay_1.pdf", 10cm, 8cm),
     plot(dat, x = :body_size, y = :fitness, Geom.point))

# モデル当てはめと、かるく様子見
model = glm(@formula(fitness ~ body_size), dat, Poisson(), LogLink())
pred = predict(model)
res = dat[:fitness] .- pred # 逸脱度残差の符号を決めるのだけに使う
draw(PDF("fig/soay_7.pdf", 10cm, 8cm),
     plot(layer(dat, x = :body_size, y = :fitness, Geom.point),
          layer(x = dat[:body_size], y = pred, Geom.line)))

# 診断プロット１：Residuals vs Fitted
# 縦軸はただの残差ではなく、逸脱度残差 deviance residual
# devresid = 2 * (y * log(y / μ) - (y - μ)) (glmtools.jl)
logpred = log.(pred)
devRes = sqrt.(devresid.(Poisson.(pred), dat[:fitness], pred)) .* sign.(res)
# または sqrt.(model.model.rr.devresid) .* sign.(res)
draw(PDF("fig/soay_8.pdf", 10cm, 8cm),
     plot(x = logpred, y = devRes, Geom.point,
          Guide.title("Residuals vs Fitted"),
          Guide.xlabel("Fitted values"),
          Guide.ylabel("Residuals")))

# 診断プロット２：Normal Q-Q plot
# 残差は標準化逸脱残差 standardized deviance residual
# 逸脱度は -log(p), p はモデル曲線による期待値からデータ点の生じる確率
# 逸脱度残差は、そのデータ点の逸脱度
N = length(devRes)
stdDevRes = copy(devRes)
stdDevRes = (devRes .- ones(N) .* mean(devRes)) ./ std(devRes)
qx = quantile.(Normal(), range(0.5, stop=(N-0.5), length=(N)) ./ (N + 1))
draw(PDF("fig/soay_9.pdf", 10cm, 8cm),
     plot(layer(x = qx, y = sort(stdDevRes), 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
# GSwR2 とちょっと違うところがある
sqrtStdRes = sqrt.(abs.(stdDevRes))
draw(PDF("fig/soay_10.pdf", 10cm, 8cm),
     plot(x = logpred, y = sqrtStdRes, Geom.point,
          Guide.title("Scale-Location"),
          Guide.xlabel("Fitted values"),
          Guide.ylabel("|Standardized residuals|<sup>1/2</sup>")))

# 診断プロット４：Residuals vs Leverage
# てこ比（レバレッジ）の計算には重み行列が必要だが、glm の返すオブジェクトの
# 中にその対角成分ベクトルがある（繰り返し計算が収束した時の重み行列、glmfit.jl
# で定義されている GlmResp 構造体の wrkwt フィールド）
# sum(leverage) は (説明変数の数 + 1) = 2.0 になる
X = [dat[:body_size] ones(N)]          # デザイン行列、変数値の列と切片用の列
W = diagm(0 => model.model.rr.wrkwt)   # 対角成分が working weight の対角行列
W2 = real.(sqrt(W))                    # W2 * W2 = W、虚部は理論的には０になる
H = W2 * X * inv(X' * W * X) * X' * W2 # ハット行列
leverage = diag(H)                     # ハット行列の対角成分がレバレッジ
draw(PDF("fig/soay_11.pdf", 10cm, 8cm),
     plot(x = leverage, y = stdDevRes, Geom.point,
          Guide.title("Residuals vs Leverage"),
          Guide.xlabel("Leverage"), Guide.ylabel("Std. Pearson resid.")))

# ということで、データとモデルと信頼区間のプロット
# まずｘ軸の範囲を決める、データの範囲を探して､50点に区切る
# データフレームに入れて、回帰時の変数名と同じ変数名を付ける
xmin = describe(dat)[2, :min]
xmax = describe(dat)[2, :max]
newx = DataFrame([range(xmin, stop = xmax, length = 50)])
names!(newx, [:body_size])
# 予測値を計算して、結果のデータフレームにｘの値を挿入
# GLM のモデルからは信頼区間は計算してくれないので計算するコードを書く
# 線形予測子の空間では残差は正規分布でその分散は一定という前提なので、
# 線形空間での残差の標準偏差/N を標準誤差とする
# ただ観測値が 0 だと対数空間では -Inf になって線形モデルの空間での信頼区間を
# 計算できないので、その点はプロットから除く
pred = DataFrame(fitness = predict(model, newx)) # 予測値ベクトル
logpred = log.(pred[:fitness])                   # 線形予測子の値ベクトル
logres = log.(dat[:fitness]) .- logpred          # 線形空間での残差ベクトル
x = newx[:body_size]                             # 説明変数値ベクトル
N = length(pred[:fitness])                       # サンプル数
k = 0
for i in 1:N
  global k += 1
  if (dat[i, :fitness] < 1)                      # 目的変数値が０だったら
    splice!(logres, k)                           # プロット対象から削除
    splice!(logpred, k)
    splice!(x, k)
    k -= 1
  end
end
newx = DataFrame(body_size = x)                  # 削除ずみ説明変数値ベクトル
pred = DataFrame(fitness = exp.(logpred))        # 削除ずみ予測値ベクトル
se = std(logres) / sqrt(length(logres))          # 線形空間での標準誤差（一定）
ci = 1.96 .* se                                  # 95% 信頼区間
insertcols!(pred, 1, body_size=newx[:body_size]) # データフレームにまとめる
insertcols!(pred, 3, lower = exp.(logpred .- ci))# 観測値空間での信頼区間下界
insertcols!(pred, 4, upper = exp.(logpred .+ ci))# 観測値空間での信頼区間上界
# データと予測曲線をプロット
# 先に置いた Geom.point のレイヤーが Geom.line の上に描画される
# ymin と ymax は Geom.ribbon が使う
draw(PDF("fig/soay_12.pdf", 10cm, 8cm),
     plot(layer(dat, x = :body_size, y = :fitness, Geom.point),
          layer(pred,x = :body_size, y = :fitness,
                     ymin = :lower, ymax = :upper,
                     Geom.line, Geom.ribbon),
          Coord.cartesian(xmin=4.5, xmax=9),
          Guide.title("ソアイ羊の生涯繁殖成功度"),
          Guide.xlabel("体の大きさ"),
          Guide.ylabel("産んだ子羊の数")))