メモ:lmの信頼区間をbroom::augumentで出す
アヒル本の図4.3をpredict()
を使わずに描くにはどうするのか気になって調べたときのメモ。ちなみにアヒル本が使っているコードはサポートページにあります。
まずはデータをとってきて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()