Upgrade to Pro — share decks privately, control downloads, hide ads and more …

TokyoR#101_RegressionAnalysis

kilometer
September 17, 2022

 TokyoR#101_RegressionAnalysis

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

kilometer

September 17, 2022
Tweet

More Decks by kilometer

Other Decks in Technology

Transcript

  1. #101 @kilometer00 2022.09.17 Regression Analysis

  2. Who!? 誰だ?

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

    医療システム⼯学 R歴: ~ 10年ぐらい 流⾏: パン切り包丁
  4. 宣伝!! (書籍の翻訳に参加しました。) 第五刷!!

  5. None
  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!!