Julia でふつうの統計解析

線形回帰 Linear Regression

前のページに戻る 
#
# GSwR2 第５章 その３ 線形回帰
#
# 植物の生育速度のデータ
#
# 2019/03/08 Daisuke TOMINGAGA

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

# プロットサイズはあまり小さくすると軸ラベルが収まらなくなり、たとえばｙ軸の
# ラベルは元々縦なのに水平になったりする
set_default_plot_size(10cm, 8cm)

# データ読み込み
dat = CSV.read("datasets/plant.growth.rate.csv", allowmissing=:none)
# カラム名に . が入ってるとあとでうまく扱えないので､名前を変える
names!(dat, [:moisture, :growth]);

# データは要するに二群、各群でのオゾン濃度の平均値と分散
describe(dat)
var(dat[:moisture])
var(dat[:growth])

p = hstack(plot(dat, y = :moisture, Geom.boxplot),
           plot(dat, y = :growth, Geom.boxplot))
draw(PDF("plant_1.pdf", 10cm, 8cm), p)
p = plot(dat, x = :moisture, y = :growth, Geom.point)
draw(PDF("plant_2.pdf", 10cm, 8cm), p)

# 線形回帰して、モデルの値（期待値）を計算し、診断プロットを描く
model = lm(@formula(growth ~ moisture), dat)
pred = predict(model)

# 診断プロット１：Residuals vs Fitted
res = dat[:growth] - pred
p = plot(x = pred, y = res, Geom.point,
         Guide.title("Residuals vs Fitted"),
         Guide.xlabel("Fitted values"), Guide.ylabel("Residuals"))
draw(PDF("plant_3.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("plant_4.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("plant_5.pdf", 10cm, 8cm), p)

# 診断プロット４：Residuals vs Leverage
x = dat[:moisture]
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("plant_6.pdf", 10cm, 8cm), p)

# で、データとモデルと信頼区間のプロット
# まずｘ軸の範囲を決める、データの範囲を探して､50点に区切る
# データフレームに入れて、回帰時の変数名と同じ変数名を付けないと､
# predict が受け取ってくれない
xmin = describe(dat)[1, :min]
xmax = describe(dat)[1, :max]
newx = DataFrame([range(xmin, stop = xmax, length = 50)])
names!(newx, [:moisture])
# で予測値を計算して、結果のデータフレームにｘの値を挿入
pred = predict(model, newx, interval=:confidence)
insertcols!(pred, 1, moisture=newx[:moisture])
# でデータと予測直線をプロット
# 先に置いた Geom.point が Geom.line の上に描画される
# ymin と ymax は Geom.ribbon が使う
p = plot(layer(dat, x = :moisture, y = :growth, Geom.point),
           layer(pred,x = :moisture, y = :prediction,
                      ymin = :lower, ymax = :upper,
                      Geom.line, Geom.ribbon),
           Coord.cartesian(xmin=0.2, xmax=2.0),
           Guide.title("土壌水分量と生育速度"),
           Guide.xlabel("土壌水分量"),
           Guide.ylabel("生育速度"))
draw(PDF("plant_7.pdf", 10cm, 8cm), p)
2019, © Daisuke TOMINAGA.