kilometer
September 17, 2022
340

# TokyoR#101_RegressionAnalysis

## kilometer

September 17, 2022

## Transcript

3. ### Who！？ 名前： 三村 @kilometer 職業： ポスドク (こうがくはくし) 専⾨： ⾏動神経科学(霊⻑類) 脳イメージング

医療システム⼯学 R歴： ~ 10年ぐらい 流⾏: パン切り包丁

11. ### library(tidyverse) library(magrittr) set.seed(1) 準備 N <- 10 # データ数 a

<- 0.1 # 傾き dat <- tibble(x = c(1:N), y = a * x + rnorm(N))
12. ### > dat # A tibble: 10 × 2 x y

<int> <dbl> 1 1 -0.526 2 2 0.384 3 3 -0.536 4 4 2.00 5 5 0.830 6 6 -0.220 7 7 1.19 8 8 1.54 9 9 1.48 10 10 0.695 準備
13. ### 𝑦 = 𝑎! + 𝑎" 𝑥 + 𝑢 𝑢 ~

!.!.#. 𝒩(0, σ\$) (単)回帰モデル
14. ### 線形モデル linear model dat <- tibble(x = c(1:N), y =

a * x + rnorm(N)) lm(y ~ x, data = dat) Call: lm(formula = y ~ x, data = dat) Coefficients: (Intercept) x -0.1688 0.1547
15. ### fit_lm <- dat %>% lm(y ~ x, data = .)

> fit_lm %>% names() [1] "coefficients" "residuals" [3] "effects" "rank" [5] "fitted.values" "assign" [7] "qr" "df.residual" [9] "xlevels" "call" [11] "terms" "model" 線形モデル linear model
16. ### > fit_lm %>% summary() Call: lm(formula = y ~ x,

data = .) Residuals: Min 1Q Median 3Q Max -0.9800 -0.6410 0.2338 0.2678 1.5452 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.16882 0.55270 -0.305 0.768 x 0.15473 0.08908 1.737 0.121 Residual standard error: 0.8091 on 8 degrees of freedom Multiple R-squared: 0.2739, Adjusted R-squared: 0.1831 F-statistic: 3.018 on 1 and 8 DF, p-value: 0.1206 線形モデル linear model

20. ### 例えば、 > fit_lm %>% summary() Call: lm(formula = y ~

x, data = .) Residuals: Min 1Q Median 3Q Max -0.9800 -0.6410 0.2338 0.2678 1.5452 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.16882 0.55270 -0.305 0.768 x 0.15473 0.08908 1.737 0.121 Residual standard error: 0.8091 on 8 degrees of freedom Multiple R-squared: 0.2739, Adjusted R-squared: 0.1831 F-statistic: 3.018 on 1 and 8 DF, p-value: 0.1206 この中に誤記があります(知ってます？)。
21. ### 今⽇のミッション > fit_lm %>% summary() Call: lm(formula = y ~

x, data = .) Residuals: Min 1Q Median 3Q Max -0.9800 -0.6410 0.2338 0.2678 1.5452 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.16882 0.55270 -0.305 0.768 x 0.15473 0.08908 1.737 0.121 Residual standard error: 0.8091 on 8 degrees of freedom Multiple R-squared: 0.2739, Adjusted R-squared: 0.1831 F-statistic: 3.018 on 1 and 8 DF, p-value: 0.1206 この中に書かれた数字を全て計算する。
22. ### > fit_lm %>% summary() Call: lm(formula = y ~ x,

data = .) Residuals: Min 1Q Median 3Q Max -0.9800 -0.6410 0.2338 0.2678 1.5452 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.16882 0.55270 -0.305 0.768 x 0.15473 0.08908 1.737 0.121 Residual standard error: 0.8091 on 8 degrees of freedom Multiple R-squared: 0.2739, Adjusted R-squared: 0.1831 F-statistic: 3.018 on 1 and 8 DF, p-value: 0.1206 1. 残差 residuals
23. ### fit_lm <- dat %>% lm(y ~ x, data = .)

> fit_lm %>% names() [1] "coefficients" "residuals" [3] "effects" "rank" [5] "fitted.values" "assign" [7] "qr" "df.residual" [9] "xlevels" "call" [11] "terms" "model" 1. 残差 residuals
24. ### res <- fit_lm\$residuals > res 1 2 3 -0.5123623 0.2430028

-0.8310012 4 5 6 1.5451761 0.2246710 -0.9800372 7 8 9 0.2731282 0.4692918 0.2520164 10 -0.6838854 1. 残差 residuals
25. ### res <- fit_lm\$residuals > res %>% median() [1] 0.2187553 >

res %>% min() [1] -2.390729 > res %>% max() [1] 1.481385 > res %>% + quantile(probs = c(0.25, 0.75)) 25% 75% -0.3329012 0.5308412 1. 残差 residuals
26. ### > fit_lm %>% summary() Call: lm(formula = y ~ x,

data = .) Residuals: Min 1Q Median 3Q Max -0.9800 -0.6410 0.2338 0.2678 1.5452 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.16882 0.55270 -0.305 0.768 x 0.15473 0.08908 1.737 0.121 Residual standard error: 0.8091 on 8 degrees of freedom Multiple R-squared: 0.2739, Adjusted R-squared: 0.1831 F-statistic: 3.018 on 1 and 8 DF, p-value: 0.1206 1. 残差 residuals
27. ### res <- fit_lm\$residuals 残差εは ε = 𝑦 − ) 𝑦

データ 推定値 計算してみましょう。 1. 残差 residuals
28. ### a0 <- fit_lm\$coefficients[1] a1 <- fit_lm\$coefficients[2] dat %>% mutate(y_hat =

a0 + a1 * x, res = y - y_hat) # A tibble: 10 × 4 x y y_hat res <int> <dbl> <dbl> <dbl> 1 1 -0.526 -0.0141 -0.512 2 2 0.384 0.141 0.243 3 3 -0.536 0.295 -0.831 1. 残差 residuals

30. ### dat_fit <- dat %>% group_nest() %>% mutate(fit = map(data, ~

lm(y ~ x, data = .))) > dat_fit # A tibble: 1 × 2 data fit <list> <list> 1 <tibble [10 × 2]> <lm> 「線形回帰？nestでしょう。」 1. 残差 residuals
31. ### dat_res <- dat_fit %>% mutate(y_hat = map( fit, predict)) %>%

select(-fit) %>% unnest(everything()) %>% mutate(res = y - y_hat) 1. 残差 residuals
32. ### > dat_res # A tibble: 10 × 4 x y

y_hat res <int> <dbl> <dbl> <dbl> 1 1 -0.526 -0.0141 -0.512 2 2 0.384 0.141 0.243 3 3 -0.536 0.295 -0.831 4 4 2.00 0.450 1.55 1. 残差 residuals
33. ### dat_res %>% summarize( across(res, list(Min = ~ min(.), Q1 =

~ quantile(., probs = c(0.25)), Median = ~ median(.), Q3 = ~ quantile(., probs = c(0.75)), Max = ~ max(.)))) # A tibble: 1 × 5 res_Min res_Q1 res_Median res_Q3 res_Max <dbl> <dbl> <dbl> <dbl> <dbl> 1 -0.980 -0.641 0.234 0.268 1.55 1. 残差 residuals
34. ### > fit_lm %>% summary() Call: lm(formula = y ~ x,

data = .) Residuals: Min 1Q Median 3Q Max -0.9800 -0.6410 0.2338 0.2678 1.5452 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.16882 0.55270 -0.305 0.768 x 0.15473 0.08908 1.737 0.121 Residual standard error: 0.8091 on 8 degrees of freedom Multiple R-squared: 0.2739, Adjusted R-squared: 0.1831 F-statistic: 3.018 on 1 and 8 DF, p-value: 0.1206 1. 残差 residuals
35. ### > dat_res # A tibble: 10 × 4 x y

y_hat res <int> <dbl> <dbl> <dbl> 1 1 -0.526 -0.0141 -0.512 2 2 0.384 0.141 0.243 3 3 -0.536 0.295 -0.831 4 4 2.00 0.450 1.55 1. 残差 residuals
36. ### ggplot(data = dat_res) + aes(x = x, y = y)

+ geom_path(aes(y = y_hat), color = "red", linetype = "dotted") + geom_segment(aes(xend = x, yend = y_hat), color = "red") + geom_point() 1. 残差 residuals
37. ### > fit_lm %>% summary() Call: lm(formula = y ~ x,

data = .) Residuals: Min 1Q Median 3Q Max -0.9800 -0.6410 0.2338 0.2678 1.5452 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.16882 0.55270 -0.305 0.768 x 0.15473 0.08908 1.737 0.121 Residual standard error: 0.8091 on 8 degrees of freedom Multiple R-squared: 0.2739, Adjusted R-squared: 0.1831 F-statistic: 3.018 on 1 and 8 DF, p-value: 0.1206 線形モデル linear model
38. ### 標準誤差 standard error 𝑆𝐸 = 𝑆𝐷 𝑛 𝑆𝐷 = 1

𝑛 − 1 ( !"# \$ (𝑦! − + 𝑦)% 標本平均の標準偏差 （の推定値）
39. ### Residual standard error: 0.8091 dat_res %>% summarize_at("res", ~ sd(.) /

sqrt(N)) # A tibble: 1 × 1 res <dbl> 1 0.241
40. ### dat_res %>% summarize_at("res", ~ sd(.) / sqrt(N)) Residual standard error:

0.8091 # A tibble: 1 × 1 res <dbl> 1 0.241 !!?
41. ### Residual standard error: 0.8091 dat_res %>% summarize_at("res", ~ sqrt(sum(.^2) /

(N - 2))) # A tibble: 1 × 1 res <dbl> 1 0.809
42. ### Residual standard error: 0.8091 dat_res %>% summarize_at("res", ~ sqrt(sum(.^2) /

(N - 2))) # A tibble: 1 × 1 res <dbl> 1 0.809 deviation
43. ### 𝑅𝑆𝐷 = 1 𝑛 − 𝑘 − 1 ) !"#

\$ (𝑦! − , 𝑦)% 残差標準偏差 Residual standard deviation 説明変数の次元（今は𝑘 = 1） 誤差𝑢の標準偏差の期待値
44. ### > ?stats::sigma Description Extract the estimated standard deviation of the

errors, the “residual standard deviation” (misnamed also “residual standard error”, e.g., in summary.lm()'s output, from a fitted model). 2. 残差標準偏差 RSD
45. ### > fit_lm %>% summary() Call: lm(formula = y ~ x,

data = .) Residuals: Min 1Q Median 3Q Max -0.9800 -0.6410 0.2338 0.2678 1.5452 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.16882 0.55270 -0.305 0.768 x 0.15473 0.08908 1.737 0.121 Residual standard error: 0.8091 on 8 degrees of freedom Multiple R-squared: 0.2739, Adjusted R-squared: 0.1831 F-statistic: 3.018 on 1 and 8 DF, p-value: 0.1206 deviation 2. 残差標準偏差 RSD
46. ### 3. 決定係数 R-squared > fit_lm %>% summary() Call: lm(formula =

y ~ x, data = .) Residuals: Min 1Q Median 3Q Max -0.9800 -0.6410 0.2338 0.2678 1.5452 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.16882 0.55270 -0.305 0.768 x 0.15473 0.08908 1.737 0.121 Residual standard error: 0.8091 on 8 degrees of freedom Multiple R-squared: 0.2739, Adjusted R-squared: 0.1831 F-statistic: 3.018 on 1 and 8 DF, p-value: 0.1206 deviation
47. ### 3つの変動 残差変動 RSS 全変動 TSS 回帰変動 ESS 𝑅𝑆𝑆 = &

!"# \$ (𝑦! − * 𝑦)% 𝑇𝑆𝑆 = & !"# \$ (𝑦! − - 𝑦)% 𝐸𝑆𝑆 = & !"# \$ (* 𝑦 − - 𝑦)% データ 推定値 データ 平均 推定値 平均
48. ### 3つの変動 残差変動 RSS 全変動 TSS 回帰変動 ESS 𝑅𝑆𝑆 = &

!"# \$ (𝑦! − * 𝑦)% 𝑇𝑆𝑆 = & !"# \$ (𝑦! − - 𝑦)% 𝐸𝑆𝑆 = & !"# \$ (* 𝑦 − - 𝑦)% データ 推定値 データ 平均 推定値 平均 = -
49. ### 𝑅2 = 𝐸𝑆𝑆 𝑇𝑆𝑆 = 1 − 𝑅𝑆𝑆 𝑇𝑆𝑆 全変動

回帰変動 残差変動 3. 決定係数 R-squared 決定係数
50. ### dat_dev <- dat_res %>% mutate(yhat_mean = mean(y_hat), RSS = y

- y_hat, ESS = y_hat - yhat_mean, TSS = y - yhat_mean) %>% summarize_at(c("RSS", "ESS", "TSS"), ~ sum(. ^ 2)) 3つの変動
51. ### > dat_dev # A tibble: 1 × 3 RSS ESS

TSS <dbl> <dbl> <dbl> 1 5.24 1.98 7.21 > dat_dev %\$% {ESS / TSS} [1] 0.2738825 3つの変動
52. ### 4. ⾃由度とf検定 > fit_lm %>% summary() Call: lm(formula = y

~ x, data = .) Residuals: Min 1Q Median 3Q Max -0.9800 -0.6410 0.2338 0.2678 1.5452 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.16882 0.55270 -0.305 0.768 x 0.15473 0.08908 1.737 0.121 Residual standard error: 0.8091 on 8 degrees of freedom Multiple R-squared: 0.2739, Adjusted R-squared: 0.1831 F-statistic: 3.018 on 1 and 8 DF, p-value: 0.1206 deviation

54. ### F分布 𝑌 = 𝑋! /𝑚 𝑋" /𝑛 互いに独⽴な𝑋! ~𝜒"(𝑚)、𝑋" ~𝜒"(𝑛)について、

なる𝑌の確率分布をF分布と呼び、 𝑌~𝐹(𝑚, 𝑛) と表す。 (𝑚を第1⾃由度,𝑛を第2⾃由度と呼ぶ) カイ2乗分布
55. ### 𝜒!分布 標準正規分布𝒩(0,1)に従う独⽴な𝑘個の 確率変数𝑋! , 𝑋" , … , 𝑋# について、

𝑍 = 1 12! 3 𝑋1 " なる 𝑍 の確率分布を 𝜒" 分布と呼び、 𝑍~𝜒"(𝑘) と表す。(𝑘を⾃由度と呼ぶ)
56. ### 𝜒!分布 正規分布𝒩(𝑢, σ")に従う独⽴な𝑘個の 確率変数𝑋! , 𝑋" , … , 𝑋#

について、 𝑍 = 1 12! 3 (𝑋1 −𝑢)" σ" なる𝑍は⾃由度𝑘 − 1の𝜒"分布に従う 𝑍~𝜒"(𝑘 − 1)
57. ### 正規分布に従う独⽴な𝑘個の確率変数 𝑋! , 𝑋" , … , 𝑋# の2乗和、すなわちコレ！？ データ数𝑘の時、残差⼆乗和は

⾃由度𝑘 − 1のχ2分布に従う？
58. ### 正規分布に従う独⽴な𝑘個の確率変数 𝑋! , 𝑋" , … , 𝑋# の2乗和、すなわちコレ！？ データ数𝑘の時、残差⼆乗和は

⾃由度𝑘 − 1のχ2分布に従う？ → 従わない

− 1のχ2分布に従う。

61. ### ⾃由度 degrees of freedom 変数のうち⾃由な値をとりうるものの数 ⺟集団 標本𝑥 (𝑛個) ⾃由に選ばれた𝑛個の標本の⾃由度は𝒏で、 標本の平均値は標本平均

̅ 𝑥。 では、標本平均が既知の値 ̅ 𝑥になるように 選ばれた𝑛個の標本の⾃由度は？
62. ### ⾃由度 degrees of freedom 変数のうち⾃由な値をとりうるものの数 ⺟集団 標本𝑥 (𝑛個) ⾃由に選ばれた𝑛個の標本の⾃由度は𝒏で、 標本の平均値は標本平均

̅ 𝑥。 では、標本平均が既知の値 ̅ 𝑥になるように 選ばれた𝑛個の標本の⾃由度は？ → 𝒏 − 𝟏
63. ### ⾃由度 degrees of freedom 変数のうち⾃由な値をとりうるものの数 標本平均が既知の値 ̅ 𝑥になるように 選ばれた𝑛個の標本の⾃由度は 𝒏

− 𝟏 ̅ 𝑥 = 1 𝑛 (𝑥! + 𝑥" + ⋯ + 𝑥#\$! + 𝑥# ) 𝑥# = 𝑛 ̅ 𝑥 − (𝑥! + 𝑥" + ⋯ + 𝑥#\$! ) ⾃由に値をとる𝒏 − 𝟏項 既知の 標本平均

66. ### 𝑦 = 𝑎! + 𝑎" 𝑥 + 𝑢 𝑢 ~

!.!.#. 𝒩(0, σ\$) (単)回帰モデル

ε 𝑢 𝑦 𝑥 𝑌
69. ### 𝑦 = 𝑎! + 𝑎" 𝑥 + 𝑢 𝑢 ~

!.!.#. 𝒩(0, σ\$) 誤差 𝑦 = ) 𝑎! + ) 𝑎" 𝑥 + ε ε = 𝑦 − ) 𝑦 真の回帰モデル 残差 最⼩⼆乗法で推定された回帰式
70. ### (正規)回帰モデル5つの仮定 1. 誤差の期待値はゼロ 2. 誤差の分散は定数 3. 誤差は互いに無相関 4. 誤差と説明変数は無相関 5.

誤差は正規分布に従う 𝑢 ~ %.%.'. 𝒩(0, σ") 𝐸[𝑢% (𝑥% − 𝐸[𝑥% ])] = 0 𝐸[𝑢% 𝑢( ] = 0 (𝑖 ≠ 𝑗) 𝑉𝑎𝑟 𝑢 = σ" 𝐸 𝑢 = 0
71. ### (正規)回帰モデル5つの仮定 1. 誤差の期待値はゼロ 2. 誤差の分散は定数 3. 誤差は互いに無相関 4. 誤差と説明変数は無相関 ガウス=マルコフの定理

仮定1-4が満たされる場合、最⼩⼆乗推定量は 線型不偏推定量の中で最⼩分散を持つ。 5. 誤差は正規分布に従う + 仮定5 最⼩⼆乗推定量は最⼩分散を持つ。
72. ### 𝑦 = 𝑎! + 𝑎" 𝑥 + 𝑢 𝑢 ~

!.!.#. 𝒩(0, σ\$) 誤差 真の回帰モデル 最⼩⼆乗法OLSでモデルの パラメータを推定する (OLS, Ordinary Last Squares) → 残差⼆乗和RSSの最⼩化 (RSS, Residual Sum of Squares)
73. ### OLSが要請する残差の性質 1. 残差の合計値はゼロ 2. 残差と説明変数は無相関 ; %)! # ε% =

0 ; %)! # 𝑥% ε% = 0
74. ### 1. 残差の合計値はゼロ 2. 残差と説明変数は無相関 ; %)! # ε% = 0

; %)! # 𝑥*% ε% = 0 (𝑘 = 1, 2, … ) 説明変数の次元数 OLSが要請する残差の性質
75. ### 𝑦 = / 𝑎& + / 𝑎# 𝑥# + ⋯

+ / 𝑎' 𝑥' + ε 残差 OLSで推定された(重)回帰式 ; %)! # ε% = 0 ; %)! # 𝑥*% ε% = 0 (𝑘 = 1, 2, … ) かつ 𝑛個の残差εは、 𝑘 + 1個の式を満たさなければならない。 残差εの⾃由度はn − 𝑘 − 1となる。 𝑘項

77. ### 𝑅𝑆𝐷 = 1 𝑛 − 𝑘 − 1 ) !"#

\$ (𝑦! − , 𝑦)% 残差標準偏差 Residual standard deviation 説明変数の次元（今は𝑘 = 1） 𝑦 = 𝑎+ + 𝑎! 𝑥 + 𝑢
78. ### 1. 残差の合計値はゼロ 2. 残差と説明変数は無相関 ; %)! # ε% = 0

; %)! # 𝑥*% ε% = 0 (𝑘 = 1, 2, … ) 説明変数の次元数 OLSが要請する残差の性質 ※ 定数項/ 𝒂𝟎 が無い回帰式では この性質は満たされない。 → 残差𝜺の⾃由度は𝐧 − 𝒌。
79. ### 𝜒!分布 正規分布に従う独⽴な𝑘個の確率変数 𝑋! , 𝑋" , … , 𝑋# について、

𝑍 = 1 12! 3 (𝑋1 − 8 𝑋)" σ" なる𝑍は⾃由度𝑘 − 1の𝜒"分布に従う 𝑍~𝜒"(𝑘 − 1)

1)

𝑦! − - 𝑦 = 0
82. ### 回帰変動は(理想的には) ⾃由度𝑘のχ2分布に従う。 𝐸𝑆𝑆 = 𝑇𝑆𝑆 − 𝑅𝑆𝑆 𝜒!(𝑛 − 1)

𝜒!(𝑛 − 𝑘 − 1) 𝜒!(𝑘) これらが互いに独⽴なので、 & !"# \$ * 𝑦! − - 𝑦 = 0
83. ### 𝑅𝑆𝑆~𝜒((𝑛 − 𝑘 − 1) 𝑇𝑆𝑆~𝜒((𝑛 − 1) 𝐸𝑆𝑆~𝜒((𝑘) OLSの変動の確率分布

残差変動 全変動 回帰変動
84. ### 回帰モデルの適合性 𝑅2 = 𝐸𝑆𝑆 𝑇𝑆𝑆 = 1 − 𝑅𝑆𝑆 𝑇𝑆𝑆

𝑅2 = 1 − 𝑅𝑆𝑆/(𝑛 − 𝑘 − 1) 𝑇𝑆𝑆/(𝑛 − 1) 決定係数 ⾃由度調整済み決定係数
85. ### > dat_dev # A tibble: 1 × 3 RSS ESS

TSS <dbl> <dbl> <dbl> 1 10.1 1.82 11.9 3つの変動 > dat_dev %\$% {ESS / TSS} [1] 0.2738825 決定係数 > dat_dev %\$% {1 - RSS * (N - 1) / TSS / (N - 1 - 1)} [1] 0.1831179
86. ### 4. ⾃由度とf検定 > fit_lm %>% summary() Call: lm(formula = y

~ x, data = .) Residuals: Min 1Q Median 3Q Max -0.9800 -0.6410 0.2338 0.2678 1.5452 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.16882 0.55270 -0.305 0.768 x 0.15473 0.08908 1.737 0.121 Residual standard error: 0.8091 on 8 degrees of freedom Multiple R-squared: 0.2739, Adjusted R-squared: 0.1831 F-statistic: 3.018 on 1 and 8 DF, p-value: 0.1206 deviation
87. ### F分布 𝑌 = 𝑋! /𝑚 𝑋" /𝑛 互いに独⽴な𝑋! ~𝜒"(𝑚)、𝑋" ~𝜒"(𝑛)について、

なる𝑌の確率分布をF分布と呼び、 𝑌~𝐹(𝑚, 𝑛) と表す。 (𝑚を第1⾃由度,𝑛を第2⾃由度と呼ぶ) カイ2乗分布
88. ### 𝑓 = 𝐸𝑆𝑆/𝑘 𝑅𝑆𝑆/(𝑛 − 𝑘 − 1) F値 無相関(傾き=0)なら0

分散0なら0 𝑓~𝐹(𝑘, 𝑛 − 𝑘 − 1) 4. ⾃由度とf検定
89. ### > dat_dev # A tibble: 1 × 3 RSS ESS

TSS <dbl> <dbl> <dbl> 1 10.1 1.82 11.9 3つの変動 > df1 <- dat %>% select(contains("x")) %>% ncol() > df2 <- N – df1 - 1 F値 > dat_dev %\$% {ESS * df2 / RSS / df1} [1] 3.017501
90. ### > df1 <- dat %>% select(contains("x")) %>% ncol() > df2

<- N – df1 - 1 > f <- dat_dev %\$% {ESS * df2 / RSS / df1} [1] 3.017501 F値 > pf(q = f, df1 = df1, df2 = df2, lower.tail = FALSE) [1] 0.264185 P値
91. ### 回帰モデルの適合性 F(1,8) f = 3.017501 p = 0.264185 tibble(x =

seq(0, 5, by = 0.01), df = df(x, df1, df2)) lower.tail=TRUE lower.tail=FALSE
92. ### 回帰 dat %>% lm(y ~ x, data = .) %>%

summary() %>% names() [1] "call" "terms" "residuals" [4] "coefficients" "aliased" "sigma" [7] "df" "r.squared" "adj.r.squared" [10] "fstatistic" "cov.unscaled"
93. ### 回帰 dat %>% lm(y ~ x, data = .) %>%

summary() %\$% fstatistic value numdf dendf 3.017501 1.000000 8.000000
94. ### 回帰 dat %>% lm(y ~ x, data = .) %>%

summary() %\$% fstatistic %>% { pf(.[1], .[2], .[3], lower.tail = FALSE) } value 0.1205754
95. ### 回帰 lm_pval <- function(fit){ summary(fit) %\$% fstatistic %>% { pf(.[1],

.[2], .[3], lower.tail = FALSE)} } dat_fit %>% mutate(pval = map_dbl(fit, lm_pval)) # A tibble: 1 × 3 data fit pval <list> <list> <dbl> 1 <tibble [10 × 2]> <lm> 0.121
96. ### > fit_lm %>% summary() Call: lm(formula = y ~ x,

data = .) Residuals: Min 1Q Median 3Q Max -0.9800 -0.6410 0.2338 0.2678 1.5452 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.16882 0.55270 -0.305 0.768 x 0.15473 0.08908 1.737 0.121 Residual standard error: 0.8091 on 8 degrees of freedom Multiple R-squared: 0.2739, Adjusted R-squared: 0.1831 F-statistic: 3.018 on 1 and 8 DF, p-value: 0.1206 deviation 4. ⾃由度とf検定
97. ### > fit_lm %>% summary() Call: lm(formula = y ~ x,

data = .) Residuals: Min 1Q Median 3Q Max -0.9800 -0.6410 0.2338 0.2678 1.5452 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.16882 0.55270 -0.305 0.768 x 0.15473 0.08908 1.737 0.121 Residual standard error: 0.8091 on 8 degrees of freedom Multiple R-squared: 0.2739, Adjusted R-squared: 0.1831 F-statistic: 3.018 on 1 and 8 DF, p-value: 0.1206 deviation 5. 回帰係数
98. ### 最⼩⼆乗法（線形回帰） argmin(-!,-") ; %)! # ε" 𝑦 = / 𝑎&

+ / 𝑎# 𝑥 + ε 残差 & 𝑎" = ∑#\$" % (𝑥# − ̅ 𝑥)(𝑦# − / 𝑦) ∑ #\$" % (𝑥# − ̅ 𝑥)! & 𝑎& = / 𝑦 − & 𝑎" ̅ 𝑥
99. ### 𝑉𝑎𝑟[𝑥] = 1 𝑛 − 1 6 )*+ , (𝑥)

− ̅ 𝑥)( 𝐶𝑜𝑣[𝑥, 𝑦] = 1 𝑛 − 1 6 )*+ , (𝑥) − ̅ 𝑥)(𝑦) − = 𝑦) 最⼩⼆乗法（線形回帰） < 𝑎! = ∑\$%! & (𝑥\$ − ̅ 𝑥)(𝑦\$ − A 𝑦) ∑ \$%! & (𝑥\$ − ̅ 𝑥)" = 𝐶𝑜𝑣[𝑥, 𝑦] 𝑉𝑎𝑟[𝑥] < 𝑎' = A 𝑦 − < 𝑎! ̅ 𝑥 = A 𝑦 − 𝐶𝑜𝑣[𝑥, 𝑦] 𝑉𝑎𝑟[𝑥] ̅ 𝑥 ただし、 分散 共分散
100. ### dat_g <- dat_res %>% mutate(x_bar = mean(x), y_bar = mean(y))

ggplot(dat_varcov) + aes(x, y) + geom_rect(data = dat_g %>% filter((x - x_bar)*(y - y_bar) >= 0), aes(xmin = x, xmax = x_bar, ymin = y, ymax = y_bar), alpha = 0.2, fill = "red") + geom_rect(data = dat_g %>% filter((x - x_bar)*(y - y_bar) < 0), aes(xmin = x, xmax = x_bar, ymin = y, ymax = y_bar), alpha = 0.2, fill = "blue") + geom_path(aes(y = y_hat)) + geom_point() 最⼩⼆乗法（線形回帰） 𝐶𝑜𝑣[𝑥, 𝑦] = 1 𝑛 − 1 6 )*+ , (𝑥) − ̅ 𝑥)(𝑦) − = 𝑦) 共分散
101. ### 最⼩⼆乗法（線形回帰） 𝐶𝑜𝑣[𝑥, 𝑦] = 1 𝑛 − 1 6 )*+

, (𝑥) − ̅ 𝑥)(𝑦) − = 𝑦) 共分散
102. ### 𝑉𝑎𝑟[𝑥] = 1 𝑛 − 1 & !"# \$ (𝑥!

− ̅ 𝑥)% 𝐶𝑜𝑣[𝑥, 𝑦] = 1 𝑛 − 1 & !"# \$ (𝑥! − ̅ 𝑥)(𝑦! − - 𝑦) # A tibble: 1 × 2 var cov <dbl> <dbl> 1 9.17 1.42 dat_varcov <- dat_res %>% mutate(x_bar = mean(x), y_bar = mean(y), x_d = x - x_bar, y_d = y - y_bar) %>% summarise(var = sum(x_d^2) / (N - 1), cov = sum(x_d * y_d) / (N - 1)) 最⼩⼆乗法（線形回帰） 分散 共分散
103. ### # A tibble: 1 × 2 var cov <dbl> <dbl>

1 9.17 1.42 dat_varcov 最⼩⼆乗法（線形回帰） > dat %\$% var(x) [1] 9.166667 > dat %\$% cov(x, y) [1] 1.418377 𝑉𝑎𝑟[𝑥] = 1 𝑛 − 1 & !"# \$ (𝑥! − ̅ 𝑥)% 𝐶𝑜𝑣[𝑥, 𝑦] = 1 𝑛 − 1 & !"# \$ (𝑥! − ̅ 𝑥)(𝑦! − - 𝑦) 分散 共分散
104. ### 最⼩⼆乗法（線形回帰） > 𝑎+ = 𝐶𝑜𝑣[𝑥, 𝑦] 𝑉𝑎𝑟[𝑥] , > 𝑎-

= = 𝑦 − 𝐶𝑜𝑣[𝑥, 𝑦] 𝑉𝑎𝑟 𝑥 ̅ 𝑥 a1 <- dat_varcov %\$% {cov / var} a0 <- dat %>% summarise(x_bar = mean(x), y_bar = mean(y)) %\$% {y_bar - a1 * x_bar} > a1 [1] 0.1547321 > a0 [1] -0.1688236
105. ### 5. 回帰係数のt検定 > fit_lm %>% summary() Call: lm(formula = y

~ x, data = .) Residuals: Min 1Q Median 3Q Max -0.9800 -0.6410 0.2338 0.2678 1.5452 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.16882 0.55270 -0.305 0.768 x 0.15473 0.08908 1.737 0.121 Residual standard error: 0.8091 on 8 degrees of freedom Multiple R-squared: 0.2739, Adjusted R-squared: 0.1831 F-statistic: 3.018 on 1 and 8 DF, p-value: 0.1206 deviation

107. ### 最⼩⼆乗法（線形回帰） 𝑦 = 2 𝑎? + 2 𝑎@ 𝑥 +

ε 回帰係数< 𝑎' と< 𝑎! の仮説検定を⾏う。 帰無仮説：< 𝑎! = 0 対⽴仮説：< 𝑎! ≠ 0 すなわち例えば< 𝑎! について、 → < 𝑎' と< 𝑎! の確率分布を知りたい。 → < 𝑎' と< 𝑎! の期待値と分散を調べる。
108. ### 𝑦 = 𝑎! + 𝑎" 𝑥 + 𝑢 𝑢 ~

!.!.#. 𝒩(0, σ\$) 誤差 𝑦 = ) 𝑎! + ) 𝑎" 𝑥 + ε ε = 𝑦 − ) 𝑦 真の回帰モデル 残差 最⼩⼆乗法で推定された回帰式
109. ### 期待と分散の基本的性質 𝐸 𝑋 + 𝑌 = 𝐸 𝑋 + 𝐸[𝑌]

𝐸 𝑐𝑋 = 𝑐𝐸 𝑋 𝑋と𝑌が独⽴の時、𝐸 𝑋𝑌 = 𝐸 𝑋 𝐸 𝑌 𝑉𝑎𝑟 𝑋 = 𝐸 (𝑋 − 𝐸[𝑋])" = 𝐸 𝑋" − (𝐸[𝑋])" 𝑉𝑎𝑟 𝑋 + 𝑌 = 𝑉𝑎𝑟 𝑋 + 𝑉𝑎𝑟 𝑌 + 2𝐶𝑜𝑣[𝑋, 𝑌] 𝑉𝑎𝑟 𝑐𝑋 = 𝑐"𝑉𝑎𝑟 𝑋 𝑉𝑎𝑟 𝑋 + 𝑐 = 𝑉𝑎𝑟 𝑋
110. ### > 𝑎+ = ∑)*+ , (𝑥) − ̅ 𝑥)(𝑦) −

= 𝑦) ∑ )*+ , (𝑥) − ̅ 𝑥)( = ∑)*+ , (𝑥) − ̅ 𝑥)( 𝑎- + 𝑎+ 𝑥) + 𝑢) − 𝑎- + 𝑎+ ̅ 𝑥 + = 𝑢 ) ∑ )*+ , (𝑥) − ̅ 𝑥)( = ∑)*+ , 𝑥) − ̅ 𝑥 (𝑎+ 𝑥) − ̅ 𝑥 + 𝑢) − = 𝑢) ∑ )*+ , (𝑥) − ̅ 𝑥)( = 𝑎+ + ∑)*+ , 𝑥) − ̅ 𝑥 𝑢) ∑ )*+ , (𝑥) − ̅ 𝑥)( + ∑)*+ , 𝑥) − ̅ 𝑥 ∑ )*+ , (𝑥) − ̅ 𝑥)( = 𝑢 𝐸[F 𝑎! ] = 𝑎! 従って、 …式(1) 5. 回帰係数のt検定
111. ### 𝑉𝑎𝑟 < 𝑎! = 𝐸 (< 𝑎! − 𝐸[< 𝑎!

])" = 𝐸 (< 𝑎! − 𝑎! )" = 𝐸 ∑!"# \$ 𝑥! − ̅ 𝑥 𝑢! ∑ !"# \$ (𝑥! − ̅ 𝑥)% % = 𝐸 ∑!"# \$ 𝑥! − ̅ 𝑥 %𝑢! % ∑ !"# \$ (𝑥! − ̅ 𝑥)% % + 𝐸 ∑!"# \$ ∑&"# \$ 𝑥! − ̅ 𝑥 𝑥& − ̅ 𝑥 𝑢! 𝑢& ∑ !"# \$ (𝑥! − ̅ 𝑥)% % 𝑖 ≠ 𝑗 = ∑!"# \$ 𝑥! − ̅ 𝑥 %𝐸 𝑢! % ∑ !"# \$ (𝑥! − ̅ 𝑥)% % + ∑!"# \$ ∑&"# \$ 𝑥! − ̅ 𝑥 𝑥& − ̅ 𝑥 𝐸[𝑢! 𝑢& ] ∑ !"# \$ (𝑥! − ̅ 𝑥)% % = 𝐸 𝑢% − (𝐸[𝑢])% ∑ !"# \$ (𝑥! − ̅ 𝑥)% = 𝑉𝑎𝑟[𝑢%] ∑ !"# \$ (𝑥! − ̅ 𝑥)% = σ% ∑ !"# \$ (𝑥! − ̅ 𝑥)% 𝐸 𝑢 = 0 𝑉𝑎𝑟 𝑎! = 𝐸 𝑎! − (𝐸 𝑎 )! 5. 回帰係数のt検定
112. ### 𝐸 < 𝑎! = 𝑎! 𝑉𝑎𝑟[< 𝑎! ] = σ"

∑ \$%! & (𝑥\$ − ̅ 𝑥)" 頑張って計算した結果、 > 𝑎+ = 𝑎+ + ∑)*+ , 𝑥) − ̅ 𝑥 ∑ )*+ , (𝑥) − ̅ 𝑥)( = 𝑢 + ∑)*+ , 𝑥) − ̅ 𝑥 𝑢) ∑ )*+ , (𝑥) − ̅ 𝑥)( また式(1) 𝑢~ 𝒩(0, σ() より * 𝑎# は正規分布に従い、 > 𝑎+ ~ 𝒩 𝑎+ , σ( ∑ )*+ , (𝑥) − ̅ 𝑥)( ⇔ > 𝑎+ − 𝑎+ σ(/ ∑ )*+ , (𝑥) − ̅ 𝑥)( ~𝒩 0,1 5. 回帰係数のt検定
113. ### F 𝑎! − 𝑎! σ"/ ∑ %)! # (𝑥% −

̅ 𝑥)" ~𝒩 0,1 コレに関する帰無仮説を検証する コレをどうにかしたい(データから推定したい) σは誤差𝑢の標準偏差 5. 回帰係数のt検定
114. ### 𝑅𝑆𝐷 = 1 𝑛 − 𝑘 − 1 ) !"#

\$ (𝑦! − , 𝑦)% 残差標準偏差 Residual standard deviation 誤差𝑢の標準偏差の期待値
115. ### F 𝑎! − 𝑎! 𝑅𝑆𝐷"/ ∑ %)! # (𝑥% −

̅ 𝑥)" = F 𝑎! − 𝑎! 𝑠𝑒(F 𝑎! ) ~ ? F 𝑎! − 𝑎! σ"/ ∑ %)! # (𝑥% − ̅ 𝑥)" ~𝒩 0,1 置き換える 確率分布は？ 𝑠𝑒 > 𝑎+ ∶= 𝑅𝑆𝐷( ∑ )*+ , (𝑥) − ̅ 𝑥)( 5. 回帰係数のt検定
116. ### F 𝑎! − 𝑎! 𝑠𝑒(F 𝑎! ) ~ ? F

𝑎! − 𝑎! σ"/ ∑ %)! # (𝑥% − ̅ 𝑥)" ~𝒩 0,1 置き換える 確率分布は？ (𝑠𝑒 * 𝑎# )%= 1 𝑛 − 𝑘 − 1 ∑!"# \$ (𝑦! − - 𝑦)% ∑!"# \$ (𝑦! − * 𝑦)% ~𝝌𝟐(𝒏 − 𝒌 − 𝟏) 5. 回帰係数のt検定
117. ### t分布 互いに独⽴な、 標準正規分布𝒩(0,1)に従う確率変数𝑋と、 カイ⼆乗分布χ"(𝑛)に従う確率変数𝑌について 𝑍 = 𝑋 𝑌/𝑛 なる𝑍の確率分布を𝑡分布と呼び、 𝑍~𝑡(𝑛)

と表す。(𝑛を⾃由度と呼ぶ)
118. ### F 𝑎! − 𝑎! 𝑠𝑒(F 𝑎! ) ~𝑡(𝑛 − 𝑘

− 1) F 𝑎! − 𝑎! σ"/𝑅𝑆𝑆 ~𝒩 0,1 置き換える 𝑎# = 0の時の左辺を計算し、分布𝑡(𝑛 − 𝑘 − 1)に照らせば良い。 t値 (t-value) 5. 回帰係数のt検定
119. ### > a1_se [1] 0.08907516 > t_value [1] 1.737096 a1_se <-

dat_res %>% summarize(RSD = sqrt(sum(res^2) / (N - 2)), vxx = sum((x - mean(x))^2)) %\$% sqrt(RSD^2 / vxx) t_value <- a1 / a1_se 5. 回帰係数のt検定
120. ### > a1_se [1] 0.08907516 > t_value [1] 1.737096 Coefficients: Estimate

Std. Error t value Pr(>|t|) (Intercept) -0.16882 0.55270 -0.305 0.768 x 0.15473 0.08908 1.737 0.121 5. 回帰係数のt検定

122. ### t分布 dat_dt <- tibble(x = seq(-5, 5, by = 0.01),

y = dt(x, df = 8)) geom_path() geom_vline( xintercept = t_value, linetype = "dotted") geom_ribbon( data = dat_dt %>% filter(x >= t_value), aes(ymin = 0, ymax = y), fill = "grey") geom_ribbon( data = dat_dt %>% filter(x <= t_value), aes(ymin = 0, ymax = y), fill = "grey")
123. ### > pt(t_value, 8, lower.tail = F) * 2 [1] 0.1205754

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.16882 0.55270 -0.305 0.768 x 0.15473 0.08908 1.737 0.121 5. 回帰係数のt検定
124. ### > fit_lm %>% summary() Call: lm(formula = y ~ x,

data = .) Residuals: Min 1Q Median 3Q Max -0.9800 -0.6410 0.2338 0.2678 1.5452 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.16882 0.55270 -0.305 0.768 x 0.15473 0.08908 1.737 0.121 Residual standard error: 0.8091 on 8 degrees of freedom Multiple R-squared: 0.2739, Adjusted R-squared: 0.1831 F-statistic: 3.018 on 1 and 8 DF, p-value: 0.1206 deviation 5. 回帰係数のt検定
125. ### F 𝑎+ − 𝑎+ 𝑠𝑒(F 𝑎+ ) ~𝑡(𝑛 − 𝑘

− 1) F 𝑎+ − 𝑎+ σ"の関数!? ~𝒩 0,1 σを𝑅𝑆𝐷に置き換える 5. 回帰係数のt検定
126. ### > 𝑎- = = 𝑦 − > 𝑎+ ̅ 𝑥

= 𝑎- + 𝑎+ ̅ 𝑥 − = 𝑢 − > 𝑎+ ̅ 𝑥 𝐸[> 𝑎- ] = 𝐸[𝑎- + 𝑎+ ̅ 𝑥 − = 𝑢 − > 𝑎+ ̅ 𝑥] = 𝑎- + 𝑎+ ̅ 𝑥 − 0 − 𝑎+ ̅ 𝑥 = 𝑎- 𝑉𝑎𝑟 > 𝑎- = 𝐸 (> 𝑎- −𝐸[> 𝑎- ])( = 𝐸 (> 𝑎- −𝑎- )( = σ( 1 𝑛 + ̅ 𝑥( ∑ )*+ , (𝑥) − ̅ 𝑥)( = … 5. 回帰係数のt検定
127. ### F 𝑎+ − 𝑎+ 𝑠𝑒(F 𝑎+ ) ~𝑡(𝑛 − 𝑘

− 1) F 𝑎+ − 𝑎+ σ" 1 𝑛 + ̅ 𝑥" ∑ %)! # (𝑥% − ̅ 𝑥)" ~𝒩 0,1 σを𝑅𝑆𝐷に置き換える 𝑠𝑒 * 𝑎4 ≔ 𝑅𝑆𝐷 1 𝑛 + ̅ 𝑥% ∑!"# \$ (𝑥! − ̅ 𝑥)% 5. 回帰係数のt検定
128. ### > a0_se [1] 0.5526968 > t_a0 [1] -0.3054542 a0_se <-

dat_res %>% summarize(RSD = sqrt(sum(res^2) / (N - 2)), vxx = sum((x - mean(x))^2), barx2 = mean(x)^2) %\$% {RSD * sqrt(1 / N + barx2 / vxx)} t_a0 <- a0 / a0_se 5. 回帰係数のt検定
129. ### > pt(t_a0, 8, lower.tail = T) * 2 1] 0.767817

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.16882 0.55270 -0.305 0.768 x 0.15473 0.08908 1.737 0.121 5. 回帰係数のt検定
130. ### 回帰を(ほぼ)完全に理解した > fit_lm %>% summary() Call: lm(formula = y ~

x, data = .) Residuals: Min 1Q Median 3Q Max -0.9800 -0.6410 0.2338 0.2678 1.5452 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.16882 0.55270 -0.305 0.768 x 0.15473 0.08908 1.737 0.121 Residual standard error: 0.8091 on 8 degrees of freedom Multiple R-squared: 0.2739, Adjusted R-squared: 0.1831 F-statistic: 3.018 on 1 and 8 DF, p-value: 0.1206 deviation 1 2 4 5 1. 残差 2. 残差標準誤差 3. 決定係数 4. f検定 5. 回帰係数 6. t検定
131. ### Pipe演算⼦: %>% Dollar演算⼦: %\$% 線形回帰: lm() names属性: names() 要約統計量の算出：summarize() 要約統計量の算出：summarize_at()

確率分布関連：pt(), dt(), pf(), df() 技術的復習