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 full-size slide

  2. Who!?
    誰だ?

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

  11. library(tidyverse)
    library(magrittr)
    set.seed(1)
    準備
    N <- 10 # データ数
    a <- 0.1 # 傾き
    dat <-
    tibble(x = c(1:N),
    y = a * x + rnorm(N))

    View full-size slide

  12. > 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 full-size slide

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

    View full-size slide

  14. 線形モデル linear model
    dat <-
    tibble(x = c(1:N),
    y = a * x + rnorm(N))
    lm(y ~ x, data = dat)
    Call:
    lm(formula = y ~ x, data = dat)
    Coefficients:
    (Intercept) x
    -0.1688 0.1547

    View full-size slide

  15. fit_lm <-
    dat %>%
    lm(y ~ x, data = .)
    > fit_lm %>% names()
    [1] "coefficients" "residuals"
    [3] "effects" "rank"
    [5] "fitted.values" "assign"
    [7] "qr" "df.residual"
    [9] "xlevels" "call"
    [11] "terms" "model"
    線形モデル linear model

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size 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 full-size 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
    1. 残差 residuals

    View full-size slide

  23. fit_lm <-
    dat %>%
    lm(y ~ x, data = .)
    > fit_lm %>% names()
    [1] "coefficients" "residuals"
    [3] "effects" "rank"
    [5] "fitted.values" "assign"
    [7] "qr" "df.residual"
    [9] "xlevels" "call"
    [11] "terms" "model"
    1. 残差 residuals

    View full-size slide

  24. res <- fit_lm$residuals
    > res
    1 2 3
    -0.5123623 0.2430028 -0.8310012
    4 5 6
    1.5451761 0.2246710 -0.9800372
    7 8 9
    0.2731282 0.4692918 0.2520164
    10
    -0.6838854
    1. 残差 residuals

    View full-size slide

  25. res <- fit_lm$residuals
    > res %>% median()
    [1] 0.2187553
    > res %>% min()
    [1] -2.390729
    > res %>% max()
    [1] 1.481385
    > res %>%
    + quantile(probs = c(0.25, 0.75))
    25% 75%
    -0.3329012 0.5308412
    1. 残差 residuals

    View full-size slide

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

    View full-size slide

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

    View full-size slide

  28. a0 <- fit_lm$coefficients[1]
    a1 <- fit_lm$coefficients[2]
    dat %>%
    mutate(y_hat = a0 + a1 * x,
    res = y - y_hat)
    # A tibble: 10 × 4
    x y y_hat res

    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 full-size slide

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

    View full-size slide

  30. 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 full-size slide

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

    View full-size slide

  32. > 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 full-size slide

  33. dat_res %>%
    summarize(
    across(res,
    list(Min = ~ min(.),
    Q1 = ~ quantile(., probs = c(0.25)),
    Median = ~ median(.),
    Q3 = ~ quantile(., probs = c(0.75)),
    Max = ~ max(.))))
    # A tibble: 1 × 5
    res_Min res_Q1 res_Median res_Q3 res_Max

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

    View full-size slide

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

    View full-size slide

  35. > 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 full-size slide

  36. ggplot(data = dat_res) +
    aes(x = x, y = y) +
    geom_path(aes(y = y_hat),
    color = "red",
    linetype = "dotted") +
    geom_segment(aes(xend = x,
    yend = y_hat),
    color = "red") +
    geom_point()
    1. 残差 residuals

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    1 0.241

    View full-size slide

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

    1 0.241 !!?

    View full-size slide

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

    1 0.809

    View full-size 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
    deviation

    View full-size slide

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

    View full-size slide

  44. > ?stats::sigma
    Description
    Extract the estimated standard deviation of
    the errors, the “residual standard deviation”
    (misnamed also “residual standard error”, e.g.,
    in summary.lm()'s output, from a fitted model).
    2. 残差標準偏差 RSD

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

  50. dat_dev <-
    dat_res %>%
    mutate(yhat_mean = mean(y_hat),
    RSS = y - y_hat,
    ESS = y_hat - yhat_mean,
    TSS = y - yhat_mean) %>%
    summarize_at(c("RSS", "ESS", "TSS"),
    ~ sum(. ^ 2))
    3つの変動

    View full-size slide

  51. > 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 full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

  85. > 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 full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    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

    View full-size slide

  90. > df1 <- dat %>%
    select(contains("x")) %>% ncol()
    > df2 <- N – df1 - 1
    > f <- dat_dev %$% {ESS * df2 / RSS / df1}
    [1] 3.017501
    F値
    > pf(q = f,
    df1 = df1,
    df2 = df2,
    lower.tail = FALSE)
    [1] 0.264185
    P値

    View full-size slide

  91. 回帰モデルの適合性
    F(1,8)
    f = 3.017501
    p = 0.264185
    tibble(x = seq(0, 5, by = 0.01),
    df = df(x, df1, df2))
    lower.tail=TRUE lower.tail=FALSE

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

  95. 回帰
    lm_pval <- function(fit){
    summary(fit) %$%
    fstatistic %>% {
    pf(.[1], .[2], .[3],
    lower.tail = FALSE)}
    }
    dat_fit %>%
    mutate(pval = map_dbl(fit, lm_pval))
    # A tibble: 1 × 3
    data fit pval

    1 0.121

    View full-size slide

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

    View full-size 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
    5. 回帰係数

    View full-size slide

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

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

    View full-size slide

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

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

  102. 𝑉𝑎𝑟[𝑥] =
    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 full-size slide

  103. # 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 full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

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

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

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

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

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

    View full-size slide

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

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

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

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

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

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

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

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

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

    View full-size slide

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

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

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

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

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

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

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

    View full-size slide

  120. > a1_se
    [1] 0.08907516
    > t_value
    [1] 1.737096
    Coefficients:
    Estimate Std. Error t value Pr(>|t|)
    (Intercept) -0.16882 0.55270 -0.305 0.768
    x 0.15473 0.08908 1.737 0.121
    5. 回帰係数のt検定

    View full-size slide

  121. 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 full-size slide

  122. > 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 full-size slide

  123. > 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 full-size slide

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

    View full-size slide

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

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

    View full-size slide

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

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

    View full-size slide

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

    View full-size slide

  128. > 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 full-size slide

  129. 回帰を(ほぼ)完全に理解した
    > 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 full-size slide

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

    View full-size slide