Julia でふつうの統計解析
一般化線形モデル GLM (Generalized Linear Model)
前のページに戻る
#
# GSwR2 第７章 一般化線形モデル
#
# 羊の生涯繁殖成功度
#
# 2021/06/02 Daisuke TOMINGAGA
using CSV, Gadfly, DataFrames, DataFramesMeta, Statistics,
HypothesisTests, Cairo, Fontconfig, GLM, Distributions, Distances,
LinearAlgebra
# プロットサイズはあまり小さくすると軸ラベルが収まらなくなり、たとえばｙ軸の
# ラベルは元々縦なのに水平になったりする
set_default_plot_size(10cm, 9cm)
# データ読み込み、データ形式の確認、要約統計量
dat = CSV.File("datasets/SoaySheepFitness.csv") |> DataFrame
describe(dat) # Int と Float の二列、それぞれ子の数と体格
rename!(dat, [:fitness, :bodysize]) # 列名からピリオド取る
# プロット、二つの実数変数だけのデータなので、散布図を見る
p = plot(dat, x = :bodysize, y = :fitness, Geom.point)
draw(PDF("fig/soay_1.pdf"), p)
model = glm(@formula(fitness ~ bodysize), dat, Poisson(), LogLink())
pred = predict(model)
res = dat[!, :fitness] .- pred # 逸脱度残差の符号を決めるのだけに使う
p = plot(layer(dat, x = :bodysize, y = :fitness, Geom.point),
layer( x = dat[!, :bodysize], y = pred, Geom.line))
draw(PDF("fig/soay_7.pdf"), p)
# 診断プロット１：Residuals vs Fitted
# 縦軸はただの残差ではなく、逸脱度残差 deviance residual
# http://juliastats.github.io/GLM.jl/dev/api/#GLM.devresid
# 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)
p = plot(x = logpred, y = devRes, Geom.point,
Guide.title("Residuals vs Fitted"),
Guide.xlabel("Fitted values"),
Guide.ylabel("Residuals"))
draw(PDF("fig/soay_8.pdf"), p)
# 診断プロット２：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))
p = 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"))
draw(PDF("fig/soay_9.pdf"), p)
# 診断プロット３：Scale-Location
sqrtStdRes = sqrt.(abs.(stdDevRes))
p = plot(x = logpred, y = sqrtStdRes, Geom.point,
Guide.title("Scale-Location"),
Guide.xlabel("Fitted values"),
Guide.ylabel("|Standardized residuals|1/2"))
draw(PDF("fig/soay_10.pdf"), p)
# 診断プロット４：Residuals vs Leverage
# てこ比（レバレッジ）の計算には重み行列が必要だが、glm の返すオブジェクトの
# 中にその対角成分ベクトルがある（繰り返し計算が収束した時の重み行列、glmfit.jl
# で定義されている GlmResp 構造体の wrkwt フィールド）
# sum(leverage) は (説明変数の数 + 1) = 2.0 になる
# http://civil.colorado.edu/~balajir/CVEN6833/lectures/
# glm-estimation-presentation.pdf
# https://forum.vsni.co.uk/viewtopic.php?t=1188
# https://uk.sagepub.com/sites/default/files/upm-binaries/
# 17840_Chapter_6.pdf
# https://courses.ms.ut.ee/MTMS.01.011/2018_spring/uploads/Main/
# GLM_slides_6_binary_response.pdf
X = [dat[!, :bodysize] 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) # ハット行列の対角成分がレバレッジ
p = plot(x = leverage, y = stdDevRes, Geom.point,
Guide.title("Residuals vs Leverage"),
Guide.xlabel("Leverage"), Guide.ylabel("Std. Pearson resid."))
draw(PDF("fig/soay_11.pdf"), p)
# ということで、データとモデルと信頼区間のプロット
# まずｘ軸の範囲を決める、データの範囲を探して､50点に区切る
# データフレームに入れて、回帰時の変数名と同じ変数名を付ける
xmin = describe(dat)[2, :min]
xmax = describe(dat)[2, :max]
newx = DataFrame([range(xmin, stop = xmax, length = 50)], :auto)
rename!(newx, [:bodysize])
# 予測値を計算して、結果のデータフレームにｘの値を挿入
# GLM のモデルからは信頼区間は計算してくれないので計算するコードを書く
# 線形予測子の空間では残差は正規分布でその分散は一定という前提なので、
# 線形空間での残差の標準偏差/N を標準誤差とする
# ただ観測値が 0 だと対数空間では -Inf になって線形モデルの空間での信頼区間を
# 計算できないので、その点はプロットから除く
pred = DataFrame(fitness = predict(model, newx)) # 予測値ベクトル
logpred = log.(pred[!, :fitness]) # 線形予測子の値ベクトル
logres = log.(dat[!, :fitness]) .- logpred # 線形空間での残差ベクトル
x = newx[!, :bodysize] # 説明変数値ベクトル
N = length(pred[!, :fitness]) # サンプル数
k = 0 # カウンター
for i in 1:N
global k += 1 # REPL v1.5 以降なら 'global' は削除可
if (dat[i, :fitness] < 1) # 目的変数値が０だったら...
splice!(logres, k) # プロット対象から削除
splice!(logpred, k)
splice!(x, k)
k -= 1
end
end
newx = DataFrame(bodysize = x) # 削除ずみ説明変数値ベクトル
pred = DataFrame(fitness = exp.(logpred)) # 削除ずみ予測値ベクトル
se = std(logres) / sqrt(length(logres)) # 線形空間での標準誤差（一定）
ci = 1.96 .* se # 95% 信頼区間
insertcols!(pred, 1, :bodysize => newx[!, :bodysize]) # データフレームにする
insertcols!(pred, 3, :lower => exp.(logpred .- ci)) # 信頼区間下界
insertcols!(pred, 4, :upper => exp.(logpred .+ ci)) # 信頼区間上界
# データと予測直線をプロット
# 先に置いた Geom.point のレイヤーが Geom.line の上に描画される
# ymin と ymax は Geom.ribbon が使う
p = plot(layer(dat, x = :bodysize, y = :fitness, Geom.point),
layer(pred,x = :bodysize, y = :fitness,
ymin = :lower, ymax = :upper,
Geom.line, Geom.ribbon),
Coord.cartesian(xmin=4.5, xmax=9),
Guide.title("ソアイ羊の生涯繁殖成功度"),
Guide.xlabel("体の大きさ"),
Guide.ylabel("産んだ子羊の数"))
draw(PDF("fig/soay_12.pdf"), p)
2021, © Daisuke TOMINAGA.