読者です 読者をやめる 読者になる 読者になる

メモ:lmの信頼区間をbroom::augumentで出す

R broom

アヒル本の図4.3をpredict()を使わずに描くにはどうするのか気になって調べたときのメモ。ちなみにアヒル本が使っているコードはサポートページにあります。

github.com

まずはデータをとってきてlm()を使います。

d <- readr::read_csv("https://raw.githubusercontent.com/MatsuuraKentaro/RStanBook/master/chap04/input/data-salary.txt")
res_lm <- lm(Y ~ X, data = d)

このlmのオブジェクトをaugument()に渡します。predict()のように渡したデータに対する予測値を出すには、newdata引数に指定します。

library(broom)

res_aug <- augment(res_lm, newdata = data.frame(X = 23:60))
res_aug
X .fitted .se.fit
23 384.0995 35.59951
24 406.0037 34.29026
25 427.9079 32.99893
26 449.8121 31.72768
27 471.7163 30.47904
28 493.6205 29.25591
29 515.5247 28.06161
30 537.4289 26.89999
31 559.3331 25.77547
32 581.2373 24.69311
33 603.1415 23.65870
34 625.0457 22.67881
35 646.9499 21.76081
36 668.8541 20.91283
37 690.7583 20.14373
38 712.6625 19.46286
39 734.5667 18.87978
40 756.4709 18.40377
41 778.3751 18.04331
42 800.2793 17.80542
43 822.1835 17.69505
44 844.0877 17.71457
45 865.9919 17.86357
46 887.8961 18.13886
47 909.8003 18.53480
48 931.7045 19.04387
49 953.6087 19.65729
50 975.5129 20.36563
51 997.4171 21.15935
52 1019.3213 22.02924
53 1041.2255 22.96663
54 1063.1297 23.96362
55 1085.0339 25.01306
56 1106.9381 26.10864
57 1128.8423 27.24479
58 1150.7466 28.41665
59 1172.6508 29.61998
60 1194.5550 30.85110

ここから信頼区間と予測区間を計算するのは、以下のSOの答えに詳しく書いてありました。

数式的な理解がまったく追い付いていませんが、とりあえずやってみたときのメモ。

# alphaは真の値が区間に入らない確率(あってる?)。qtはt分布の分位点を求める関数、dfは自由度を指定する引数
alpha <- qt((1-0.95)/2, df = res_lm$df.residual)

library(ggplot2)
library(tidyverse)

# alphaを予測値の標準誤差.se.fitにかける
res_aug %>%
  mutate(upper = .fitted - .se.fit * alpha,
         lower = .fitted + .se.fit * alpha) %>%
  ggplot(aes(X)) +
  geom_ribbon(aes(ymax = upper, ymin = lower), fill = "blue", alpha = 0.2) +
  geom_line(aes(y = .fitted)) +
  theme_minimal()

f:id:yutannihilation:20170121184636p:plain