Julia でふつうの統計解析
一般化線形モデル GLM (Generalized Linear Model)
前のページに戻る
#
# GSwR2 第７章 一般化線形モデル
#
# 羊の生涯繁殖成功度
#
# 2019/03/10 Daisuke TOMINGAGA
using CSV, Gadfly, DataFrames, DataFramesMeta, Statistics, HypothesisTests,
Fontconfig, GLM, Distributions, Distances, LinearAlgebra, Cairo
# プロットサイズはあまり小さくすると軸ラベルが収まらなくなり、たとえばｙ軸の
# ラベルは元々縦なのに水平になったりする
set_default_plot_size(10cm, 9cm)
# データ読み込み、データ形式の確認、要約統計量
dat = CSV.read("datasets/SoaySheepFitness.csv", allowmissing=:none)
describe(dat)
names!(dat, [:fitness, :body_size])
# プロット、二つの実数変数だけのデータなので、散布図を見る
p = plot(dat, x = :body_size, y = :fitness, Geom.point)
draw(PDF("soay_1.pdf", 10cm, 8cm), p)
# まず、線形回帰して診断プロットを見て、前提条件が破綻している様子を見る
# まずモデリング、ｐ値は悪くないが
model = lm(@formula(fitness ~ body_size), dat)
pred = predict(model)
# 診断プロット１：Residuals vs Fitted
res = dat[:fitness] - pred
p = plot(x = pred, y = res, Geom.point,
Guide.title("Residuals vs Fitted"),
Guide.xlabel("Fitted values"), Guide.ylabel("Residuals"))
draw(PDF("soay_2.pdf", 10cm, 8cm), 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("soay_3.pdf", 10cm, 8cm), p)
# 診断プロット３：Scale-Location
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("soay_4.pdf", 10cm, 9cm), p)
# 診断プロット４：Residuals vs Leverage
x = dat[:body_size]
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("soay_5.pdf", 10cm, 8cm), p)
# データとモデルのプロット
# まずｘ軸の範囲を決める、データの範囲を探して､50点に区切る
# データフレームに入れて、回帰時の変数名と同じ変数名を付ける
# モデルは明らかにデータに当てはまってない（右に行くほど）
xmin = describe(dat)[2, :min]
xmax = describe(dat)[2, :max]
newx = DataFrame([range(xmin, stop = xmax, length = 50)])
names!(newx, [:body_size])
# で予測値を計算して、結果のデータフレームにｘの値を挿入
pred = predict(model, newx, interval = :confidence)
insertcols!(pred, 1, body_size=newx[:body_size])
# でデータと予測直線をプロット
# 先に置いた Geom.point が Geom.line の上に描画される
# ymin と ymax は Geom.ribbon が使う
p = plot(layer(dat, x = :body_size, y = :fitness, Geom.point),
layer(pred,x = :body_size, y = :prediction,
ymin = :lower, ymax = :upper,
Geom.line, Geom.ribbon))
draw(PDF("soay_6.pdf", 10cm, 8cm), p)
# ＧＬＭ ######################################################################
model = glm(@formula(fitness ~ body_size), dat, Poisson(), LogLink())
pred = predict(model)
res = dat[:fitness] .- pred # 逸脱度残差の符号を決めるのだけに使う
p = plot(layer(dat, x = :body_size, y = :fitness, Geom.point),
layer(x = dat[:body_size], y = pred, Geom.line))
draw(PDF("soay_7.pdf", 10cm, 8cm), 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("soay_8.pdf", 10cm, 8cm), 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("soay_9.pdf", 10cm, 8cm), p)
# 診断プロット３：Scale-Location
# GSwR2 とちょっと違うところがある
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("soay_10.pdf", 10cm, 8cm), 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[: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) # ハット行列の対角成分がレバレッジ
p = plot(x = leverage, y = stdDevRes, Geom.point,
Guide.title("Residuals vs Leverage"),
Guide.xlabel("Leverage"), Guide.ylabel("Std. Pearson resid."))
draw(PDF("soay_11.pdf", 10cm, 8cm), p)
# ということで、データとモデルと信頼区間のプロット
# まずｘ軸の範囲を決める、データの範囲を探して､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 が使う
p = 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("生育速度"))
draw(PDF("soay_12.pdf", 10cm, 8cm), p)
2019, © Daisuke TOMINAGA.