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

    View Slide

  2. Who!?
    誰だ?

    View Slide

  3. Who!?
    名前: 三村 @kilometer
    職業: ポスドク (こうがくはくし)
    専⾨: ⾏動神経科学(霊⻑類)
    脳イメージング
    医療システム⼯学
    R歴: ~ 10年ぐらい
    流⾏: パン切り包丁

    View Slide

  4. 宣伝!!
    (書籍の翻訳に参加しました。)
    第五刷!!

    View Slide

  5. View Slide

  6. ⼀度は学んでおきたい
    線形回帰分析

    View Slide

  7. 今⽇の⽬標
    「回帰分析?
    ああ、知ってるよ。」
    と⾔えるようになる。

    View Slide

  8. 数式を使わずに直感的に理解できる説明をします!
    数式を(たくさん)使います。
    背景原理はここでは説明しません。
    背景原理も(出来るだけ)説明します。
    サルでもわかる!タコでも!ネコでも!!
    (認知神経科学ガチ勢としてはちょっとコレは保証できかねます)

    View Slide

  9. 前提
    ・Rがある程度使えること
    (応⽤セッションなので%$%とか普通に使います)
    ・確率, 確率分布に関する基礎知識
    (#96 Tokyo.Rのトークを参照)
    ・回帰分析を実⾏した経験が
    全くないと多分、分からない

    View Slide

  10. 回帰していますか?
    みなさん、

    View Slide

  11. あ、今⽇はまだ
    回帰していません!
    えっ!
    今すぐ回帰しましょう。

    View Slide

  12. library(tidyverse)
    library(magrittr)
    set.seed(1)
    準備
    N a dat tibble(x = c(1:N),
    y = a * x + rnorm(N))

    View Slide

  13. > dat
    # A tibble: 10 × 2
    x y

    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
    準備

    View Slide

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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  18. かいき?
    ああ、SNSで話題だよね。
    知ってる知ってる。

    View Slide

  19. 結果のsummaryに表⽰される
    全ての数字を
    みたことがある
    使ったことがある
    正確な意味を知っている
    値を計算で求めることができる
    他⼈に説明したことがある
    説明できる資料がある

    View Slide

  20. 結果のsummaryに表⽰される
    全ての数字を
    みたことがある
    使ったことがある
    正確な意味を知っている
    値を計算で求めることができる
    他⼈に説明したことがある
    説明できる資料がある(壇上にどうぞ)

    View Slide

  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
    この中に誤記があります(知ってます?)。

    View Slide

  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
    この中に書かれた数字を全て計算する。

    View Slide

  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

    View Slide

  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

    View Slide

  25. res > 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

    View Slide

  26. res > 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

    View Slide

  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

    View Slide

  28. res 残差εは
    ε = 𝑦 − )
    𝑦
    データ 推定値
    計算してみましょう。
    1. 残差 residuals

    View Slide

  29. a0 a1 dat %>%
    mutate(y_hat = a0 + a1 * x,
    res = y - y_hat)
    # A tibble: 10 × 4
    x y y_hat res

    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

    View Slide

  30. コレでいいんですが、
    こうやって回帰係数を引っ張るのは
    説明変数が増えたり減ったりすると
    とてつもなく⾯倒臭くなるので、
    出来るだけ格好よくやりたい。
    fit_lm$coefficients[1]
    1. 残差 residuals

    View Slide

  31. dat_fit dat %>%
    group_nest() %>%
    mutate(fit = map(data,
    ~ lm(y ~ x, data = .)))
    > dat_fit
    # A tibble: 1 × 2
    data fit

    1
    「線形回帰?nestでしょう。」
    1. 残差 residuals

    View Slide

  32. dat_res dat_fit %>%
    mutate(y_hat = map(
    fit, predict)) %>%
    select(-fit) %>%
    unnest(everything()) %>%
    mutate(res = y - y_hat)
    1. 残差 residuals

    View Slide

  33. > dat_res
    # A tibble: 10 × 4
    x y y_hat res

    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

    View Slide

  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

    1 -0.980 -0.641 0.234 0.268 1.55
    1. 残差 residuals

    View Slide

  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

    View Slide

  36. > dat_res
    # A tibble: 10 × 4
    x y y_hat res

    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

    View Slide

  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

    View Slide

  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

    View Slide

  39. 標準誤差 standard error
    𝑆𝐸 =
    𝑆𝐷
    𝑛
    𝑆𝐷 =
    1
    𝑛 − 1
    (
    !"#
    $
    (𝑦!
    − +
    𝑦)%
    標本平均の標準偏差
    (の推定値)

    View Slide

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

    1 0.241

    View Slide

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

    1 0.241 !!?

    View Slide

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

    1 0.809

    View Slide

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

    1 0.809
    deviation

    View Slide

  44. 𝑅𝑆𝐷 =
    1
    𝑛 − 𝑘 − 1
    )
    !"#
    $
    (𝑦!
    − ,
    𝑦)%
    残差標準偏差
    Residual standard deviation
    説明変数の次元(今は𝑘 = 1)
    誤差𝑢の標準偏差の期待値

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

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

    View Slide

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

    View Slide

  50. 𝑅2 =
    𝐸𝑆𝑆
    𝑇𝑆𝑆
    = 1 −
    𝑅𝑆𝑆
    𝑇𝑆𝑆
    全変動
    回帰変動 残差変動
    3. 決定係数 R-squared
    決定係数

    View Slide

  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つの変動

    View Slide

  52. > dat_dev
    # A tibble: 1 × 3
    RSS ESS TSS

    1 5.24 1.98 7.21
    > dat_dev %$% {ESS / TSS}
    [1] 0.2738825
    3つの変動

    View Slide

  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

    View Slide

  54. ぴーち?
    ああ、あれは美味しいよね。
    夏のデザートにぴったり。
    ん?えふち?
    知ってる知ってる。
    あれだよね、あれ、はは

    View Slide

  55. F分布
    𝑌 =
    𝑋!
    /𝑚
    𝑋"
    /𝑛
    互いに独⽴な𝑋!
    ~𝜒"(𝑚)、𝑋"
    ~𝜒"(𝑛)について、
    なる𝑌の確率分布をF分布と呼び、
    𝑌~𝐹(𝑚, 𝑛)
    と表す。
    (𝑚を第1⾃由度,𝑛を第2⾃由度と呼ぶ)
    カイ2乗分布

    View Slide

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

    View Slide

  57. 𝜒!分布
    正規分布𝒩(𝑢, σ")に従う独⽴な𝑘個の
    確率変数𝑋!
    , 𝑋"
    , … , 𝑋#
    について、
    𝑍 = 1
    12!
    3
    (𝑋1
    −𝑢)"
    σ"
    なる𝑍は⾃由度𝑘 − 1の𝜒"分布に従う
    𝑍~𝜒"(𝑘 − 1)

    View Slide

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

    View Slide

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

    View Slide

  60. 正規分布に従う独⽴な𝑘個の確率変数
    𝑋!
    , 𝑋"
    , … , 𝑋#
    の2乗和、すなわちコレ!?
    コレは⾃由度𝑘 − 1のχ2分布に従う。

    View Slide

  61. ⾃由度 degrees of freedom
    変数のうち⾃由な値をとりうるものの数

    View Slide

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

    View Slide

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

    View Slide

  64. ⾃由度 degrees of freedom
    変数のうち⾃由な値をとりうるものの数
    標本平均が既知の値 ̅
    𝑥になるように
    選ばれた𝑛個の標本の⾃由度は 𝒏 − 𝟏
    ̅
    𝑥 =
    1
    𝑛
    (𝑥!
    + 𝑥"
    + ⋯ + 𝑥#$!
    + 𝑥#
    )
    𝑥#
    = 𝑛 ̅
    𝑥 − (𝑥!
    + 𝑥"
    + ⋯ + 𝑥#$!
    )
    ⾃由に値をとる𝒏 − 𝟏項
    既知の
    標本平均

    View Slide

  65. ここまでおさえておいて、

    View Slide

  66. 回帰モデルの考え⽅を再確認。

    View Slide

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

    View Slide

  68. データ
    ⺟集団の(真の)
    回帰直線
    誤差 𝑢
    𝑥
    𝑦
    𝑌

    View Slide

  69. データ
    データから
    推定された
    回帰直線
    ⺟集団の(真の)
    回帰直線
    誤差
    残差
    &
    𝑦
    ε
    𝑢
    𝑦
    𝑥
    𝑌

    View Slide

  70. 𝑦 = 𝑎!
    + 𝑎"
    𝑥 + 𝑢
    𝑢 ~ !.!.#.
    𝒩(0, σ$)
    誤差
    𝑦 = )
    𝑎!
    + )
    𝑎"
    𝑥 + ε
    ε = 𝑦 − )
    𝑦
    真の回帰モデル
    残差
    最⼩⼆乗法で推定された回帰式

    View Slide

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

    View Slide

  72. (正規)回帰モデル5つの仮定
    1. 誤差の期待値はゼロ
    2. 誤差の分散は定数
    3. 誤差は互いに無相関
    4. 誤差と説明変数は無相関
    ガウス=マルコフの定理
    仮定1-4が満たされる場合、最⼩⼆乗推定量は
    線型不偏推定量の中で最⼩分散を持つ。
    5. 誤差は正規分布に従う
    + 仮定5
    最⼩⼆乗推定量は最⼩分散を持つ。

    View Slide

  73. 𝑦 = 𝑎!
    + 𝑎"
    𝑥 + 𝑢
    𝑢 ~ !.!.#.
    𝒩(0, σ$)
    誤差
    真の回帰モデル
    最⼩⼆乗法OLSでモデルの
    パラメータを推定する
    (OLS, Ordinary Last Squares)
    → 残差⼆乗和RSSの最⼩化
    (RSS, Residual Sum of Squares)

    View Slide

  74. OLSが要請する残差の性質
    1. 残差の合計値はゼロ
    2. 残差と説明変数は無相関
    ;
    %)!
    #
    ε%
    = 0
    ;
    %)!
    #
    𝑥%
    ε%
    = 0

    View Slide

  75. 1. 残差の合計値はゼロ
    2. 残差と説明変数は無相関
    ;
    %)!
    #
    ε%
    = 0
    ;
    %)!
    #
    𝑥*%
    ε%
    = 0 (𝑘 = 1, 2, … )
    説明変数の次元数
    OLSが要請する残差の性質

    View Slide

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

    View Slide

  77. 残差εの⾃由度はn − 𝑘 − 1となる。

    View Slide

  78. 𝑅𝑆𝐷 =
    1
    𝑛 − 𝑘 − 1
    )
    !"#
    $
    (𝑦!
    − ,
    𝑦)%
    残差標準偏差
    Residual standard deviation
    説明変数の次元(今は𝑘 = 1)
    𝑦 = 𝑎+
    + 𝑎!
    𝑥 + 𝑢

    View Slide

  79. 1. 残差の合計値はゼロ
    2. 残差と説明変数は無相関
    ;
    %)!
    #
    ε%
    = 0
    ;
    %)!
    #
    𝑥*%
    ε%
    = 0 (𝑘 = 1, 2, … )
    説明変数の次元数
    OLSが要請する残差の性質
    ※ 定数項/
    𝒂𝟎
    が無い回帰式では
    この性質は満たされない。
    → 残差𝜺の⾃由度は𝐧 − 𝒌。

    View Slide

  80. 𝜒!分布
    正規分布に従う独⽴な𝑘個の確率変数
    𝑋!
    , 𝑋"
    , … , 𝑋#
    について、
    𝑍 = 1
    12!
    3
    (𝑋1
    − 8
    𝑋)"
    σ"
    なる𝑍は⾃由度𝑘 − 1の𝜒"分布に従う
    𝑍~𝜒"(𝑘 − 1)

    View Slide

  81. データ数𝑛の時、残差平⽅和は(理想的には)
    ⾃由度𝑛 − 𝑘 − 1のχ2分布に従う。
    𝑅𝑆𝑆~𝜒"(𝑛 − 𝑘 − 1)

    View Slide

  82. データ数𝑛の時、全変動は(理想的には)
    ⾃由度𝑛 − 1のχ2分布に従う。
    𝑇𝑆𝑆~𝜒"(𝑛 − 1)
    &
    !"#
    $
    𝑦!
    − -
    𝑦 = 0

    View Slide

  83. 回帰変動は(理想的には)
    ⾃由度𝑘のχ2分布に従う。
    𝐸𝑆𝑆 = 𝑇𝑆𝑆 − 𝑅𝑆𝑆
    𝜒!(𝑛 − 1) 𝜒!(𝑛 − 𝑘 − 1)
    𝜒!(𝑘)
    これらが互いに独⽴なので、
    &
    !"#
    $
    *
    𝑦! − -
    𝑦 = 0

    View Slide

  84. 𝑅𝑆𝑆~𝜒((𝑛 − 𝑘 − 1) 𝑇𝑆𝑆~𝜒((𝑛 − 1) 𝐸𝑆𝑆~𝜒((𝑘)
    OLSの変動の確率分布
    残差変動 全変動 回帰変動

    View Slide

  85. 回帰モデルの適合性
    𝑅2 =
    𝐸𝑆𝑆
    𝑇𝑆𝑆
    = 1 −
    𝑅𝑆𝑆
    𝑇𝑆𝑆
    𝑅2 = 1 −
    𝑅𝑆𝑆/(𝑛 − 𝑘 − 1)
    𝑇𝑆𝑆/(𝑛 − 1)
    決定係数
    ⾃由度調整済み決定係数

    View Slide

  86. > dat_dev
    # A tibble: 1 × 3
    RSS ESS TSS

    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

    View Slide

  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

    View Slide

  88. F分布
    𝑌 =
    𝑋!
    /𝑚
    𝑋"
    /𝑛
    互いに独⽴な𝑋!
    ~𝜒"(𝑚)、𝑋"
    ~𝜒"(𝑛)について、
    なる𝑌の確率分布をF分布と呼び、
    𝑌~𝐹(𝑚, 𝑛)
    と表す。
    (𝑚を第1⾃由度,𝑛を第2⾃由度と呼ぶ)
    カイ2乗分布

    View Slide

  89. 𝑓 =
    𝐸𝑆𝑆/𝑘
    𝑅𝑆𝑆/(𝑛 − 𝑘 − 1)
    F値
    無相関(傾き=0)なら0
    分散0なら0
    𝑓~𝐹(𝑘, 𝑛 − 𝑘 − 1)
    4. ⾃由度とf検定

    View Slide

  90. > dat_dev
    # A tibble: 1 × 3
    RSS ESS TSS

    1 10.1 1.82 11.9
    3つの変動
    > df1 %
    select(contains("x")) %>% ncol()
    > df2 F値
    > dat_dev %$% {ESS * df2 / RSS / df1}
    [1] 3.017501

    View Slide

  91. > df1 %
    select(contains("x")) %>% ncol()
    > df2 > f [1] 3.017501
    F値
    > pf(q = f,
    df1 = df1,
    df2 = df2,
    lower.tail = FALSE)
    [1] 0.264185
    P値

    View Slide

  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

    View Slide

  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"

    View Slide

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

    View Slide

  95. 回帰
    dat %>%
    lm(y ~ x, data = .) %>%
    summary() %$%
    fstatistic %>% {
    pf(.[1], .[2], .[3],
    lower.tail = FALSE)
    }
    value
    0.1205754

    View Slide

  96. 回帰
    lm_pval 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

    1 0.121

    View Slide

  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検定

    View Slide

  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. 回帰係数

    View Slide

  99. 最⼩⼆乗法(線形回帰)
    argmin(-!,-")
    ;
    %)!
    #
    ε"
    𝑦 = /
    𝑎&
    + /
    𝑎#
    𝑥 + ε 残差
    &
    𝑎"
    =
    ∑#$"
    % (𝑥#
    − ̅
    𝑥)(𝑦#
    − /
    𝑦)

    #$"
    % (𝑥#
    − ̅
    𝑥)!
    &
    𝑎&
    = /
    𝑦 − &
    𝑎"
    ̅
    𝑥

    View Slide

  100. 𝑉𝑎𝑟[𝑥] =
    1
    𝑛 − 1
    6
    )*+
    ,
    (𝑥)
    − ̅
    𝑥)(
    𝐶𝑜𝑣[𝑥, 𝑦] =
    1
    𝑛 − 1
    6
    )*+
    ,
    (𝑥)
    − ̅
    𝑥)(𝑦)
    − =
    𝑦)
    最⼩⼆乗法(線形回帰)
    <
    𝑎!
    =
    ∑$%!
    & (𝑥$
    − ̅
    𝑥)(𝑦$
    − A
    𝑦)

    $%!
    & (𝑥$
    − ̅
    𝑥)"
    =
    𝐶𝑜𝑣[𝑥, 𝑦]
    𝑉𝑎𝑟[𝑥]
    <
    𝑎'
    = A
    𝑦 − <
    𝑎!
    ̅
    𝑥 = A
    𝑦 −
    𝐶𝑜𝑣[𝑥, 𝑦]
    𝑉𝑎𝑟[𝑥]
    ̅
    𝑥
    ただし、
    分散
    共分散

    View Slide

  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
    )*+
    ,
    (𝑥)
    − ̅
    𝑥)(𝑦)
    − =
    𝑦)
    共分散

    View Slide

  102. 最⼩⼆乗法(線形回帰)
    𝐶𝑜𝑣[𝑥, 𝑦] =
    1
    𝑛 − 1
    6
    )*+
    ,
    (𝑥)
    − ̅
    𝑥)(𝑦)
    − =
    𝑦)
    共分散

    View Slide

  103. 𝑉𝑎𝑟[𝑥] =
    1
    𝑛 − 1
    &
    !"#
    $
    (𝑥!
    − ̅
    𝑥)% 𝐶𝑜𝑣[𝑥, 𝑦] =
    1
    𝑛 − 1
    &
    !"#
    $
    (𝑥!
    − ̅
    𝑥)(𝑦!
    − -
    𝑦)
    # A tibble: 1 × 2
    var cov

    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))
    最⼩⼆乗法(線形回帰)
    分散 共分散

    View Slide

  104. # A tibble: 1 × 2
    var cov

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

    View Slide

  105. 最⼩⼆乗法(線形回帰)
    >
    𝑎+
    =
    𝐶𝑜𝑣[𝑥, 𝑦]
    𝑉𝑎𝑟[𝑥]
    , >
    𝑎-
    = =
    𝑦 −
    𝐶𝑜𝑣[𝑥, 𝑦]
    𝑉𝑎𝑟 𝑥
    ̅
    𝑥
    a1 a0 dat %>%
    summarise(x_bar = mean(x),
    y_bar = mean(y)) %$%
    {y_bar - a1 * x_bar}
    > a1
    [1] 0.1547321
    > a0
    [1] -0.1688236

    View Slide

  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

    View Slide

  107. ぴーち?
    ああ、あれは...

    View Slide

  108. 最⼩⼆乗法(線形回帰)
    𝑦 = 2
    𝑎?
    + 2
    𝑎@
    𝑥 + ε
    回帰係数<
    𝑎'
    と<
    𝑎!
    の仮説検定を⾏う。
    帰無仮説:<
    𝑎!
    = 0
    対⽴仮説:<
    𝑎!
    ≠ 0
    すなわち例えば<
    𝑎!
    について、
    → <
    𝑎'
    と<
    𝑎!
    の確率分布を知りたい。
    → <
    𝑎'
    と<
    𝑎!
    の期待値と分散を調べる。

    View Slide

  109. 𝑦 = 𝑎!
    + 𝑎"
    𝑥 + 𝑢
    𝑢 ~ !.!.#.
    𝒩(0, σ$)
    誤差
    𝑦 = )
    𝑎!
    + )
    𝑎"
    𝑥 + ε
    ε = 𝑦 − )
    𝑦
    真の回帰モデル
    残差
    最⼩⼆乗法で推定された回帰式

    View Slide

  110. 期待と分散の基本的性質
    𝐸 𝑋 + 𝑌 = 𝐸 𝑋 + 𝐸[𝑌]
    𝐸 𝑐𝑋 = 𝑐𝐸 𝑋
    𝑋と𝑌が独⽴の時、𝐸 𝑋𝑌 = 𝐸 𝑋 𝐸 𝑌
    𝑉𝑎𝑟 𝑋 = 𝐸 (𝑋 − 𝐸[𝑋])" = 𝐸 𝑋" − (𝐸[𝑋])"
    𝑉𝑎𝑟 𝑋 + 𝑌 = 𝑉𝑎𝑟 𝑋 + 𝑉𝑎𝑟 𝑌 + 2𝐶𝑜𝑣[𝑋, 𝑌]
    𝑉𝑎𝑟 𝑐𝑋 = 𝑐"𝑉𝑎𝑟 𝑋
    𝑉𝑎𝑟 𝑋 + 𝑐 = 𝑉𝑎𝑟 𝑋

    View Slide

  111. >
    𝑎+
    =
    ∑)*+
    , (𝑥)
    − ̅
    𝑥)(𝑦)
    − =
    𝑦)

    )*+
    , (𝑥)
    − ̅
    𝑥)(
    =
    ∑)*+
    , (𝑥)
    − ̅
    𝑥)( 𝑎-
    + 𝑎+
    𝑥)
    + 𝑢)
    − 𝑎-
    + 𝑎+
    ̅
    𝑥 + =
    𝑢 )

    )*+
    , (𝑥)
    − ̅
    𝑥)(
    =
    ∑)*+
    , 𝑥)
    − ̅
    𝑥 (𝑎+
    𝑥)
    − ̅
    𝑥 + 𝑢)
    − =
    𝑢)

    )*+
    , (𝑥)
    − ̅
    𝑥)(
    = 𝑎+
    +
    ∑)*+
    , 𝑥)
    − ̅
    𝑥 𝑢)

    )*+
    , (𝑥)
    − ̅
    𝑥)(
    +
    ∑)*+
    , 𝑥)
    − ̅
    𝑥

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

    View Slide

  112. 𝑉𝑎𝑟 <
    𝑎!
    = 𝐸 (<
    𝑎!
    − 𝐸[<
    𝑎!
    ])" = 𝐸 (<
    𝑎!
    − 𝑎!
    )"
    = 𝐸
    ∑!"#
    $ 𝑥!
    − ̅
    𝑥 𝑢!

    !"#
    $ (𝑥!
    − ̅
    𝑥)%
    %
    = 𝐸
    ∑!"#
    $ 𝑥!
    − ̅
    𝑥 %𝑢!
    %

    !"#
    $ (𝑥!
    − ̅
    𝑥)% %
    + 𝐸
    ∑!"#
    $ ∑&"#
    $ 𝑥!
    − ̅
    𝑥 𝑥&
    − ̅
    𝑥 𝑢!
    𝑢&

    !"#
    $ (𝑥!
    − ̅
    𝑥)% %
    𝑖 ≠ 𝑗
    =
    ∑!"#
    $ 𝑥!
    − ̅
    𝑥 %𝐸 𝑢!
    %

    !"#
    $ (𝑥!
    − ̅
    𝑥)% %
    +
    ∑!"#
    $ ∑&"#
    $ 𝑥!
    − ̅
    𝑥 𝑥&
    − ̅
    𝑥 𝐸[𝑢!
    𝑢&
    ]

    !"#
    $ (𝑥!
    − ̅
    𝑥)% %
    =
    𝐸 𝑢% − (𝐸[𝑢])%

    !"#
    $ (𝑥!
    − ̅
    𝑥)%
    =
    𝑉𝑎𝑟[𝑢%]

    !"#
    $ (𝑥!
    − ̅
    𝑥)%
    =
    σ%

    !"#
    $ (𝑥!
    − ̅
    𝑥)%
    𝐸 𝑢 = 0 𝑉𝑎𝑟 𝑎! = 𝐸 𝑎! − (𝐸 𝑎 )!
    5. 回帰係数のt検定

    View Slide

  113. 𝐸 <
    𝑎!
    = 𝑎!
    𝑉𝑎𝑟[<
    𝑎!
    ] =
    σ"

    $%!
    & (𝑥$
    − ̅
    𝑥)"
    頑張って計算した結果、
    >
    𝑎+
    = 𝑎+
    +
    ∑)*+
    , 𝑥)
    − ̅
    𝑥

    )*+
    , (𝑥)
    − ̅
    𝑥)(
    =
    𝑢 +
    ∑)*+
    , 𝑥)
    − ̅
    𝑥 𝑢)

    )*+
    , (𝑥)
    − ̅
    𝑥)(
    また式(1)
    𝑢~ 𝒩(0, σ()
    より *
    𝑎#
    は正規分布に従い、
    >
    𝑎+
    ~ 𝒩 𝑎+
    ,
    σ(

    )*+
    , (𝑥)
    − ̅
    𝑥)(

    >
    𝑎+
    − 𝑎+
    σ(/ ∑
    )*+
    , (𝑥)
    − ̅
    𝑥)(
    ~𝒩 0,1
    5. 回帰係数のt検定

    View Slide

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

    View Slide

  115. 𝑅𝑆𝐷 =
    1
    𝑛 − 𝑘 − 1
    )
    !"#
    $
    (𝑦!
    − ,
    𝑦)%
    残差標準偏差
    Residual standard deviation
    誤差𝑢の標準偏差の期待値

    View Slide

  116. F
    𝑎!
    − 𝑎!
    𝑅𝑆𝐷"/ ∑
    %)!
    # (𝑥%
    − ̅
    𝑥)"
    =
    F
    𝑎!
    − 𝑎!
    𝑠𝑒(F
    𝑎!
    )
    ~ ?
    F
    𝑎!
    − 𝑎!
    σ"/ ∑
    %)!
    # (𝑥%
    − ̅
    𝑥)"
    ~𝒩 0,1
    置き換える
    確率分布は?
    𝑠𝑒 >
    𝑎+
    ∶=
    𝑅𝑆𝐷(

    )*+
    , (𝑥)
    − ̅
    𝑥)(
    5. 回帰係数のt検定

    View Slide

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

    View Slide

  118. t分布
    互いに独⽴な、
    標準正規分布𝒩(0,1)に従う確率変数𝑋と、
    カイ⼆乗分布χ"(𝑛)に従う確率変数𝑌について
    𝑍 =
    𝑋
    𝑌/𝑛
    なる𝑍の確率分布を𝑡分布と呼び、
    𝑍~𝑡(𝑛)
    と表す。(𝑛を⾃由度と呼ぶ)

    View Slide

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

    View Slide

  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 5. 回帰係数のt検定

    View Slide

  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検定

    View Slide

  122. t分布

    View Slide

  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")

    View Slide

  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検定

    View Slide

  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検定

    View Slide

  126. F
    𝑎+
    − 𝑎+
    𝑠𝑒(F
    𝑎+
    )
    ~𝑡(𝑛 − 𝑘 − 1)
    F
    𝑎+
    − 𝑎+
    σ"の関数!?
    ~𝒩 0,1
    σを𝑅𝑆𝐷に置き換える
    5. 回帰係数のt検定

    View Slide

  127. >
    𝑎-
    = =
    𝑦 − >
    𝑎+
    ̅
    𝑥
    = 𝑎-
    + 𝑎+
    ̅
    𝑥 − =
    𝑢 − >
    𝑎+
    ̅
    𝑥
    𝐸[>
    𝑎-
    ] = 𝐸[𝑎-
    + 𝑎+
    ̅
    𝑥 − =
    𝑢 − >
    𝑎+
    ̅
    𝑥]
    = 𝑎-
    + 𝑎+
    ̅
    𝑥 − 0 − 𝑎+
    ̅
    𝑥
    = 𝑎-
    𝑉𝑎𝑟 >
    𝑎-
    = 𝐸 (>
    𝑎-
    −𝐸[>
    𝑎-
    ])( = 𝐸 (>
    𝑎-
    −𝑎-
    )(
    = σ(
    1
    𝑛
    +
    ̅
    𝑥(

    )*+
    , (𝑥)
    − ̅
    𝑥)(
    = …
    5. 回帰係数のt検定

    View Slide

  128. F
    𝑎+
    − 𝑎+
    𝑠𝑒(F
    𝑎+
    )
    ~𝑡(𝑛 − 𝑘 − 1)
    F
    𝑎+
    − 𝑎+
    σ"
    1
    𝑛
    +
    ̅
    𝑥"

    %)!
    # (𝑥%
    − ̅
    𝑥)"
    ~𝒩 0,1
    σを𝑅𝑆𝐷に置き換える
    𝑠𝑒 *
    𝑎4 ≔ 𝑅𝑆𝐷
    1
    𝑛
    +
    ̅
    𝑥%
    ∑!"#
    $ (𝑥! − ̅
    𝑥)%
    5. 回帰係数のt検定

    View Slide

  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 5. 回帰係数のt検定

    View Slide

  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検定

    View Slide

  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検定

    View Slide

  132. Pipe演算⼦: %>%
    Dollar演算⼦: %$%
    線形回帰: lm()
    names属性: names()
    要約統計量の算出:summarize()
    要約統計量の算出:summarize_at()
    確率分布関連:pt(), dt(), pf(), df()
    技術的復習

    View Slide

  133. Enjoy!!

    View Slide