TokyoR#101_RegressionAnalysis

#101
@kilometer00
2022.09.17
Regression Analysis

View Slide

Who！？
誰だ？

View Slide

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

View Slide

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

View Slide

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

> 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

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

View Slide

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

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"

View Slide

> 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

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

View Slide

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

View Slide

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

View Slide

例えば、
Call:
Residuals:
-0.9800 -0.6410 0.2338 0.2678 1.5452
Coefficients:
(Intercept) -0.16882 0.55270 -0.305 0.768
x 0.15473 0.08908 1.737 0.121
この中に誤記があります(知ってます？)。

View Slide

今⽇のミッション
Call:
Residuals:
-0.9800 -0.6410 0.2338 0.2678 1.5452
Coefficients:
(Intercept) -0.16882 0.55270 -0.305 0.768
x 0.15473 0.08908 1.737 0.121
この中に書かれた数字を全て計算する。

View Slide

Call:
Residuals:
-0.9800 -0.6410 0.2338 0.2678 1.5452
Coefficients:
(Intercept) -0.16882 0.55270 -0.305 0.768
x 0.15473 0.08908 1.737 0.121
1. 残差 residuals

View Slide

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

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

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

Call:
Residuals:
-0.9800 -0.6410 0.2338 0.2678 1.5452
Coefficients:
(Intercept) -0.16882 0.55270 -0.305 0.768
x 0.15473 0.08908 1.737 0.121
1. 残差 residuals

View Slide

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

View Slide

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

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

View Slide

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

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

View Slide

> 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

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

Call:
Residuals:
-0.9800 -0.6410 0.2338 0.2678 1.5452
Coefficients:
(Intercept) -0.16882 0.55270 -0.305 0.768
x 0.15473 0.08908 1.737 0.121
1. 残差 residuals

View Slide

> 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

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

Call:
Residuals:
-0.9800 -0.6410 0.2338 0.2678 1.5452
Coefficients:
(Intercept) -0.16882 0.55270 -0.305 0.768
x 0.15473 0.08908 1.737 0.121

View Slide

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

View Slide

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

1 0.241

View Slide

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

1 0.241 !!?

View Slide

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

1 0.809

View Slide

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

1 0.809
deviation

View Slide

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

View Slide

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

Call:
Residuals:
-0.9800 -0.6410 0.2338 0.2678 1.5452
Coefficients:
(Intercept) -0.16882 0.55270 -0.305 0.768
x 0.15473 0.08908 1.737 0.121
deviation
2. 残差標準偏差 RSD

View Slide

3. 決定係数 R-squared
Call:
Residuals:
-0.9800 -0.6410 0.2338 0.2678 1.5452
Coefficients:
(Intercept) -0.16882 0.55270 -0.305 0.768
x 0.15473 0.08908 1.737 0.121
deviation

View Slide

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

View Slide

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

View Slide

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

View Slide

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

> 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

4. ⾃由度とf検定
Call:
Residuals:
-0.9800 -0.6410 0.2338 0.2678 1.5452
Coefficients:
(Intercept) -0.16882 0.55270 -0.305 0.768
x 0.15473 0.08908 1.737 0.121
deviation

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

𝑋!
, 𝑋"
, … , 𝑋#
データ数𝑘の時、残差⼆乗和は
⾃由度𝑘 − 1のχ2分布に従う？ → 従わない

View Slide

𝑋!
, 𝑋"
, … , 𝑋#
コレは⾃由度𝑘 − 1のχ2分布に従う。

View Slide

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

View Slide

⺟集団
標本𝑥
(𝑛個)
⾃由に選ばれた𝑛個の標本の⾃由度は𝒏で、
標本の平均値は標本平均 ̅
𝑥。
では、標本平均が既知の値 ̅
𝑥になるように
選ばれた𝑛個の標本の⾃由度は？

View Slide

⺟集団
標本𝑥
(𝑛個)
⾃由に選ばれた𝑛個の標本の⾃由度は𝒏で、
標本の平均値は標本平均 ̅
𝑥。
では、標本平均が既知の値 ̅
選ばれた𝑛個の標本の⾃由度は？ → 𝒏 − 𝟏

View Slide

標本平均が既知の値 ̅
選ばれた𝑛個の標本の⾃由度は 𝒏 − 𝟏
̅
𝑥 =
1
𝑛
(𝑥!
+ 𝑥"
+ ⋯ + 𝑥#$!
+ 𝑥#
)
𝑥#
= 𝑛 ̅
𝑥 − (𝑥!
+ 𝑥"
+ ⋯ + 𝑥#$!
)
⾃由に値をとる𝒏 − 𝟏項
既知の
標本平均

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

;
%)!
#
ε%
= 0
;
%)!
#
𝑥*%
ε%
= 0 (𝑘 = 1, 2, … )
説明変数の次元数

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

𝜒!分布
𝑋!
, 𝑋"
, … , 𝑋#
について、
𝑍 = 1
12!
3
(𝑋1
− 8
𝑋)"
σ"
なる𝑍は⾃由度𝑘 − 1の𝜒"分布に従う
𝑍~𝜒"(𝑘 − 1)

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

> 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

Call:
Residuals:
-0.9800 -0.6410 0.2338 0.2678 1.5452
Coefficients:
(Intercept) -0.16882 0.55270 -0.305 0.768
x 0.15473 0.08908 1.737 0.121
deviation

View Slide

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

View Slide

𝑓 =
𝐸𝑆𝑆/𝑘
𝑅𝑆𝑆/(𝑛 − 𝑘 − 1)
F値
無相関(傾き=0)なら0
分散0なら0
𝑓~𝐹(𝑘, 𝑛 − 𝑘 − 1)

View Slide

> 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

> 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

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

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

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

View Slide

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

View Slide

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

Call:
Residuals:
-0.9800 -0.6410 0.2338 0.2678 1.5452
Coefficients:
(Intercept) -0.16882 0.55270 -0.305 0.768
x 0.15473 0.08908 1.737 0.121
deviation

View Slide

Call:
Residuals:
-0.9800 -0.6410 0.2338 0.2678 1.5452
Coefficients:
(Intercept) -0.16882 0.55270 -0.305 0.768
x 0.15473 0.08908 1.737 0.121
deviation
5. 回帰係数

View Slide

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

View Slide

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

View Slide

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

1
𝑛 − 1
6
)*+
,
(𝑥)
− ̅
𝑥)(𝑦)
− =
𝑦)
共分散

View Slide

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

# 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

>
𝑎+
=
𝑉𝑎𝑟[𝑥]
, >
𝑎-
= =
𝑦 −
𝑉𝑎𝑟 𝑥
̅
𝑥
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

5. 回帰係数のt検定
Call:
Residuals:
-0.9800 -0.6410 0.2338 0.2678 1.5452
Coefficients:
(Intercept) -0.16882 0.55270 -0.305 0.768
x 0.15473 0.08908 1.737 0.121
deviation

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

>
𝑎+
=
∑)*+
, (𝑥)
− ̅
𝑥)(𝑦)
− =
𝑦)
∑
)*+
, (𝑥)
− ̅
𝑥)(
=
∑)*+
, (𝑥)
− ̅
𝑥)( 𝑎-
+ 𝑎+
𝑥)
+ 𝑢)
− 𝑎-
+ 𝑎+
̅
𝑥 + =
𝑢 )
∑
)*+
, (𝑥)
− ̅
𝑥)(
=
∑)*+
, 𝑥)
− ̅
𝑥 (𝑎+
𝑥)
− ̅
𝑥 + 𝑢)
− =
𝑢)
∑
)*+
, (𝑥)
− ̅
𝑥)(
= 𝑎+
+
∑)*+
, 𝑥)
− ̅
𝑥 𝑢)
∑
)*+
, (𝑥)
− ̅
𝑥)(
+
∑)*+
, 𝑥)
− ̅
𝑥
∑
)*+
, (𝑥)
− ̅
𝑥)(
=
𝑢
𝐸[F
𝑎!
] = 𝑎!
従って、
…式(1)

View Slide

𝑉𝑎𝑟 <
𝑎!
= 𝐸 (<
𝑎!
− 𝐸[<
𝑎!
])" = 𝐸 (<
𝑎!
− 𝑎!
)"
= 𝐸
∑!"#
$ 𝑥!
− ̅
𝑥 𝑢!
∑
!"#
$ (𝑥!
− ̅
𝑥)%
%
= 𝐸
∑!"#
$ 𝑥!
− ̅
𝑥 %𝑢!
%
∑
!"#
$ (𝑥!
− ̅
𝑥)% %
+ 𝐸
∑!"#
$ ∑&"#
$ 𝑥!
− ̅
𝑥 𝑥&
− ̅
𝑥 𝑢!
𝑢&
∑
!"#
$ (𝑥!
− ̅
𝑥)% %
𝑖 ≠ 𝑗
=
∑!"#
$ 𝑥!
− ̅
𝑥 %𝐸 𝑢!
%
∑
!"#
$ (𝑥!
− ̅
𝑥)% %
+
∑!"#
$ ∑&"#
$ 𝑥!
− ̅
𝑥 𝑥&
− ̅
𝑥 𝐸[𝑢!
𝑢&
]
∑
!"#
$ (𝑥!
− ̅
𝑥)% %
=
𝐸 𝑢% − (𝐸[𝑢])%
∑
!"#
$ (𝑥!
− ̅
𝑥)%
=
𝑉𝑎𝑟[𝑢%]
∑
!"#
$ (𝑥!
− ̅
𝑥)%
=
σ%
∑
!"#
$ (𝑥!
− ̅
𝑥)%
𝐸 𝑢 = 0 𝑉𝑎𝑟 𝑎! = 𝐸 𝑎! − (𝐸 𝑎 )!

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

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

View Slide

> 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

> a1_se
[1] 0.08907516
> t_value
[1] 1.737096
Coefficients:
(Intercept) -0.16882 0.55270 -0.305 0.768
x 0.15473 0.08908 1.737 0.121

View Slide

t分布

View Slide

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

> pt(t_value, 8, lower.tail = F) * 2
[1] 0.1205754
Coefficients:
(Intercept) -0.16882 0.55270 -0.305 0.768
x 0.15473 0.08908 1.737 0.121

View Slide

Call:
Residuals:
-0.9800 -0.6410 0.2338 0.2678 1.5452
Coefficients:
(Intercept) -0.16882 0.55270 -0.305 0.768
x 0.15473 0.08908 1.737 0.121
deviation

View Slide

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

View Slide

>
𝑎-
= =
𝑦 − >
𝑎+
̅
𝑥
= 𝑎-
+ 𝑎+
̅
𝑥 − =
𝑢 − >
𝑎+
̅
𝑥
𝐸[>
𝑎-
] = 𝐸[𝑎-
+ 𝑎+
̅
𝑥 − =
𝑢 − >
𝑎+
̅
𝑥]
= 𝑎-
+ 𝑎+
̅
𝑥 − 0 − 𝑎+
̅
𝑥
= 𝑎-
𝑉𝑎𝑟 >
𝑎-
= 𝐸 (>
𝑎-
−𝐸[>
𝑎-
])( = 𝐸 (>
𝑎-
−𝑎-
)(
= σ(
1
𝑛
+
̅
𝑥(
∑
)*+
, (𝑥)
− ̅
𝑥)(
= …

View Slide

F
𝑎+
− 𝑎+
𝑠𝑒(F
𝑎+
)
~𝑡(𝑛 − 𝑘 − 1)
F
𝑎+
− 𝑎+
σ"
1
𝑛
+
̅
𝑥"
∑
%)!
# (𝑥%
− ̅
𝑥)"
~𝒩 0,1
σを𝑅𝑆𝐷に置き換える
𝑠𝑒 *
𝑎4 ≔ 𝑅𝑆𝐷
1
𝑛
+
̅
𝑥%
∑!"#
$ (𝑥! − ̅
𝑥)%

View Slide

> 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

> pt(t_a0, 8, lower.tail = T) * 2
1] 0.767817
Coefficients:
(Intercept) -0.16882 0.55270 -0.305 0.768
x 0.15473 0.08908 1.737 0.121

View Slide

回帰を(ほぼ)完全に理解した
Call:
Residuals:
-0.9800 -0.6410 0.2338 0.2678 1.5452
Coefficients:
(Intercept) -0.16882 0.55270 -0.305 0.768
x 0.15473 0.08908 1.737 0.121
deviation
1
2
4 5
1. 残差
2. 残差標準誤差
3. 決定係数 4. f検定
5. 回帰係数 6. t検定

View Slide

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

View Slide

Enjoy!!

View Slide

TokyoR#101_RegressionAnalysis

TokyoR#101_RegressionAnalysis

kilometer

More Decks by kilometer

Other Decks in Technology

Featured

Transcript