TokyoR#101_RegressionAnalysis

第101回Tokyo.Rでトークした際のスライドです。

Transcript

1. #101 @kilometer00 2022.09.17 Regression Analysis
2. Who！？ 誰だ？
3. Who！？ 名前： 三村 @kilometer 職業： ポスドク (こうがくはくし) 専⾨： ⾏動神経科学(霊⻑類) 脳イメージング 医療システム⼯学 R歴： ~ 10年ぐらい 流⾏: パン切り包丁
4. 宣伝!! （書籍の翻訳に参加しました。） 第五刷！！
6. ⼀度は学んでおきたい 線形回帰分析
7. 今⽇の⽬標 「回帰分析？ ああ、知ってるよ。」 と⾔えるようになる。
8. 数式を使わずに直感的に理解できる説明をします！ 数式を(たくさん)使います。 背景原理はここでは説明しません。 背景原理も(出来るだけ)説明します。 サルでもわかる！タコでも！ネコでも！！ (認知神経科学ガチ勢としてはちょっとコレは保証できかねます)
9. 前提 ・Rがある程度使えること (応⽤セッションなので%$%とか普通に使います) ・確率, 確率分布に関する基礎知識 (#96 Tokyo.Rのトークを参照) ・回帰分析を実⾏した経験が 全くないと多分、分からない
10. 回帰していますか？ みなさん、
11. あ、今⽇はまだ 回帰していません！ えっ！ 今すぐ回帰しましょう。
12. library(tidyverse) library(magrittr) set.seed(1) 準備 N <- 10 # データ数 a <- 0.1 # 傾き dat <- tibble(x = c(1:N), y = a * x + rnorm(N))
13. > 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 準備
14. 𝑦 = 𝑎! + 𝑎" 𝑥 + 𝑢 𝑢 ~ !.!.#. 𝒩(0, σ$) (単)回帰モデル
15. 線形モデル 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
16. 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
17. > 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
18. かいき？ ああ、SNSで話題だよね。 知ってる知ってる。
19. 結果のsummaryに表⽰される 全ての数字を みたことがある 使ったことがある 正確な意味を知っている 値を計算で求めることができる 他⼈に説明したことがある 説明できる資料がある
20. 結果のsummaryに表⽰される 全ての数字を みたことがある 使ったことがある 正確な意味を知っている 値を計算で求めることができる 他⼈に説明したことがある 説明できる資料がある(壇上にどうぞ)
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 この中に書かれた数字を全て計算する。
23. > 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
24. 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
25. 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
26. 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
27. > 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
28. res <- fit_lm$residuals 残差εは ε = 𝑦 − ) 𝑦 データ 推定値 計算してみましょう。 1. 残差 residuals
29. 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. コレでいいんですが、 こうやって回帰係数を引っ張るのは 説明変数が増えたり減ったりすると とてつもなく⾯倒臭くなるので、 出来るだけ格好よくやりたい。 fit_lm$coefficients[1] 1. 残差 residuals
31. 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
32. dat_res <- dat_fit %>% mutate(y_hat = map( fit, predict)) %>% select(-fit) %>% unnest(everything()) %>% mutate(res = y - y_hat) 1. 残差 residuals
33. > 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
34. 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
35. > 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
36. > 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
37. 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
38. > 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
39. 標準誤差 standard error 𝑆𝐸 = 𝑆𝐷 𝑛 𝑆𝐷 = 1 𝑛 − 1 ( !"# $ (𝑦! − + 𝑦)% 標本平均の標準偏差 （の推定値）
40. Residual standard error: 0.8091 dat_res %>% summarize_at("res", ~ sd(.) / sqrt(N)) # A tibble: 1 × 1 res <dbl> 1 0.241
41. dat_res %>% summarize_at("res", ~ sd(.) / sqrt(N)) Residual standard error: 0.8091 # A tibble: 1 × 1 res <dbl> 1 0.241 !!?
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
43. 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
44. 𝑅𝑆𝐷 = 1 𝑛 − 𝑘 − 1 ) !"# $ (𝑦! − , 𝑦)% 残差標準偏差 Residual standard deviation 説明変数の次元（今は𝑘 = 1） 誤差𝑢の標準偏差の期待値
45. > ?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
46. > 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
47. 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
48. 3つの変動 残差変動 RSS 全変動 TSS 回帰変動 ESS 𝑅𝑆𝑆 = & !"# $ (𝑦! − * 𝑦)% 𝑇𝑆𝑆 = & !"# $ (𝑦! − - 𝑦)% 𝐸𝑆𝑆 = & !"# $ (* 𝑦 − - 𝑦)% データ 推定値 データ 平均 推定値 平均
49. 3つの変動 残差変動 RSS 全変動 TSS 回帰変動 ESS 𝑅𝑆𝑆 = & !"# $ (𝑦! − * 𝑦)% 𝑇𝑆𝑆 = & !"# $ (𝑦! − - 𝑦)% 𝐸𝑆𝑆 = & !"# $ (* 𝑦 − - 𝑦)% データ 推定値 データ 平均 推定値 平均 = -
50. 𝑅2 = 𝐸𝑆𝑆 𝑇𝑆𝑆 = 1 − 𝑅𝑆𝑆 𝑇𝑆𝑆 全変動 回帰変動 残差変動 3. 決定係数 R-squared 決定係数
51. 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つの変動
52. > 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つの変動
53. 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. ぴーち？ ああ、あれは美味しいよね。 夏のデザートにぴったり。 ん？えふち？ 知ってる知ってる。 あれだよね、あれ、はは
55. F分布 𝑌 = 𝑋! /𝑚 𝑋" /𝑛 互いに独⽴な𝑋! ~𝜒"(𝑚)、𝑋" ~𝜒"(𝑛)について、 なる𝑌の確率分布をF分布と呼び、 𝑌~𝐹(𝑚, 𝑛) と表す。 (𝑚を第1⾃由度,𝑛を第2⾃由度と呼ぶ) カイ2乗分布
56. 𝜒!分布 標準正規分布𝒩(0,1)に従う独⽴な𝑘個の 確率変数𝑋! , 𝑋" , … , 𝑋# について、 𝑍 = 1 12! 3 𝑋1 " なる 𝑍 の確率分布を 𝜒" 分布と呼び、 𝑍~𝜒"(𝑘) と表す。(𝑘を⾃由度と呼ぶ)
57. 𝜒!分布 正規分布𝒩(𝑢, σ")に従う独⽴な𝑘個の 確率変数𝑋! , 𝑋" , … , 𝑋# について、 𝑍 = 1 12! 3 (𝑋1 −𝑢)" σ" なる𝑍は⾃由度𝑘 − 1の𝜒"分布に従う 𝑍~𝜒"(𝑘 − 1)
58. 正規分布に従う独⽴な𝑘個の確率変数 𝑋! , 𝑋" , … , 𝑋# の2乗和、すなわちコレ！？ データ数𝑘の時、残差⼆乗和は ⾃由度𝑘 − 1のχ2分布に従う？
59. 正規分布に従う独⽴な𝑘個の確率変数 𝑋! , 𝑋" , … , 𝑋# の2乗和、すなわちコレ！？ データ数𝑘の時、残差⼆乗和は ⾃由度𝑘 − 1のχ2分布に従う？ → 従わない
60. 正規分布に従う独⽴な𝑘個の確率変数 𝑋! , 𝑋" , … , 𝑋# の2乗和、すなわちコレ！？ コレは⾃由度𝑘 − 1のχ2分布に従う。
61. ⾃由度 degrees of freedom 変数のうち⾃由な値をとりうるものの数
62. ⾃由度 degrees of freedom 変数のうち⾃由な値をとりうるものの数 ⺟集団 標本𝑥 (𝑛個) ⾃由に選ばれた𝑛個の標本の⾃由度は𝒏で、 標本の平均値は標本平均 ̅ 𝑥。 では、標本平均が既知の値 ̅ 𝑥になるように 選ばれた𝑛個の標本の⾃由度は？
63. ⾃由度 degrees of freedom 変数のうち⾃由な値をとりうるものの数 ⺟集団 標本𝑥 (𝑛個) ⾃由に選ばれた𝑛個の標本の⾃由度は𝒏で、 標本の平均値は標本平均 ̅ 𝑥。 では、標本平均が既知の値 ̅ 𝑥になるように 選ばれた𝑛個の標本の⾃由度は？ → 𝒏 − 𝟏
64. ⾃由度 degrees of freedom 変数のうち⾃由な値をとりうるものの数 標本平均が既知の値 ̅ 𝑥になるように 選ばれた𝑛個の標本の⾃由度は 𝒏 − 𝟏 ̅ 𝑥 = 1 𝑛 (𝑥! + 𝑥" + ⋯ + 𝑥#$! + 𝑥# ) 𝑥# = 𝑛 ̅ 𝑥 − (𝑥! + 𝑥" + ⋯ + 𝑥#$! ) ⾃由に値をとる𝒏 − 𝟏項 既知の 標本平均
65. ここまでおさえておいて、
66. 回帰モデルの考え⽅を再確認。
67. 𝑦 = 𝑎! + 𝑎" 𝑥 + 𝑢 𝑢 ~ !.!.#. 𝒩(0, σ$) (単)回帰モデル
68. データ ⺟集団の(真の) 回帰直線 誤差 𝑢 𝑥 𝑦 𝑌
69. データ データから 推定された 回帰直線 ⺟集団の(真の) 回帰直線 誤差 残差 & 𝑦 ε 𝑢 𝑦 𝑥 𝑌
70. 𝑦 = 𝑎! + 𝑎" 𝑥 + 𝑢 𝑢 ~ !.!.#. 𝒩(0, σ$) 誤差 𝑦 = ) 𝑎! + ) 𝑎" 𝑥 + ε ε = 𝑦 − ) 𝑦 真の回帰モデル 残差 最⼩⼆乗法で推定された回帰式
71. (正規)回帰モデル5つの仮定 1. 誤差の期待値はゼロ 2. 誤差の分散は定数 3. 誤差は互いに無相関 4. 誤差と説明変数は無相関 5. 誤差は正規分布に従う 𝑢 ~ %.%.'. 𝒩(0, σ") 𝐸[𝑢% (𝑥% − 𝐸[𝑥% ])] = 0 𝐸[𝑢% 𝑢( ] = 0 (𝑖 ≠ 𝑗) 𝑉𝑎𝑟 𝑢 = σ" 𝐸 𝑢 = 0
72. (正規)回帰モデル5つの仮定 1. 誤差の期待値はゼロ 2. 誤差の分散は定数 3. 誤差は互いに無相関 4. 誤差と説明変数は無相関 ガウス=マルコフの定理 仮定1-4が満たされる場合、最⼩⼆乗推定量は 線型不偏推定量の中で最⼩分散を持つ。 5. 誤差は正規分布に従う + 仮定5 最⼩⼆乗推定量は最⼩分散を持つ。
73. 𝑦 = 𝑎! + 𝑎" 𝑥 + 𝑢 𝑢 ~ !.!.#. 𝒩(0, σ$) 誤差 真の回帰モデル 最⼩⼆乗法OLSでモデルの パラメータを推定する (OLS, Ordinary Last Squares) → 残差⼆乗和RSSの最⼩化 (RSS, Residual Sum of Squares)
74. OLSが要請する残差の性質 1. 残差の合計値はゼロ 2. 残差と説明変数は無相関 ; %)! # ε% = 0 ; %)! # 𝑥% ε% = 0
75. 1. 残差の合計値はゼロ 2. 残差と説明変数は無相関 ; %)! # ε% = 0 ; %)! # 𝑥*% ε% = 0 (𝑘 = 1, 2, … ) 説明変数の次元数 OLSが要請する残差の性質
76. 𝑦 = / 𝑎& + / 𝑎# 𝑥# + ⋯ + / 𝑎' 𝑥' + ε 残差 OLSで推定された(重)回帰式 ; %)! # ε% = 0 ; %)! # 𝑥*% ε% = 0 (𝑘 = 1, 2, … ) かつ 𝑛個の残差εは、 𝑘 + 1個の式を満たさなければならない。 残差εの⾃由度はn − 𝑘 − 1となる。 𝑘項
77. 残差εの⾃由度はn − 𝑘 − 1となる。
78. 𝑅𝑆𝐷 = 1 𝑛 − 𝑘 − 1 ) !"# $ (𝑦! − , 𝑦)% 残差標準偏差 Residual standard deviation 説明変数の次元（今は𝑘 = 1） 𝑦 = 𝑎+ + 𝑎! 𝑥 + 𝑢
79. 1. 残差の合計値はゼロ 2. 残差と説明変数は無相関 ; %)! # ε% = 0 ; %)! # 𝑥*% ε% = 0 (𝑘 = 1, 2, … ) 説明変数の次元数 OLSが要請する残差の性質 ※ 定数項/ 𝒂𝟎 が無い回帰式では この性質は満たされない。 → 残差𝜺の⾃由度は𝐧 − 𝒌。
80. 𝜒!分布 正規分布に従う独⽴な𝑘個の確率変数 𝑋! , 𝑋" , … , 𝑋# について、 𝑍 = 1 12! 3 (𝑋1 − 8 𝑋)" σ" なる𝑍は⾃由度𝑘 − 1の𝜒"分布に従う 𝑍~𝜒"(𝑘 − 1)
81. データ数𝑛の時、残差平⽅和は(理想的には) ⾃由度𝑛 − 𝑘 − 1のχ2分布に従う。 𝑅𝑆𝑆~𝜒"(𝑛 − 𝑘 − 1)
82. データ数𝑛の時、全変動は(理想的には) ⾃由度𝑛 − 1のχ2分布に従う。 𝑇𝑆𝑆~𝜒"(𝑛 − 1) & !"# $ 𝑦! − - 𝑦 = 0
83. 回帰変動は(理想的には) ⾃由度𝑘のχ2分布に従う。 𝐸𝑆𝑆 = 𝑇𝑆𝑆 − 𝑅𝑆𝑆 𝜒!(𝑛 − 1) 𝜒!(𝑛 − 𝑘 − 1) 𝜒!(𝑘) これらが互いに独⽴なので、 & !"# $ * 𝑦! − - 𝑦 = 0
84. 𝑅𝑆𝑆~𝜒((𝑛 − 𝑘 − 1) 𝑇𝑆𝑆~𝜒((𝑛 − 1) 𝐸𝑆𝑆~𝜒((𝑘) OLSの変動の確率分布 残差変動 全変動 回帰変動
85. 回帰モデルの適合性 𝑅2 = 𝐸𝑆𝑆 𝑇𝑆𝑆 = 1 − 𝑅𝑆𝑆 𝑇𝑆𝑆 𝑅2 = 1 − 𝑅𝑆𝑆/(𝑛 − 𝑘 − 1) 𝑇𝑆𝑆/(𝑛 − 1) 決定係数 ⾃由度調整済み決定係数
86. > 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
87. 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
88. F分布 𝑌 = 𝑋! /𝑚 𝑋" /𝑛 互いに独⽴な𝑋! ~𝜒"(𝑚)、𝑋" ~𝜒"(𝑛)について、 なる𝑌の確率分布をF分布と呼び、 𝑌~𝐹(𝑚, 𝑛) と表す。 (𝑚を第1⾃由度,𝑛を第2⾃由度と呼ぶ) カイ2乗分布
89. 𝑓 = 𝐸𝑆𝑆/𝑘 𝑅𝑆𝑆/(𝑛 − 𝑘 − 1) F値 無相関(傾き=0)なら0 分散0なら0 𝑓~𝐹(𝑘, 𝑛 − 𝑘 − 1) 4. ⾃由度とf検定
90. > 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
91. > 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値
92. 回帰モデルの適合性 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
93. 回帰 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"
94. 回帰 dat %>% lm(y ~ x, data = .) %>% summary() %$% fstatistic value numdf dendf 3.017501 1.000000 8.000000
95. 回帰 dat %>% lm(y ~ x, data = .) %>% summary() %$% fstatistic %>% { pf(.[1], .[2], .[3], lower.tail = FALSE) } value 0.1205754
96. 回帰 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
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 4. ⾃由度とf検定
98. > 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. 回帰係数
99. 最⼩⼆乗法（線形回帰） argmin(-!,-") ; %)! # ε" 𝑦 = / 𝑎& + / 𝑎# 𝑥 + ε 残差 & 𝑎" = ∑#$" % (𝑥# − ̅ 𝑥)(𝑦# − / 𝑦) ∑ #$" % (𝑥# − ̅ 𝑥)! & 𝑎& = / 𝑦 − & 𝑎" ̅ 𝑥
100. 𝑉𝑎𝑟[𝑥] = 1 𝑛 − 1 6 )*+ , (𝑥) − ̅ 𝑥)( 𝐶𝑜𝑣[𝑥, 𝑦] = 1 𝑛 − 1 6 )*+ , (𝑥) − ̅ 𝑥)(𝑦) − = 𝑦) 最⼩⼆乗法（線形回帰） < 𝑎! = ∑$%! & (𝑥$ − ̅ 𝑥)(𝑦$ − A 𝑦) ∑ $%! & (𝑥$ − ̅ 𝑥)" = 𝐶𝑜𝑣[𝑥, 𝑦] 𝑉𝑎𝑟[𝑥] < 𝑎' = A 𝑦 − < 𝑎! ̅ 𝑥 = A 𝑦 − 𝐶𝑜𝑣[𝑥, 𝑦] 𝑉𝑎𝑟[𝑥] ̅ 𝑥 ただし、 分散 共分散
101. 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 )*+ , (𝑥) − ̅ 𝑥)(𝑦) − = 𝑦) 共分散
102. 最⼩⼆乗法（線形回帰） 𝐶𝑜𝑣[𝑥, 𝑦] = 1 𝑛 − 1 6 )*+ , (𝑥) − ̅ 𝑥)(𝑦) − = 𝑦) 共分散
103. 𝑉𝑎𝑟[𝑥] = 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)) 最⼩⼆乗法（線形回帰） 分散 共分散
104. # 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 & !"# $ (𝑥! − ̅ 𝑥)(𝑦! − - 𝑦) 分散 共分散
105. 最⼩⼆乗法（線形回帰） > 𝑎+ = 𝐶𝑜𝑣[𝑥, 𝑦] 𝑉𝑎𝑟[𝑥] , > 𝑎- = = 𝑦 − 𝐶𝑜𝑣[𝑥, 𝑦] 𝑉𝑎𝑟 𝑥 ̅ 𝑥 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
106. 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. ぴーち？ ああ、あれは...
108. 最⼩⼆乗法（線形回帰） 𝑦 = 2 𝑎? + 2 𝑎@ 𝑥 + ε 回帰係数< 𝑎' と< 𝑎! の仮説検定を⾏う。 帰無仮説：< 𝑎! = 0 対⽴仮説：< 𝑎! ≠ 0 すなわち例えば< 𝑎! について、 → < 𝑎' と< 𝑎! の確率分布を知りたい。 → < 𝑎' と< 𝑎! の期待値と分散を調べる。
109. 𝑦 = 𝑎! + 𝑎" 𝑥 + 𝑢 𝑢 ~ !.!.#. 𝒩(0, σ$) 誤差 𝑦 = ) 𝑎! + ) 𝑎" 𝑥 + ε ε = 𝑦 − ) 𝑦 真の回帰モデル 残差 最⼩⼆乗法で推定された回帰式
110. 期待と分散の基本的性質 𝐸 𝑋 + 𝑌 = 𝐸 𝑋 + 𝐸[𝑌] 𝐸 𝑐𝑋 = 𝑐𝐸 𝑋 𝑋と𝑌が独⽴の時、𝐸 𝑋𝑌 = 𝐸 𝑋 𝐸 𝑌 𝑉𝑎𝑟 𝑋 = 𝐸 (𝑋 − 𝐸[𝑋])" = 𝐸 𝑋" − (𝐸[𝑋])" 𝑉𝑎𝑟 𝑋 + 𝑌 = 𝑉𝑎𝑟 𝑋 + 𝑉𝑎𝑟 𝑌 + 2𝐶𝑜𝑣[𝑋, 𝑌] 𝑉𝑎𝑟 𝑐𝑋 = 𝑐"𝑉𝑎𝑟 𝑋 𝑉𝑎𝑟 𝑋 + 𝑐 = 𝑉𝑎𝑟 𝑋
111. > 𝑎+ = ∑)*+ , (𝑥) − ̅ 𝑥)(𝑦) − = 𝑦) ∑ )*+ , (𝑥) − ̅ 𝑥)( = ∑)*+ , (𝑥) − ̅ 𝑥)( 𝑎- + 𝑎+ 𝑥) + 𝑢) − 𝑎- + 𝑎+ ̅ 𝑥 + = 𝑢 ) ∑ )*+ , (𝑥) − ̅ 𝑥)( = ∑)*+ , 𝑥) − ̅ 𝑥 (𝑎+ 𝑥) − ̅ 𝑥 + 𝑢) − = 𝑢) ∑ )*+ , (𝑥) − ̅ 𝑥)( = 𝑎+ + ∑)*+ , 𝑥) − ̅ 𝑥 𝑢) ∑ )*+ , (𝑥) − ̅ 𝑥)( + ∑)*+ , 𝑥) − ̅ 𝑥 ∑ )*+ , (𝑥) − ̅ 𝑥)( = 𝑢 𝐸[F 𝑎! ] = 𝑎! 従って、 …式(1) 5. 回帰係数のt検定
112. 𝑉𝑎𝑟 < 𝑎! = 𝐸 (< 𝑎! − 𝐸[< 𝑎! ])" = 𝐸 (< 𝑎! − 𝑎! )" = 𝐸 ∑!"# $ 𝑥! − ̅ 𝑥 𝑢! ∑ !"# $ (𝑥! − ̅ 𝑥)% % = 𝐸 ∑!"# $ 𝑥! − ̅ 𝑥 %𝑢! % ∑ !"# $ (𝑥! − ̅ 𝑥)% % + 𝐸 ∑!"# $ ∑&"# $ 𝑥! − ̅ 𝑥 𝑥& − ̅ 𝑥 𝑢! 𝑢& ∑ !"# $ (𝑥! − ̅ 𝑥)% % 𝑖 ≠ 𝑗 = ∑!"# $ 𝑥! − ̅ 𝑥 %𝐸 𝑢! % ∑ !"# $ (𝑥! − ̅ 𝑥)% % + ∑!"# $ ∑&"# $ 𝑥! − ̅ 𝑥 𝑥& − ̅ 𝑥 𝐸[𝑢! 𝑢& ] ∑ !"# $ (𝑥! − ̅ 𝑥)% % = 𝐸 𝑢% − (𝐸[𝑢])% ∑ !"# $ (𝑥! − ̅ 𝑥)% = 𝑉𝑎𝑟[𝑢%] ∑ !"# $ (𝑥! − ̅ 𝑥)% = σ% ∑ !"# $ (𝑥! − ̅ 𝑥)% 𝐸 𝑢 = 0 𝑉𝑎𝑟 𝑎! = 𝐸 𝑎! − (𝐸 𝑎 )! 5. 回帰係数のt検定
113. 𝐸 < 𝑎! = 𝑎! 𝑉𝑎𝑟[< 𝑎! ] = σ" ∑ $%! & (𝑥$ − ̅ 𝑥)" 頑張って計算した結果、 > 𝑎+ = 𝑎+ + ∑)*+ , 𝑥) − ̅ 𝑥 ∑ )*+ , (𝑥) − ̅ 𝑥)( = 𝑢 + ∑)*+ , 𝑥) − ̅ 𝑥 𝑢) ∑ )*+ , (𝑥) − ̅ 𝑥)( また式(1) 𝑢~ 𝒩(0, σ() より * 𝑎# は正規分布に従い、 > 𝑎+ ~ 𝒩 𝑎+ , σ( ∑ )*+ , (𝑥) − ̅ 𝑥)( ⇔ > 𝑎+ − 𝑎+ σ(/ ∑ )*+ , (𝑥) − ̅ 𝑥)( ~𝒩 0,1 5. 回帰係数のt検定
114. F 𝑎! − 𝑎! σ"/ ∑ %)! # (𝑥% − ̅ 𝑥)" ~𝒩 0,1 コレに関する帰無仮説を検証する コレをどうにかしたい(データから推定したい) σは誤差𝑢の標準偏差 5. 回帰係数のt検定
115. 𝑅𝑆𝐷 = 1 𝑛 − 𝑘 − 1 ) !"# $ (𝑦! − , 𝑦)% 残差標準偏差 Residual standard deviation 誤差𝑢の標準偏差の期待値
116. F 𝑎! − 𝑎! 𝑅𝑆𝐷"/ ∑ %)! # (𝑥% − ̅ 𝑥)" = F 𝑎! − 𝑎! 𝑠𝑒(F 𝑎! ) ~ ? F 𝑎! − 𝑎! σ"/ ∑ %)! # (𝑥% − ̅ 𝑥)" ~𝒩 0,1 置き換える 確率分布は？ 𝑠𝑒 > 𝑎+ ∶= 𝑅𝑆𝐷( ∑ )*+ , (𝑥) − ̅ 𝑥)( 5. 回帰係数のt検定
117. F 𝑎! − 𝑎! 𝑠𝑒(F 𝑎! ) ~ ? F 𝑎! − 𝑎! σ"/ ∑ %)! # (𝑥% − ̅ 𝑥)" ~𝒩 0,1 置き換える 確率分布は？ (𝑠𝑒 * 𝑎# )%= 1 𝑛 − 𝑘 − 1 ∑!"# $ (𝑦! − - 𝑦)% ∑!"# $ (𝑦! − * 𝑦)% ~𝝌𝟐(𝒏 − 𝒌 − 𝟏) 5. 回帰係数のt検定
118. t分布 互いに独⽴な、 標準正規分布𝒩(0,1)に従う確率変数𝑋と、 カイ⼆乗分布χ"(𝑛)に従う確率変数𝑌について 𝑍 = 𝑋 𝑌/𝑛 なる𝑍の確率分布を𝑡分布と呼び、 𝑍~𝑡(𝑛) と表す。(𝑛を⾃由度と呼ぶ)
119. F 𝑎! − 𝑎! 𝑠𝑒(F 𝑎! ) ~𝑡(𝑛 − 𝑘 − 1) F 𝑎! − 𝑎! σ"/𝑅𝑆𝑆 ~𝒩 0,1 置き換える 𝑎# = 0の時の左辺を計算し、分布𝑡(𝑛 − 𝑘 − 1)に照らせば良い。 t値 (t-value) 5. 回帰係数のt検定
120. > 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検定
121. > 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分布
123. 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")
124. > 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検定
125. > 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検定
126. F 𝑎+ − 𝑎+ 𝑠𝑒(F 𝑎+ ) ~𝑡(𝑛 − 𝑘 − 1) F 𝑎+ − 𝑎+ σ"の関数!? ~𝒩 0,1 σを𝑅𝑆𝐷に置き換える 5. 回帰係数のt検定
127. > 𝑎- = = 𝑦 − > 𝑎+ ̅ 𝑥 = 𝑎- + 𝑎+ ̅ 𝑥 − = 𝑢 − > 𝑎+ ̅ 𝑥 𝐸[> 𝑎- ] = 𝐸[𝑎- + 𝑎+ ̅ 𝑥 − = 𝑢 − > 𝑎+ ̅ 𝑥] = 𝑎- + 𝑎+ ̅ 𝑥 − 0 − 𝑎+ ̅ 𝑥 = 𝑎- 𝑉𝑎𝑟 > 𝑎- = 𝐸 (> 𝑎- −𝐸[> 𝑎- ])( = 𝐸 (> 𝑎- −𝑎- )( = σ( 1 𝑛 + ̅ 𝑥( ∑ )*+ , (𝑥) − ̅ 𝑥)( = … 5. 回帰係数のt検定
128. F 𝑎+ − 𝑎+ 𝑠𝑒(F 𝑎+ ) ~𝑡(𝑛 − 𝑘 − 1) F 𝑎+ − 𝑎+ σ" 1 𝑛 + ̅ 𝑥" ∑ %)! # (𝑥% − ̅ 𝑥)" ~𝒩 0,1 σを𝑅𝑆𝐷に置き換える 𝑠𝑒 * 𝑎4 ≔ 𝑅𝑆𝐷 1 𝑛 + ̅ 𝑥% ∑!"# $ (𝑥! − ̅ 𝑥)% 5. 回帰係数のt検定
129. > 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検定
130. > 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検定
131. 回帰を(ほぼ)完全に理解した > 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検定
132. Pipe演算⼦: %>% Dollar演算⼦: %$% 線形回帰: lm() names属性: names() 要約統計量の算出：summarize() 要約統計量の算出：summarize_at() 確率分布関連：pt(), dt(), pf(), df() 技術的復習
133. Enjoy!!