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

アヒル本の図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()