TokyoR#101_RegressionAnalysis

Slide 1

Slide 1 text

#101 @kilometer00 2022.09.17 Regression Analysis

Slide 2

Slide 2 text

Who！？誰だ？

Slide 3

Slide 3 text

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

Slide 4

Slide 4 text

宣伝!! （書籍の翻訳に参加しました。）第五刷！！

Slide 5

Slide 5 text

No content

Slide 6

Slide 6 text

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

Slide 7

Slide 7 text

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

Slide 8

Slide 8 text

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

Slide 9

Slide 9 text

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

Slide 10

Slide 10 text

回帰していますか？みなさん、

Slide 11

Slide 11 text

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

Slide 12

Slide 12 text

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

Slide 13

Slide 13 text

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

Slide 14

Slide 14 text

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

Slide 15

Slide 15 text

線形モデル 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

Slide 16

Slide 16 text

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

Slide 17

Slide 17 text

> 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

Slide 18

Slide 18 text

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

Slide 19

Slide 19 text

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

Slide 20

Slide 20 text

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

Slide 21

Slide 21 text

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

Slide 22

Slide 22 text

今⽇のミッション > 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 この中に書かれた数字を全て計算する。

Slide 23

Slide 23 text

Slide 24

Slide 24 text

Slide 25

Slide 25 text

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

Slide 26

Slide 26 text

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

Slide 27

Slide 27 text

Slide 28

Slide 28 text

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

Slide 29

Slide 29 text

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

Slide 30

Slide 30 text

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

Slide 31

Slide 31 text

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

Slide 32

Slide 32 text

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

Slide 33

Slide 33 text

> 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

Slide 34

Slide 34 text

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

Slide 35

Slide 35 text

Slide 36

Slide 36 text

> 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

Slide 37

Slide 37 text

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

Slide 38

Slide 38 text

Slide 39

Slide 39 text

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

Slide 40

Slide 40 text

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

Slide 41

Slide 41 text

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

Slide 42

Slide 42 text

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

Slide 43

Slide 43 text

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

Slide 44

Slide 44 text

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

Slide 45

Slide 45 text

> ?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

Slide 46

Slide 46 text

Slide 47

Slide 47 text

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

Slide 48

Slide 48 text

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

Slide 49

Slide 49 text

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

Slide 50

Slide 50 text

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

Slide 51

Slide 51 text

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

Slide 52

Slide 52 text

> dat_dev # A tibble: 1 × 3 RSS ESS TSS 1 5.24 1.98 7.21 > dat_dev %$% {ESS / TSS} [1] 0.2738825 3つの変動

Slide 53

Slide 53 text

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

Slide 54

Slide 54 text

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

Slide 55

Slide 55 text

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

Slide 56

Slide 56 text

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

Slide 57

Slide 57 text

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

Slide 58

Slide 58 text

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

Slide 59

Slide 59 text

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

Slide 60

Slide 60 text

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

Slide 61

Slide 61 text

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

Slide 62

Slide 62 text

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

Slide 63

Slide 63 text

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

Slide 64

Slide 64 text

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

Slide 65

Slide 65 text

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

Slide 66

Slide 66 text

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

Slide 67

Slide 67 text

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

Slide 68

Slide 68 text

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

Slide 69

Slide 69 text

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

Slide 70

Slide 70 text

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

Slide 71

Slide 71 text

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

Slide 72

Slide 72 text

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

Slide 73

Slide 73 text

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

Slide 74

Slide 74 text

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

Slide 75

Slide 75 text

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

Slide 76

Slide 76 text

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

Slide 77

Slide 77 text

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

Slide 78

Slide 78 text

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

Slide 79

Slide 79 text

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

Slide 80

Slide 80 text

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

Slide 81

Slide 81 text

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

Slide 82

Slide 82 text

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

Slide 83

Slide 83 text

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

Slide 84

Slide 84 text

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

Slide 85

Slide 85 text

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

Slide 86

Slide 86 text

> 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

Slide 87

Slide 87 text

Slide 88

Slide 88 text

Slide 89

Slide 89 text

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

Slide 90

Slide 90 text

> 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

Slide 91

Slide 91 text

> 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値

Slide 92

Slide 92 text

回帰モデルの適合性 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

Slide 93

Slide 93 text

回帰 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"

Slide 94

Slide 94 text

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

Slide 95

Slide 95 text

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

Slide 96

Slide 96 text

回帰 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

Slide 97

Slide 97 text

Slide 98

Slide 98 text

Slide 99

Slide 99 text

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

Slide 100

Slide 100 text

𝑉𝑎𝑟[𝑥] = 1 𝑛 − 1 6 )*+ , (𝑥) − ̅ 𝑥)( 𝐶𝑜𝑣[𝑥, 𝑦] = 1 𝑛 − 1 6 )*+ , (𝑥) − ̅ 𝑥)(𝑦) − = 𝑦) 最⼩⼆乗法（線形回帰） < 𝑎! = ∑$%! & (𝑥$ − ̅ 𝑥)(𝑦$ − A 𝑦) ∑ $%! & (𝑥$ − ̅ 𝑥)" = 𝐶𝑜𝑣[𝑥, 𝑦] 𝑉𝑎𝑟[𝑥] < 𝑎' = A 𝑦 − < 𝑎! ̅ 𝑥 = A 𝑦 − 𝐶𝑜𝑣[𝑥, 𝑦] 𝑉𝑎𝑟[𝑥] ̅ 𝑥 ただし、分散共分散

Slide 101

Slide 101 text

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

Slide 102

Slide 102 text

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

Slide 103

Slide 103 text

𝑉𝑎𝑟[𝑥] = 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)) 最⼩⼆乗法（線形回帰）分散共分散

Slide 104

Slide 104 text

# 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 & !"# $ (𝑥! − ̅ 𝑥)(𝑦! − - 𝑦) 分散共分散

Slide 105

Slide 105 text

最⼩⼆乗法（線形回帰） > 𝑎+ = 𝐶𝑜𝑣[𝑥, 𝑦] 𝑉𝑎𝑟[𝑥] , > 𝑎- = = 𝑦 − 𝐶𝑜𝑣[𝑥, 𝑦] 𝑉𝑎𝑟 𝑥 ̅ 𝑥 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

Slide 106

Slide 106 text

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

Slide 107

Slide 107 text

ぴーち？ああ、あれは...

Slide 108

Slide 108 text

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

Slide 109

Slide 109 text

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

Slide 110

Slide 110 text

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

Slide 111

Slide 111 text

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

Slide 112

Slide 112 text

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

Slide 113

Slide 113 text

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

Slide 114

Slide 114 text

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

Slide 115

Slide 115 text

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

Slide 116

Slide 116 text

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

Slide 117

Slide 117 text

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

Slide 118

Slide 118 text

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

Slide 119

Slide 119 text

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

Slide 120

Slide 120 text

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

Slide 121

Slide 121 text

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

Slide 122

Slide 122 text

t分布

Slide 123

Slide 123 text

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

Slide 124

Slide 124 text

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

Slide 125

Slide 125 text

Slide 126

Slide 126 text

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

Slide 127

Slide 127 text

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

Slide 128

Slide 128 text

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

Slide 129

Slide 129 text

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

Slide 130

Slide 130 text

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

Slide 131

Slide 131 text

回帰を(ほぼ)完全に理解した > 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検定

Slide 132

Slide 132 text

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

Slide 133

Slide 133 text

Enjoy!!