Upgrade to PRO for Only $50/Year—Limited-Time Offer! 🔥

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. library(tidyverse) library(magrittr) set.seed(1) 準備 N <- 10 # データ数 a

    <- 0.1 # 傾き dat <- tibble(x = c(1:N), y = a * x + rnorm(N))
  2. > 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 準備
  3. 𝑦 = 𝑎! + 𝑎" 𝑥 + 𝑢 𝑢 ~

    !.!.#. 𝒩(0, σ$) (単)回帰モデル
  4. 線形モデル 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
  5. 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
  6. > 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
  7. 例えば、 > 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 この中に誤記があります(知ってます?)。
  8. 今⽇のミッション > 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 この中に書かれた数字を全て計算する。
  9. > 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
  10. 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
  11. 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
  12. 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
  13. > 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
  14. res <- fit_lm$residuals 残差εは ε = 𝑦 − ) 𝑦

    データ 推定値 計算してみましょう。 1. 残差 residuals
  15. 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
  16. 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
  17. dat_res <- dat_fit %>% mutate(y_hat = map( fit, predict)) %>%

    select(-fit) %>% unnest(everything()) %>% mutate(res = y - y_hat) 1. 残差 residuals
  18. > 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
  19. 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
  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 1. 残差 residuals
  21. > 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
  22. 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
  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 線形モデル linear model
  24. 標準誤差 standard error 𝑆𝐸 = 𝑆𝐷 𝑛 𝑆𝐷 = 1

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

    0.8091 # A tibble: 1 × 1 res <dbl> 1 0.241 !!?
  26. 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
  27. 𝑅𝑆𝐷 = 1 𝑛 − 𝑘 − 1 ) !"#

    $ (𝑦! − , 𝑦)% 残差標準偏差 Residual standard deviation 説明変数の次元(今は𝑘 = 1) 誤差𝑢の標準偏差の期待値
  28. > ?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
  29. > 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
  30. 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
  31. 3つの変動 残差変動 RSS 全変動 TSS 回帰変動 ESS 𝑅𝑆𝑆 = &

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

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

    回帰変動 残差変動 3. 決定係数 R-squared 決定係数
  34. 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つの変動
  35. > 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つの変動
  36. 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
  37. F分布 𝑌 = 𝑋! /𝑚 𝑋" /𝑛 互いに独⽴な𝑋! ~𝜒"(𝑚)、𝑋" ~𝜒"(𝑛)について、

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

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

    について、 𝑍 = 1 12! 3 (𝑋1 −𝑢)" σ" なる𝑍は⾃由度𝑘 − 1の𝜒"分布に従う 𝑍~𝜒"(𝑘 − 1)
  40. ⾃由度 degrees of freedom 変数のうち⾃由な値をとりうるものの数 ⺟集団 標本𝑥 (𝑛個) ⾃由に選ばれた𝑛個の標本の⾃由度は𝒏で、 標本の平均値は標本平均

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

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

    !.!.#. 𝒩(0, σ$) (単)回帰モデル
  43. 𝑦 = 𝑎! + 𝑎" 𝑥 + 𝑢 𝑢 ~

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

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

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

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

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

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

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

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

    𝑍 = 1 12! 3 (𝑋1 − 8 𝑋)" σ" なる𝑍は⾃由度𝑘 − 1の𝜒"分布に従う 𝑍~𝜒"(𝑘 − 1)
  52. 回帰変動は(理想的には) ⾃由度𝑘のχ2分布に従う。 𝐸𝑆𝑆 = 𝑇𝑆𝑆 − 𝑅𝑆𝑆 𝜒!(𝑛 − 1)

    𝜒!(𝑛 − 𝑘 − 1) 𝜒!(𝑘) これらが互いに独⽴なので、 & !"# $ * 𝑦! − - 𝑦 = 0
  53. 回帰モデルの適合性 𝑅2 = 𝐸𝑆𝑆 𝑇𝑆𝑆 = 1 − 𝑅𝑆𝑆 𝑇𝑆𝑆

    𝑅2 = 1 − 𝑅𝑆𝑆/(𝑛 − 𝑘 − 1) 𝑇𝑆𝑆/(𝑛 − 1) 決定係数 ⾃由度調整済み決定係数
  54. > 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
  55. 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
  56. F分布 𝑌 = 𝑋! /𝑚 𝑋" /𝑛 互いに独⽴な𝑋! ~𝜒"(𝑚)、𝑋" ~𝜒"(𝑛)について、

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

    分散0なら0 𝑓~𝐹(𝑘, 𝑛 − 𝑘 − 1) 4. ⾃由度とf検定
  58. > 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
  59. > 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値
  60. 回帰モデルの適合性 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
  61. 回帰 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"
  62. 回帰 dat %>% lm(y ~ x, data = .) %>%

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

    summary() %$% fstatistic %>% { pf(.[1], .[2], .[3], lower.tail = FALSE) } value 0.1205754
  64. 回帰 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
  65. > 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検定
  66. > 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. 回帰係数
  67. 最⼩⼆乗法(線形回帰) argmin(-!,-") ; %)! # ε" 𝑦 = / 𝑎&

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

    − ̅ 𝑥)( 𝐶𝑜𝑣[𝑥, 𝑦] = 1 𝑛 − 1 6 )*+ , (𝑥) − ̅ 𝑥)(𝑦) − = 𝑦) 最⼩⼆乗法(線形回帰) < 𝑎! = ∑$%! & (𝑥$ − ̅ 𝑥)(𝑦$ − A 𝑦) ∑ $%! & (𝑥$ − ̅ 𝑥)" = 𝐶𝑜𝑣[𝑥, 𝑦] 𝑉𝑎𝑟[𝑥] < 𝑎' = A 𝑦 − < 𝑎! ̅ 𝑥 = A 𝑦 − 𝐶𝑜𝑣[𝑥, 𝑦] 𝑉𝑎𝑟[𝑥] ̅ 𝑥 ただし、 分散 共分散
  69. 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 )*+ , (𝑥) − ̅ 𝑥)(𝑦) − = 𝑦) 共分散
  70. 𝑉𝑎𝑟[𝑥] = 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)) 最⼩⼆乗法(線形回帰) 分散 共分散
  71. # 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 & !"# $ (𝑥! − ̅ 𝑥)(𝑦! − - 𝑦) 分散 共分散
  72. 最⼩⼆乗法(線形回帰) > 𝑎+ = 𝐶𝑜𝑣[𝑥, 𝑦] 𝑉𝑎𝑟[𝑥] , > 𝑎-

    = = 𝑦 − 𝐶𝑜𝑣[𝑥, 𝑦] 𝑉𝑎𝑟 𝑥 ̅ 𝑥 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
  73. 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
  74. 最⼩⼆乗法(線形回帰) 𝑦 = 2 𝑎? + 2 𝑎@ 𝑥 +

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

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

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

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

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

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

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

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

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

    𝑎! − 𝑎! σ"/ ∑ %)! # (𝑥% − ̅ 𝑥)" ~𝒩 0,1 置き換える 確率分布は? (𝑠𝑒 * 𝑎# )%= 1 𝑛 − 𝑘 − 1 ∑!"# $ (𝑦! − - 𝑦)% ∑!"# $ (𝑦! − * 𝑦)% ~𝝌𝟐(𝒏 − 𝒌 − 𝟏) 5. 回帰係数のt検定
  84. F 𝑎! − 𝑎! 𝑠𝑒(F 𝑎! ) ~𝑡(𝑛 − 𝑘

    − 1) F 𝑎! − 𝑎! σ"/𝑅𝑆𝑆 ~𝒩 0,1 置き換える 𝑎# = 0の時の左辺を計算し、分布𝑡(𝑛 − 𝑘 − 1)に照らせば良い。 t値 (t-value) 5. 回帰係数のt検定
  85. > 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検定
  86. > 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検定
  87. 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")
  88. > 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検定
  89. > 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検定
  90. F 𝑎+ − 𝑎+ 𝑠𝑒(F 𝑎+ ) ~𝑡(𝑛 − 𝑘

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

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

    − 1) F 𝑎+ − 𝑎+ σ" 1 𝑛 + ̅ 𝑥" ∑ %)! # (𝑥% − ̅ 𝑥)" ~𝒩 0,1 σを𝑅𝑆𝐷に置き換える 𝑠𝑒 * 𝑎4 ≔ 𝑅𝑆𝐷 1 𝑛 + ̅ 𝑥% ∑!"# $ (𝑥! − ̅ 𝑥)% 5. 回帰係数のt検定
  93. > 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検定
  94. > 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検定
  95. 回帰を(ほぼ)完全に理解した > 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検定