応用確率統計学(第3回) - STReP in IRIDeS and CNEAS

第5回 回帰モデル
•回帰モデル
–単回帰モデル、重回帰モデル
–非線形式の回帰分析
•最小二乗法(記述統計)
–回帰式の適合度のチェック(分散分析)
•最小二乗推定量の性質(推測統計)
–パラメータ推定量の不偏性
–パラメータ推定量の確率分布とt検定
1
回帰モデル
• いくつかの変数間に相関関係が存在
• ある変数の値を、別の変数を用いて説明
従属変数、目的変数 説明式を作成
独立変数、説明変数
x
被説明変数
変数Y,実現値yi
推計値yi=f(xi)
変数X,実現値xi
i
Y
Y
yi
yˆi
Yˆi  f ( X i )
Yˆi  f ( X i , Zi )
Z
xi
X
X
2
記述統計学と推測統計学
記述統計学
多数データの
数学的要約
・記述
母集団の
データ
推測統計学
(仮想的)
母集団
「標本統計量」
無作為
抽出
標本集団
のデータ
少数データの
数学的要約
・記述
確率的推測・記述
「母集団の確率分布・母数」
「標本統計量の確率分布」
3
目的変数を[記述]する
従属変数、目的変数 説明式を作成
独立変数、説明変数
被説明変数
変数Y,実現値yi
推計値yi=f(xi)
変数X,実現値xi
Y
yi
yˆi

xi
記述統計学の立場
関数形f(X)をうまく決めて、
Yˆi  f ( X i ) Yをうまく記述したい
(説明できなかった部分の)
残差・誤差
残差を小さくしたい
ˆ
i  yi  yi (2乗和で評価するので)
残差の二乗和を小さくしたい
X
最小二乗法
min  i2   ( yi  yˆi )2   ( yi  f ( X i )) 2
i
i
i
4
最小二乗法(単回帰モデル)
• 線形モデル yi  a  bxi
n
n
2
2
S



(
y

a

bx
)
• 残差平方和 E  i  i
i
i 1
– 未知数で微分すると
S E
 2 ( yi  a  bxi )  0
a
i 1
S E
 2 ( yi  a  bxi ) xi  0
b
• 正規方程式 上式を整理して
y
i
 a n  b xi  0
 xi yi  a xi  b xi  0
2
– それらを解いて
aˆ 
 xi
 y x x y
n x  ( x )
2
i
2
i
i
i
2
i
i
n  xi yi   xi  yi
ˆ
b
2
n xi  ( xi ) 2
5
最小二乗法(重回帰モデル)
yi  a  bwi  czi
• 線形モデル
• 残差平方和
–
n
n
S E    i   ( yi  a  bwi  czi ) 2
2
i 1
未知数で微分すると
i 1
SE a  2( yi  a  bwi  czi )  0
SE b  2( yi  a  bwi  czi )wi  0
SE c  2( yi  a  bwi  czi ) zi  0
• 正規方程式 上式を整理して
 y  a n  b w  c  z  0
 w y  a  w  b w  c  w z  0
 z y  a  z  b z w  c  z  0
i
i
i
2
i
i
i
i
i
i
i
i
i
i i
2
i
– 単回帰モデルと同様の連立一次方程式が出来る
6
最小二乗法(重回帰モデル)行列表示
yˆi  a  bwi  czi
• 線形モデル yˆ  X
  y  X
• 残差
n
• 残差平方和 S E    i 2   t  ( y  X )t ( y X )
i 1
– 未知数で微分すると
 2 X ( y  Xˆ )  0
• 正規方程式
t
t
t
1
t
ˆ
ˆ
X y  X X より、   ( X X ) ( X y)
t
7
行列を用いた重回帰式の計算例
lm(formula = 強さ ~ 薬剤A + 薬剤B)
Coefficients:
(Intercept)
薬剤A
薬剤B
39.728
6.433
1.750
8
Rによる計算例
Call:
#データフレームを作成する
lm(formula = 強さ ~ 薬剤A)
concrete <- data.frame(Call:
+ 強度 =c(40.0,45.8,53.0,42.3,46.7,54.1,44.2,47.1,58.0),
lm(formula
Residuals: = 強さ ~ 薬剤A + 薬剤B)
+ 混和A =c(0,1,2,0,1,2,0,1,2),
Min
1Q Median
3Q
Max
Residuals:
-2.1111 -1.3444 -0.8111 0.8222 3.6556
+ 混和B =c(0,0,0,1,1,1,2,2,2)
Min
1Q Median
3Q
Max
)
-2.5611
-0.3611 0.2722 0.8222 1.9056
Coefficients:
#データフレームから変数に代入(付置)する
Estimate Std. Error t value Pr(>|t|)
強さ <-concrete$強度 Coefficients:
(Intercept)
41.478
1.128 36.761 2.86e-09 ***
Estimate
Std.
Error
薬剤A <- concrete$混和A
薬剤A
6.433
0.874 t value
7.361 Pr(>|t|)
0.000154 ***
(Intercept)
39.7278
1.0076 39.426 1.78e-08 ***
--薬剤B <- concrete$混和B
薬剤A
0.6171
4.56e-05
Signif. codes:6.4333
0 ‘***’ 0.001
‘**’ 10.426
0.01 ‘*’
0.05 ‘.’***
0.1 ‘ ’
#単回帰分析の結果をkekka1に入れ,その概要を出力
薬剤B
1.7500
0.6171
2.836
0.0297
*
1
kekka1 <- lm(強さ~薬剤A)
--summary(kekka1)
Signif.
‘***’ 0.001
‘*’ 0.05
‘.’ 0.1 ‘ ’
Residualcodes:
standard0 error:
2.141‘**’
on 0.01
7 degrees
of freedom
1Multiple R-squared: 0.8856,
#重回帰分析の結果をkekka2に入れ,その概要を出力
Adjusted R-squared: 0.8692
F-statistic: 54.18 on 1 and 7 DF, p-value: 0.0001545
kekka2 <- lm(強さ~薬剤A+薬剤B)
Residual standard error: 1.511 on 6 degrees of freedom
summary(kekka2)
Multiple R-squared: 0.9511,
Adjusted R-squared: 0.9348
9
F-statistic: 58.37 on 2 and 6 DF, p-value: 0.0001168
重共線性の問題
• 多数の説明変数の間に相関がある場合
– 目的変数への効果を一意に分離できない
– 係数の推計値が安定しない(直感に反する符
号を取るなど)
すべての観測値が
ほぼ一直線上にある
Y
この直線を含むような平面で
あれば、どの式を使っても当て
はまりにはほとんど差はない
Z
X
直線上にない場所のYの予測
値には大きな差がでる
10
目的変数をどの程度[記述]出来たか?
Xiによる説明式がない場合
yiの推計値
Y
として、
平均値 y を
y
使うしかない
Yˆi  f ( X i ) なる説明式がある場合
回帰(Xiで説明
できた y から
のずれ)
Y
y
yˆi  y
yˆi
残差・誤差
xi
回帰平方和
平均値周りの
バラツキ
(全平方和)
ST   ( yi  y )
i
X
S R   ( yˆi  y ) 2
i
2
残差平方和
S E   i2   ( yi  yˆi ) 2
i
i
 i  yi  yˆi
決定係数
R2 
SR
S
 1 E
ST
ST
説明できた
平方和の割合
11
決定係数
平均値周りの
バラツキ
(全平方和)
ST   ( yi  y ) 2
i
回帰平方和
S R   ( yˆi  y ) 2
i
残差平方和
S E   i2   ( yi  yˆi ) 2
i
i
S
S
決定係数 R squared
説明変数を増やせ
R2  R  1 E
ST
ST
説明できた平方和の割合
ば1に近づく
本当はあまり意味のない変数を増やしてしまう危険性あり
S E /(n  k )
n 1
~2
自由度調整済み決定係数
R  1
 1
(1  R 2 )
ST /(n  1)
nk
Adjusted R squared
Residual standard error: 1.511 on 6 degrees of freedom
Multiple R-squared: 0.9511, Adjusted R-squared: 0.9348
12
回帰による記述(説明)力の検定
– 回帰平方和が統計的に大きな意味を持っている
か?
– 分散分析表を作り、F検定を行って判断する。
• 帰無仮説:回帰平方和は誤差平方和と同程度の大きさ
(回帰式は、誤差に比べて大きな説明力はない)
• 対立仮説:回帰平方和は誤差平方和より大きい
(回帰式によって誤差よりもかなり大きい部分が説明でき
た)
13
回帰による記述(説明)力の検定例
• 帰無仮説:回帰平方和は誤差平方和と同程度の大きさ
(回帰式は、誤差に比べて大きな説明力はない)→棄却
• 対立仮説:回帰平方和は誤差平方和より大きい
(回帰式により、誤差よりもかなり大きい部分が説明できた)
Multiple R-squared: 0.9511, Adjusted R-squared: 0.9348
F-statistic: 58.37 on 2 and 6 DF, p-value: 0.0001168
14
推測統計学の導入
誤差(残差)が確率分布 (正規分布)をもつと仮定
 i  yi  f ( X i ) ~ N(0, )
2
もともとの観測値は、真の値を中心に確率的に
ばらついていたと考える
yi  f ( X i )   i
Y
Y
yi
yˆi
Yˆi  f ( X i )
残差・誤差
 i  yi  yˆi
xi
X
X
関数形 f (X ) を決めると、誤差の実現値 1 ,,  n が決まる。
それが正規分布のもとで発生する確率をチェックできる。 15
最尤推定法としての最小二乗法
誤差(残差)が確率分布 (正規分布)をもつと仮定
 i  yi  f ( X i ) ~ N(0, )
2
残差の実現値の同時発生確率(正規分布確率密
度の積)が最大になるように、関数 f (X ) を決める。
Y
L
1
e
( i /  ) 2
2
n
n
2
L  log L   log(2 )  n log    i /  2
2
i 1
i 1, n
2
この最大化のために最終項の残差二乗和
n
n
2
S E    i  ( yi  f ( X i ))2 を最小化するように、
X
i 1
i 1
関数 f (X ) を定める必要がある
(記述統計学で使われてきた)最小二乗法は、
(推測統計学の考え方に従う)最尤推定法でもある。
16
AIC 赤池の情報量基準
• AIC =-2(最大対数尤度)+2パラメータ数
• パラメータ数=k+1(回帰係数と残差の分散)
• 最大対数尤度=  n log( 2ˆ 2 )  n
2
2
であるから,
2
ˆ
AIC  n log(2 )  n  2(k  1)
この値が小さいほうが,無駄の
ない効率的なモデル式ができ
ていることを意味している.
> AIC(kekka1)
[1] 42.98059
> AIC(kekka2)
[1] 37.32718
17
パラメータ推定量の確率分布
• サンプルの観測値は、真の回帰式から誤差の分
だけ、上下にずれている。
• 誤差を含む観測値から計算したパラメータの推
定値も、真の値からずれている。
Y
Y
切片の
推定値
真の
切片は0
傾きの推定値
真の
傾きは0
X
X
各サンプルの誤差が平均0分散 2 の正規分布に従うこと
を用い、パラメータの推定量 ˆk が従う分布を計算すれば、
18
真の値からかけ離れている確率をチェックできる。
最小二乗推定量の確率分布
• 最小二乗推定量の行列表示 ˆ  ( X t X )1 X t y
y  X  
• 観測値の誤差構造
– これらより、
ˆ  ( X t X ) 1 X t ( X   )
 ( X t X ) 1 ( X t X )   ( X t X ) 1 X t 
   ( X t X ) 1 X t 
– ここで、 ( X t X )1 X t は確率変数ではない定数。
• 不偏性
– 誤差は正と負を同程度にとるから、E[ε]=0
E[ˆ ]    ( X t X )1 X t E[ ]    ( X t X )1 X t 0  
– 多数のサンプルを取りながら繰り返して推定す
ると、その平均値は真の値に一致
19
最小二乗推定量の確率分布
(2)
• 最小二乗推定量の確率分布 ˆ    ( X t X )1 X t
• 観測値の誤差が、分散  2 の正規分布に従うなら、
それに定数項をかけて足したものは正規分布に従う。
– パラメータの分散共分散行列は、
E[(ˆ   )(ˆ   )t ]
 E[( X t X ) 1 X t  t X ( X t X ) 1 ]
  2 ( X t X ) 1 X t X ( X t X ) 1   2 ( X t X ) 1
• j番目のパラメータの分散:  V jj j番目の対角要素
• 誤差の母分散 2 はわからないから、残差平方和
2
から求めた分散で代用
ˆ
(
y

y
)
SE
 i i
2
2
ˆ 
nk

nk
20
•
最小二乗推定量の確率分布
(2)
正規分布に従う変数の「ずれ」の評価
– 母分散がわからないので、
サンプル(残差平方和)からの
推定値を用いる→t分布
• t統計量の作成
ˆ 2 
tj 
SE

nk
2
ˆ
(
y

y
)
 i i
nk
ˆ j   j
– 自由度n-kのt分布に従う
ˆ 2V jj
– 通常、帰無仮説として  j  0
(j番目の要因の影響はない)とおけば、
ˆ 2V jj を用いて検定できる。
検定統計量 t j  ˆ j
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 39.7278
1.0076 39.426 1.78e-08 ***
薬剤A
6.4333
0.6171 10.426 4.56e-05 ***
薬剤B
1.7500
0.6171 2.836 0.0297 *
21
回帰式のチェック
plot(強さ,薬剤A)
plot(kekka2$fitted.values,強さ,asp=1)
あるいは,
plot(kekka2)によって,4種類の結果診断グラ
フが出力される.
22
AICに基づくステップワイズ法に
よる変数選択
• すべての変数を含むモデルをlm()で推定
• その結果が入ったオブジェクトに,step()関
数を適用する > step(kekka2)
Start: AIC=9.79
強さ ~ 薬剤A + 薬剤B
Df Sum of Sq
RSS
AIC
<none>
13.707 9.786
- 薬剤B 1
18.375 32.082 15.440
- 薬剤A 1 248.327 262.034 34.341
Call:
lm(formula = 強さ ~ 薬剤A + 薬剤B)
Coefficients:
(Intercept)
薬剤A
薬剤B
39.728
6.433
1.750
23
非線形式の回帰分析
• 非線形の関係式でも、置き換えによって線
形化できる場合がある
y  1e2 x
log y  log 1   2 x
y  1  2 x  3 x
2
z  0   2 x
y  1  2 x  3w
24