回帰分析 - 東北大学病院

2015/10/8 第2回
医学統計勉強会
東北大学病院循環器内科・東北大学病院臨床研究推進センター共催
東北大学大学院医学系研究科EBM開発学寄附講座
宮田敏
回帰分析
回帰分析 (regression analysis) は，一つの連続数（実
数）の値を複数の変数によって説明，予測する多変量
解析 (multivariate analysis) の一つ．
Y   0  1 x1     k xk  
Y : response variable, 従属変数，被説明変数
x1 ,  , xk : independent variable, 独立変数，説明変数
 0 , 1 , ,  k : regression coefficient, 回帰係数
 : error term, 撹乱項、誤差項
2015/10/8
東北大学医学統計勉強会
2
Example 1: 親子の身長
68
child: 子の身長
parent: 両親の身長の平均
62
64
66
child
70
72
74
205組の夫婦から生まれた，
928人の成人した子供の身
長（インチ）．
64
66
68
70
72
74
データは0.1インチ刻みに
丸められている．そのため
データ点が重なっている．
parent
Galton, F. (1886). Regression Towards Mediocrity in Hereditary Stature Journal of the Anthropological Institute, 15, 246‐263 2015/10/8
東北大学医学統計勉強会
3
データの要約
データが手に入ったら，まずデータを要約して，その
傾向，特徴を把握する．
個々の変数の数量的要約：Galtonの親子の身長
Min.
parent 64
child
61.7
1st Qu. Median Mean
67.50
68.5
68.31
66.20
68.2
68.09
3rd Qu. Max.
69.5
73
70.2 73.7
SD
IQR
1.79
2
2.52
4
• 平均，中央値はほぼ等しく，分布の位置は同じ．
• 標準偏差（SD），四分位点間距離（IQR）は child の方が大きく，子の分布の散らばりの方が大きい．
2015/10/8
東北大学医学統計勉強会
4
個々の変数の視覚的要約：Galtonの親子の身長
74
200
72
100
62
64
66
68
70
72
74
72
74
66
68
parent
0
62
50
64
150
Histogram of child
Frequency
70
0
Frequency
Histogram of parent
62
64
66
68
70
parent
child
child
グラフを比較するときは，軸を揃えることがコツ．
2015/10/8
東北大学医学統計勉強会
5
標本共分散
二つの変数 x と y が互いに影響し合っているとき，
x と y が如何に強く関係しあっているか知りたい．
このとき x と yの共分散(covariance)を以下で定義する．
1 n
Cov X , Y   i 1 ( X i  X )(Yi  Y )
n
x と yの間に正の(負の)相関があるとき、Cov(X, Y)は
それぞれ正(負)になる．
2015/10/8
東北大学医学統計勉強会
6
正の相関の場合
第１象限
3
第２象限
xi  x & yi  y
2
  xi  x   0 &  yi  y   0
xi  x & yi  y
-1
  xi  x  yi  y   0
-2
y
y
0
1
  xi  x  yi  y   0
xi  x  yi  y   0
-3
  xi  x   0 &  yi  y   0
正の相関の場合，
データの大部分は
第1, 3象限にある。
第３象限
共分散
-2
x
-1
0
1
第４象限
1
2
n
x
n
 x  x  y
i 1
i
i
 y
は正になる。
東北大学医学統計勉強会
2015/10/8
7
標本相関係数
二つの変数x と yの間の相関係数 (correlation coefficient) を，以下で定義する．相関係数はx と yの線形関係の強
さを測る量である．相関係数は単位に依存しない．
Corr  X , Y   r 

n
i 1
( X i  X )(Yi  Y )
2
(
X

X
)
i 1 i
n
2
(
Y

Y
)
i 1 i
n
•相関係数r はx と yの線形関係の強さを測る．
•-1 ≦ r ≦ 1
•r = +(-) 1 : 正(負)の完全な相関，線形関係
2015/10/8
東北大学医学統計勉強会
8
回帰分析
二つの変数x と yの関係が線形（直線）で近似できると
する．このときx と yの関係を以下の回帰式 (regression equation) でモデル化する．
Y   0  1 x  
Y : response variable, 従属変数，被説明変数
x : independent variable, 独立変数，説明変数
 0 , 1: regression coefficient, 回帰係数
 : error term, 撹乱項、誤差項
2015/10/8
東北大学医学統計勉強会
9
被説明変数yの観測値
yi   0  1 xi   i
（観察できる）
xiにおける真の回帰関数の値
yi   0  1 xi
撹乱項  i
（観察できない）
真の回帰関数
y   0  1 x
（観察できない）
（観察できない）
xi
（観察できる）
観察可能な xi , yi  から、y   0  1 x を推定したい．
2015/10/8
東北大学医学統計勉強会
10
回帰分析のモデルの仮定
• 線形性 (Linearity)：Y   0  1 x  
被説明変数yと説明変数xの関係は直線で近似できる．
xi , Yi i 1 は互いに独立である．
• 独立性 (Independence)：
あるサンプルの値が他のサンプルの値に影響しない．
n


• 正規性 (Normality)： i ~ N 0,  2 , iid
撹乱項  i は正規分布に従う．
• 等分散性 (homoskedasticity)：σ2. 分散一定
東北大学医学統計勉強会
2015/10/8
11
正規分布：Normal distribution
f ( x;  ,  2 ) 
•


1
exp  ( x   ) 2 /(2 2 ) ,  x  
2 
E  X    ,V  X    2 .
• 分布の形状は  と  2 によって特徴づけ
(parameterized) られる。
• 標準正規分布 (standard normal distribution)
f ( z;0,1) 
1 z2 / 2
e
,  z  ,
2
z
 ( z )  P ( Z  z )   f ( y;0,1)dy.

• 釣り鐘型 (bell-shaped) で、左右対称な分布．
2015/10/8
東北大学医学統計勉強会
12
0.0
0.1
0.2
0.3
0.4
Standerd Normal Distribution
-3
-2
-1
0
1
2
3
東北大学医学統計勉強会
2015/10/8
13
回帰係数の推定
最小二乗法 (ordinary least squares estimation, OLSE)
min    y  
1
i 1
i
 1 xi 
2
n
0,
0
n
 
2
  i 1 yi   0  1 xi   0 n 0  1 n xi  n yi
 0

i 1
i 1


n
n
n
2
n
2
 
 0 i 1 xi  1 i 1 xi  i 1 xi yi




y


x
0



0
1 i
 1 i 1 i
ˆ1  n  xi  x  yi  y  n  xi  x 2
i 1
i 1

ˆ0  y  ˆ1 x



 

推定された回帰直線： yˆ  ˆ0  ˆ1 x
2015/10/8
東北大学医学統計勉強会
14
最小二乗推定量の性質
 
 
 x
 varˆ  
n ( x  x )
 E ˆ0   0 , E ˆ1  1 , 不偏性

2
2
ˆ
0
2
i
0
2
,
i

2
ˆ
1
 
 var ˆ1 

 ˆ0 ~ N  0 ,  2ˆ
0
2
,
( xi  x )
, ˆ1 ~ N 1 ,  2ˆ ,
2


1

正規性
( yi  yˆ i ) 2 SSE
 s  ˆ 

, E ˆ 2   2 ,
n2
n2
2
 
2
東北大学医学統計勉強会
2015/10/8
15
 信頼区間 (Confidence Interval) :
s
ˆ  1
~ t n  2 , sˆ 
T 1
1
2
sˆ
(
)
x
x


i
1


 C.I. of 1 : ˆ1  sˆ t / 2,n  2, ˆ1  sˆ t / 2,n  2, ,

1
1
同様に, C.I. of  0 : ˆ0  sˆ t / 2,n  2, ˆ0  sˆ t / 2,n  2,
0
0

(Hypothesis Testing) :
 仮説検定　H 0 : 1  10 vs. H a : 1  (or  or ) 10
ˆ1  10
検定統計量 : T 
~ t n  2 under H 0
sˆ
1
Ha : β1 > β10
Ha : β1 < β10
Ha : β1 ≠β10
2015/10/8
t > tα,n-2
t < -tα,,n-2
|t| > tα/2,n-2
東北大学医学統計勉強会
16
回帰分析の結果
•
•
•
•
•
回帰係数の推定値 ˆ0  0.9183, ˆ1  0.6904
回帰係数の有意性検定のp値
決定係数（被説明変数の変動のうち回帰によって説明され
た変動の割合） yの変動の40.19%が説明された
Model utility test（回帰モデル全体の有意性検定。後でも
う一度触れます）のp値 p=0.049
撹乱項の標準誤差 s=0.6744
検定の p‐value
H 0 :  j  0 vs. H1 :  j  0
回帰係数の推定値
ˆ  0.9183
̂ j の標準誤差
（CIに使う）
0
ˆ1  0.6904
決定係数
2015/10/8
s  ˆ
Model utility test p‐value
東北大学医学統計勉強会
17
多変量回帰分析

Y   0  1 x1     k xk   ,  ~ N 0,  2

Y : response variable, 従属変数，被説明変数
x1 ,  , xk : independent variable, 独立変数，説明変数
 0 , 1 , ,  k : regression coefficient, 回帰係数
 : error term, 撹乱項、誤差項
モデルの仮定
1. パラメターに関する線形性 (Linearity)
2. 撹乱項の独立性 (Independence)
3. 撹乱項の正規性 (Normality)
4. 撹乱項の等分散性 (homoscedasticity) 2015/10/8
東北大学医学統計勉強会
18
最小二乗法によるパラメターの推定
  0 , 1 ,  ,  k は最小二乗法で推定される:
min
 0 ,,  k
 y  
n
j 1
j
0  1 x1 j     k xkj 
2
ˆk 1   X ' X 1 X ' Y ,
where X  [ xij ]n( k 1) , Y  (Y1 ,  , Yn )'
 E ( ˆ j )   j , j  1,  , k : unbiased.
1
 V ( ˆ )   2  X ' X  .
1
 ˆ ~ N  ,  2  X ' X 


東北大学医学統計勉強会
2015/10/8
19
パラメターの推測、信頼区間、検定
( yi  yˆ i ) 2
SSE
.E ( s 2 )   2 .
 s  ˆ 

n  (k  1) n  (k  1)
ˆ j   j
1/ 2
~ t n ( k 1) , sˆ  [ s 2 ( X ' X ) 1 ] jj , j  1,  , k .
T 
j
sˆ
2
2


j
 C.I. of  j : ( ˆ j  sˆ t / 2,n ( k 1) , ˆ j  sˆ t / 2,n ( k 1) )
j
j
 Hypothesis testing : H 0 :  j   j 0 vs. H a :  j   j 0 ,
Test statistic : T 
2015/10/8
ˆ j   j 0
sˆ
~ t n ( k 1) under H 0 .
j
東北大学医学統計勉強会
20
決定係数 (coefficient of determination)
回帰係数の有意性検定は個々の係数の検定．回帰モデ
ル全体のパフォーマンスを測る方法が必要．
SST (Total Sum of Squares):
 y  y
yの全変動
  y  yˆ 
SSE (Error Sum of Squares): 回帰で説明されなかった変動
SSR (Regression sum of Squares):   yˆ  y 
回帰で説明された変動
2
n
i 1
i
2
n
i 1
i
2
n
i 1
R2=決定係数=
2015/10/8
SSR
SST
i
i
=回帰によって説明された変動の割合
東北大学医学統計勉強会
21
Model utility test
回帰モデル全体の有意性を検定するために，以下の
Model Utility Test が知られている．
H 0 : 1     k  0 vs. H1 : not H 0
検定統計量：
R2 k
SSR k

~ Fk ,n k 1 under H 0
F
2
1  R  n  k  1 SSE n  k  1
棄却域：Reject H0 if F  F ,k ,n k 1
H0は、Y   0   と同義．回帰が全く無効．
2015/10/8
東北大学医学統計勉強会
22
回帰分析の結果のチェックポイント
•
•
•
•
回帰係数の推定値
回帰係数の有意性検定のp値
回帰係数の信頼区間
決定係数（被説明変数の変動のうち回帰によって説
明された変動の割合）
 決定係数は、いくつ以上なら良いか?
もちろん、R2は大きいに越したことはない。 R2が小さいということは、
yに影響する、今のモデルに含まれない変動要因が存在する、というこ
と。もし、 R2が小さくてもなお有意な説明変数があれば、上記の
limitationの下で、回帰は有効である。
• Model utility test（回帰モデル全体の有意性検
定）のp値
• 撹乱項の標準誤差
2015/10/8
東北大学医学統計勉強会
23
回帰診断
回帰モデルの仮定の確認：
• 線形性: Y と x の間の線形関係． Y と xの間には非線
形な関係がない． x同士の間には線形関係がない．
–Multiple scatter plots.
• 独立性: {(xi, Yi)}i=1n は互いに独立
–residual vs. fitted value plot
• 正規性: ε~N(0, σ2)
–残差のQQ-norm plot （後述する）
• 等分散性: Var(ε)= σ2
–residual vs. fitted value plot
2015/10/8
東北大学医学統計勉強会
24
回帰診断（続き）
独立性，正規性，等分散性の仮定は，いずれも撹乱項
についての仮定．撹乱項そのものは観察できないため，
残差 (residuals) : ei  yi  yˆ i をレプリカとして使う．
residuals
y
0
-6
-4
-4
-2
残差
-2
-1
0
1
2
3
x
2015/10/8
-1
0
1
2
fitted values
分散増大傾向がある
• 残差が均一ならば、
等分散性の仮定は満
たされる．
-2
0
2
2
4
4
Residuals vs. fitted value plot
残差も増大傾向がある
• 独立性の仮定が満た
される場合，残差プ
ロットには特異なパ
ターンがない．
東北大学医学統計勉強会
25
分布の正規性の確認
標本分布の正規性の確認は，適切なモデルを選択する
上で重要．
Definition: n個の標本を大きさ順に並べたとき，i番目に
小さな標本は [100(i-.5)/n] 標本パーセント点 (sample percentile) であるという．
例えば標本が，正規分布など特定の確率分布から抽出
されたとする．このとき，その特定の分布の理論上の
[100(i-.5)/n]パーセント点は，データの[100(i-.5)/n]標本
パーセント点の近くにあるはずである．
2015/10/8
東北大学医学統計勉強会
26
QQ‐norm plot (Normal probability plot)
Definition: n 個の標本が得られたとき，標準正規分布の
[100(i-.5)/n] パーセント点と，i番目に小さな観測値=
[100(i-.5)/n] 標本パーセント点のプロットを，QQ‐norm plotという。
正規分布
• 元のデータが正規分布
から得られた場合，
QQ‐norm plotは直線上
にプロットされる．
非正規分布
Normal Q-Q Plot
10
Sample Quantiles
0
• 元データが正規分布に
従わない場合、直線か
ら外れる．（右図）
-2
5
-1
Sample Quantiles
1
15
2
Normal Q-Q Plot
-2
-1
0
1
-2
2
-1
0
1
2
Theoretical Quantiles
Theoretical Quantiles
東北大学医学統計勉強会
2015/10/8
209のコンピュータのCPUの持つ，
性能と各種特性値．
CPUデータ
1500
0
30000
0
20
40
1000
0 500
27
1500
0 400
perf
15000
0 500
syct
30000
0
mmin
250
0
mmax
40
0
100
cach
'name' Manufacturer and model
'syct' cycle time in nanoseconds
'mmin' minimum main memory in kilobytes
'mmax' maximum main memory in kilobytes
'cach' cache size in kilobytes
'chmin' minimum number of channels
'chmax' maximum number of channels
'perf‘ published performance on a benchmark mix relative to an IBM 370/158‐3
chmax
0 400
1000
0
15000
0
100
250
0 50
0 50
150
0
20
chmin
150
P. Ein‐Dor and J. Feldmesser (1987) Attributes of the performance of central processing units: a relative performance prediction model. Comm. ACM. 30, 308–317.
2015/10/8
東北大学医学統計勉強会
28
CPUデータ（続き）
1500
0
30000
0
20
40
1000
0 500
1500
0 400
perf
15000
0 500
syct
30000
0
mmin
• yとxの間に，非線形な関係が存在しな
いか？
– perf と syct の間に，明らかな非
線形関係がある．
• 誤差項の分散は一定か？
– mmin 等に，明らかな分差増大傾
向がある．
250
0
mmax
40
0
100
cach
• x同士の間に，線形な関係が存在しな
いか？
– 例えば，mmin と mmax の間に
明らかに線形関係がある．
chmax
0 400
1000
0
15000
0
100
250
0 50
0 50
150
0
20
chmin
150
P. Ein‐Dor and J. Feldmesser (1987) Attributes of the performance of central processing units: a relative performance prediction model. Comm. ACM. 30, 308–317.
2015/10/8
予備的な視覚的要約の段階でも，線形回
帰モデルを当てはめるのは不適切である
ことがわかる．
でも，ともかく回帰モデルを当てはめて
みる．
東北大学医学統計勉強会
29
CPUデータ（元データの回帰分析）
100
0
Sample Quantiles
200
300
400
Normal Q-Q Plot
-200
-100
非正規性
-3
-2
-1
0
1
2
3
Theoretical Quantiles
100
非線形性
0
• 決定係数：R2=0.8649 被説明変数の変動の
86.49％が説明できた．
residuals
200
300
400
• 説明変数の有意性検定は，chmin を除き，ほ
とんどが強く有意．
-100
分散増大
-200
• Model utility test: p-value < 2.2×10-16
0
回帰分析は，成功しているとしかいいようがない．
回帰診断の結果，モデルの仮定は破綻している。
2015/10/8
東北大学医学統計勉強会
200
400
600
800
fitted values
30
変数変換
線形回帰モデルの仮定（線形性, 正規性, 等分散性）
が満たされないとき，変数に何らかの変換を施すこ
とで，モデルを改善できる場合がある．
例えば，撹乱項の分散が説明変数の値とともに大き
くなる場合，logarithmic/power 変換が有効であること
が多い．
被説明変数の予測値を得るには，まず変換された変
数に対して線形回帰モデルを当てはめ，次にもとの
モデルに逆変換する．最もよい変換を選ぶため，い
くつかの変換を試してみる必要がある．
2015/10/8
東北大学医学統計勉強会
31
Box‐Cox変換
対数変換，冪変換を組み合わせたBox‐Cox変換により，
分散の安定化と正規性の改善を同時に達成できる場
合がある．
( y   1)  :   0
Box‐Cox変換： y (  )  

log( y ) :
 0
Box‐Cox変換は，パラメターλによって特徴付けられる．
パラメターλは，モデルの適合度を最適化するように，
ソフトウエアにより自動的に選択される．
（統計解析ソフトRなどが， Box‐Cox変換を実装している）
2015/10/8
東北大学医学統計勉強会
32
CPUデータ（Box‐Cox変換後）
2
0
-4
-2
Sample Quantiles
4
6
Normal Q-Q Plot
-3
-2
-1
0
1
2
3
6
Theoretical Quantiles
residuals
0
-4
• Model utility test: p-value < 2.2×10-16
-2
• 決定係数：R2=0.8821 被説明変数の変動の
88.21 (> 86.49) ％が説明できた．
2
4
• chmin，chmax は有意ではない．
回帰分析は，成功している．
2015/10/8
5
10
15
20
25
fitted values
東北大学医学統計勉強会
33
変数選択
多変量回帰モデルを考えるとき，多数の説明変数の
候補が考えられる場合がある．その時，被説明変数
の変動を説明する，説明変数の最適な組み合わせを
探索する必要がある．
• モデル選択基準（赤池の情報量基準 (AIC), ベイズ
情報量基準 (BIC), etc.）を最小にする，説明変数の
組み合わせを求める．
• 当初のモデルに新たな変数を投入もしくは除去し
たとき，モデルの適合度の差が有意なものか検定
する。（ワルド検定，尤度比検定，ほか）
2015/10/8
東北大学医学統計勉強会
34
ステップワイズ法
説明変数の候補がk個あるとき，可能な変数の組み合
わせは2k通り．総当たりの探索は無理．
• 変数増加法 (forward selection)：定数項のみのモデ
ルから出発して，有用な変数を加えていく方法．
• 変数減少法 (backward elimination)：すべての変数
を用いたモデルから出発して、不要な変数を除去
していく方法．
• 変数増減法 (stepwise procedure)：増加法と減少法
の組み合わせ．
2015/10/8
東北大学医学統計勉強会
35
変数選択の実際（私の場合）
医学統計の性質上，説明変数は取りこぼしなく広め
に選択したい．しかし，候補の数がサンプル数より
大きければ，最初に全ての変数を用いたモデルは推
定できない．また，欠測が多いとサンプル数が減る．
1. まず単変量解析で，変数ごとのp値を出す．
2. ある程度p値が小さい変数に，候補を絞る．
（p < 0.2 程度で，あまりきつく絞らない）
3. 絞った候補から，変数減少法で選択する．（割合
多めの変数が選択される）Rでは stepAIC()コマン
ドを用いる（MASSパッケージが必要）．
2015/10/8
東北大学医学統計勉強会
36
変数選択の問題
• 変数選択の結果選ばれたモデルに，有意でない変
数が含まれる．
変数選択は，被説明変数の変動をもっともよく説明
する変数の組み合わせを探索する．個々の変数が有
意でなくても，組み合わせ全体として最適であると
解釈する．
• 変数選択の結果，興味のある変数がモデルから除
かれてしまった．
上記に矛盾するようだが，回帰の目的は被説明変数
を説明することだけではない．変数間の関係を推測
するため，興味ある変数を強制投入してもよい．
東北大学医学統計勉強会
2015/10/8
37
Take Home Message
1．回帰分析
2．共分散と相関係数
3．線形回帰モデル
• 回帰係数の推定．最小二乗推定量の性質
4．回帰診断：回帰モデルの仮定の確認
• 散布図：線形性の確認
• QQ‐normプロット：残差の正規性の確認
• 残差プロット：等分散性，独立性の確認
5．Box‐Cox変換：分散の安定化と正規性の向上
6．変数選択
2015/10/8
東北大学医学統計勉強会
38

Download Report