τ 推定量に基づくロバスト・リッジ回帰の研究
M2012MM045 塚原一翔
指導教員：木村美善

はじめに

12.2

回帰モデルにおいて, 最小 2 乗推定量は標準的仮定の下では望ましい推定量である. しかし, 外れ値や多重共線性が存在する場合にはその値に対して敏感であり, 推定量の正確性を失ってしまうことが知られている. 外れ値が存在するとき, 外れ値に対する影響を受けにくいロバスト回帰を用いることが望ましい. また, 多重共線性の問題に対しては, リッジ回帰が広く用いられている手法の一つである. しかし, このリッジ回帰は, 外れ値に対して有効に対処できず, 外れ値の影響を受けやすいという欠点がある.
また, 実際の分析に用いられるデータには外れ値と多重共線性が混在している場合が多くある. このような場合には, ロバスト回帰とリッジ回帰を組み合わせたロバスト・リッジ回帰という手法を用いることで, 外れ値と多重共線性という 2 つの問題に対して同時に対処することが可能となる. 本研究の目的は，ロバスト・リッジ回帰の理論と推定量が持つ性質を整理し, 様々なロバスト推定量に基づ
くロバスト・リッジ回帰推定量, 特に τ 推定量に基づくロバスト・リッジ回帰推定量に着目し, その有効性について調べることである.

2回帰分析

2.1

線形回帰モデル

目的変数を y, p 個の説明変数を x1 , x2 , · · ·, xp , 回帰係数を β0 , β1 , · · · , βp としたとき, 次のような線形回帰モ
デル

yi = β0 + β1 xi1 + · · · + βp xip + εi ,i = 1, · · · , n (1)を考える. ここで, ε は近似によって生じる誤差を表す. このモデルを行列で表記すると

...
...
...

x1p
x2p .. ,. xnp

ε1
 ε2 
ε= . , .. εnである.

最小 2 乗法は, (1) 式のモデルにおける残差平方和 (RSS)RSS = ||ε||2 = (y − Xβ)′ (y − Xβ)(2)を最小にするような β を求める手法であり, 回帰分析手法の中で最も基本的かつ最も広く用いられているものである. (2) 式は β の 2 次関数になっており, これを β で偏微分したものを 0 とすることにより, 最小 2 乗推定量はβ̂ LS = (X ′ X)−1 X ′ Y(3)となる.
モデルにおける誤差が等分散性, 不偏性, 無相関性の仮定を満たすとき, 最小 2 乗法推定量はすべての線形不偏推
定量のなかで最も小さい分散をもつ最良線形不偏推定量となり, さらに正規分布が仮定されている場合には最良不偏推定量となる. しかし, 実際のデータは, 外れ値や多重共線性が存在したりするなど, これらの仮定をすべて満たすような場合は稀である.

3リッジ回帰推定量3.1

多重共線性とは重回帰分析において説明変数の間に強い相関関係が存在する場合, これらの説明変数の間には多重共線性があるという. この場合, 回帰分析により得られる結果に悪い影響が出ることがある. 具体的には, 同時に用いる説明変数の増減により回帰式の係数が大きく変化したり, 決定係数が高い一方で t 値が低く, 有効な推定結果が得られなかったり, 通常考えられる符号と異なる結果が得られる, などの症状が生じる.

3.2

リッジ回帰推定量

X ′ X の固有値を λ1 ≥ λ2 ≥ · · · ≥ λp+1 , 回帰係数 β の最小 2 乗推定量 (LS 推定量)β̂ LS は不偏推定量であるので, その平均 2 乗誤差は

y = Xβ + εとなる. ここで
1 x11 1 x21
X=
..
 ..
.
 .
1 xn1

最小 2 乗法y1 y2 
y= . , .. ynβ0
 β1 
β= .  .. βp

M SE[β̂] = E[(β̂ − β)′ (β̂ − β)]= σ2p+1∑λ−1i(4)i=1となる. データに多重共線性が存在する場合に, X ′ X の固有値 λ には 0 に極めて近いものが存在するため, (4) のLS 推定量の平均 2 乗誤差は大きく発散してしまう可能性がある. Hoerl and Kennard(1970) はこのようなとき
でも平均 2 乗誤差を小さく抑えるための手法としてリッジ回帰推定量を提案した. リッジ回帰推定量はモデルにリッジ・パラメータとよばれる定数 k ≥ 0 を取り入れることで LS 推定量 β̂ LS を縮小させたものであり
LSβ̂(k)= (X ′ X + kI)−1 X ′ Y(5)LSにより定義される. 特に k = 0 のとき, β̂(k)は最小 2 乗推定量に等しくなる. しかし, リッジ回帰推定量は不偏推定量ではない. したがって k の増加に伴い, 偏りも大きくなってしまう.

4ロバスト回帰推定量

4.1

M 推定量

4.5

τ 推定量

τ 推定量は Yohai と Zamar(1988) により提案されたロバスト回帰推定量であり,β̂ τ = arg minp τ (β)β∈Rによって定義される. ここで,τ (β) はM 推定量は,Huber (1964) によって提案されたロバスト推定量である. ロバスト推定量の中でも最も一般的なものであり, 微分可能な偶関数 ρ を用いて
β̂ M = arg minβn∑i=1(6)と定義される. 関数 ρ はこれまでに様々なものが提案されているが, Huber (1964) によるものが最も一般的である. また, (6) 式からもわかるように, ρ(ri ) = ri2 とすると, これは最小 2 乗推定量に等しいことがわかる.

LM S 推定量β̂

= arg min med{r1 (β) , · · · , rn (β) }2β2(7)として定義される.

4.3

LT S 推定量(12)

多重共線性の存在する人工データ作成法

5

金・田中 (1993) による多重共線性の存在するデータの作成手順は次の通りである.

5.1

手順

1. 変数の数 (p) と標本の大きさ (n) を固定する.
(a) 線形独立な p 次元ベクトル {ei }p1 を生成する.
(b) {ei }p1 をグラム・シュミッドの直交化法を用いて, ベクトルのノルムが 1 であるような正規直交ベクトル {vi }p1 に変換し, それを直交行列 Vにする.
3. 対角行列 Dp×p を作る:

LT S 推定量は Rousseeuw (1984) によって提案された手法であり, 残差平方を昇順に並び替えた順序統計量のm 番目までの和を最小にする
β̂ LT S = arg minβm∑2r(i)(β)(8)(a) condition index κ1 , κ2 , · · · , κp と分散の和 c(=∑pj=1 λj ) を指定する. 指定された conditionindex と分散の和 c に基づき, 固有値 λi =∑pc/(κi j=1 κ−1j ) を計算する.
1/2i=122として定義される. ここで r(1)(β) ≤ r(2)(β) ≤ · · · ≤2r(n) (β) である. 破綻点は ([n/2]-p + 2)/n であり, n → 0のとき 1/2 となる. LT S 推定量は y 方向だけでなく X 方向に対してもロバストであるが, 漸近効率は高くない.

4.4()
n1 ∑ri (β)ρ2nb2 i=1s(β)

2. 直交行列 Vp×p を作る:

LM S 推定量は, Rousseeuw(1984) により提案されたロバスト回帰推定量であり,LM S

τ 2 (β) = s2 (β)

であり, 尺度 s(β) は (10) により定義されるものである.
ρ2 は ρ と同じ条件を満たす関数である.

ρ(ri (β)),

ri (β) = yi − (β0 + βi xi1 + · · · + βp xip )4.2(11)

S 推定量(b) 求めた各 λiを作る.

を対角要素にする対角行列 Dp×p

4. 行列 Un×p を作る:
(a) N (0, I) に従う p 変量正規乱数 {yi }n1 を発生する.
(b) {yi }n1 の平均ベクトル ȳ と分散行列 S を計算
する.

M 推定量の柔軟性と漸近的性質の良さを保持しながら, 高い破綻点をもち, LM S 推定量や LT S 推定量よりも高い漸近効率をもつことを狙ったのが Rousseeuw とYohai(1984) による S 推定量である. S 推定量は(d) 各 yi を次のように変換する.

(β̂ s ) = arg min s(β)

zi = G− 2 Q′ (yi − ȳ),β(9)n0≤b≤11i = 1, 2, · · · , n (13)(e) 各 zi′ を行とする Un×p を作る.

を満たすものとして定義される. ここで s(β) は1 ∑ ri (β)) = b,ρ(n i=1 s(β)(c) S のスペクトル分解 S=QGQ′ を行う.

(10)を満たすものである. ρ は (−∞, ∞) 上の有界な関数であり, 原点対称, (0, ∞) 上で非減少, ρ(0)=0, b はある定数である. すなわち尺度 s(β) を推定した後, この s(β) を最小にする β を推定量とするものである.

5. データ Xn×p を作る:
上 で 求 め た 行 列 V ,D,U を 用 い て, 行 列 の 積
U DV ′ = X を計算し, 人工データを作る.

シミュレーション

620 組のデータ

ここでは多重共線性のデータを用いて, 次の 3 通りの場合についてシミュレーションを行い, ロバスト・リッジ回帰の有効性を調べる.

(yi∗ , xi1 , xi2 , xi3 , xi4 , xi5 ),にモデル (15) を当てはめ, 回帰係数 β=(β1 , β2 , β3 ,β4 , β5 )′ の推定量 β̂ を求め, (β̂ − 1)′ (β̂ − 1) を計算
する.

1. 外れ値がない場合
2. X 方向に外れ値がある場合• この一連の作業を 30 回繰り返し, 第 i 回目で得られる β̂ を
β̂j , j = 1, 2, ..., 30

3. X と y 方向の両方に外れ値がある場合
6.1

シミュレーションの手順

線形回帰モデルとして

ε ∼ N (0, 1)(14)をモデルとして考え, 回帰係数の真値を β0 = β1 = β2 =β3 = β4 = 1 とする.
• N (0,1) に従い, 多重共線性をもつ 20 組のデータ

301 ∑\MSE(β̂) ={(β̂j − 1)′ (β̂j − 1)}
30 j=1y = β0 + β1 x1 + β2 x2 + β3 x3 + β4 x4 + ε,(xi1 , xi2 , xi3 , xi4 ),

i = 1, 2, · · · , 20

i = 1, 2, · · · , 20を求める.

• LS 推定量, LM S 推定量, S 推定量, τ 推定量について k=0, k=0.01, k=0.05 の 3 つの場合を求める.
6.2を金・田中 (1993) の方法により作成し,シミュレーションによる M SE(β̂) を下の表に示す. 値が小さいほど精度良く推定できていることになる.

yi = 1 + xi1 + xi2 + xi3 + xi4 + εiとする.

表 1 k=0 の場合のシミュレーション結果• 20 組のデータ(yi , xi1 , xi2 , xi3 , xi4 ),

シミュレーション結果

i = 1, 2, · · · , 20にモデル (14) を当てはめ, 回帰係数 β=(β1 , β2 , β3 ,β4 )′ の推定量 β̂ を求め, (β̂ − 1)′ (β̂ − 1) を計算する.

LS 推定
LM S 推定
S 推定
τ 推定

外れ値なし
31.59
287.20531.9047.27

X 方向
114.86
80.02
38.50
49.12

両方
6.63
6.49
6.35
5.70

• X 方向への外れ値としてx5 ∼ (1 − η)N (0, 1) + ηN (8, 3)を考え (η =0.2), ỹi (i=1, 2,…,20) を次のように定義する.

ỹi = 1 + xi1 + xi2 + xi3 + xi4 + x̃i5 + εi .
ここで,xi5 が外れ値でないとき x˜i5 =xi5 , xi5 が外れ値のとき x̃i5 =0 とする. 20 組のデータ(y˜i , xi1 , xi2 , xi3 , xi4 , xi5 ),表 2 k=0.01 の場合のシミュレーション結果

LS 推定
LM S 推定
S 推定
τ 推定

外れ値なし
3.83
7.32
6.15
5.44

X 方向
4.90
7.84
6.18
7.83

両方
3.56
3.81
3.63
3.78

i = 1, 2, · · · , 20
表 3 k=0.05 の場合のシミュレーション結果にモデル

y = β0 +β1 x1 +β2 x2 +β3 x3 +β4 x4 +β5 x5 +ε, ε ∼ N (0, 1) LS 推定(15)
LM S 推定′を当てはめ, 回帰係数 β=(β1 , β2 , β3 , β4 , β5 ) の推
S 推定
定量 β̂ を求め, (β̂ − 1)′ (β̂ − 1) を計算する.
τ 推定

外れ値なし
3.43
5.02
4.42
4.48

X 方向
4.54
5.93
4.94
6.03

両方
3.56
3.74
3.593.66• y 方向への外れ値として ε∗ をε∗ ∼ (1 − η)N (0, 1) + ηN (8, 3)6.3を考え (η = 0.2), yi∗ (i=1, 2, · · · , 20) を次のように定義する.

yi∗ = 1 + xi1 + xi2 + xi3 + xi4 + x̃i5 + ε∗i(16)

外れ値なしの場合の結果と考察

分析結果より, 外れ値なしの場合は k の値に関わらず,>LS 推定が一番精度が良い結果となった. 特に, k=0 のときに, 他の推定量との差が非常に大きくなっている. k の値が大きくなるにつれて他の推定量との差はあまりなく

−0.5−1.00.0−2.0−1.5

Ridge estimates

0.00.51.0なっている. また, 図 1 のように LS 推定では各変数の安定する k の値が非常に小さくなっていることが多かった.
これより, 外れ値が存在しない場合には, ロバスト・リッ
ジ回帰推定量よりリッジ回帰推定量のほうが適しているということがわかった.

0.01

0.02

0.03

0.04

0.05

−1.5

−1.0

Ridge parameter (k)図 3 τ 推定のリッジ・トレース図−2.5−2.0

Ridge estimates−0.50.00

−3.0

6.6
0.00

0.01

0.02

0.030.04

図 1 LS 推定のリッジ・トレース図

6.4

X 方向に外れ値がある場合の結果と考察

20−4−2

Ridge estimates

46X 方向へ外れ値を入れた場合は, 外れ値なしの場合とは逆にロバスト・リッジ推定量のほうが精度が良かった.
ロバスト・リッジ推定量の中でも特に S 推定量が一番良かったが, τ 推定量も同じくらいの精度であった. k=0 のときは 3 つのロバスト・リッジ推定量が LS 推定量より精度が良いが, k=0.01, k=0.05 のときは LS 推定量の各変数の推定量が急速に収束しており, LS 推定量のほうがわずかに精度が良くなっている.

0.00

0.01

0.02

0.03

0.04

0.05

Ridge parameter (k)図 2 S 推定のリッジ・トレース図

6.5

各結果を踏まえたまとめ

0.05

Ridge parameter (k)X と Y 方向の両方に外れ値がある場合の結果と考察

X 方向と y 方向の両方に外れ値を入れた場合は, τ 推定, S 推定, LM S 推定, LS 推定の順番で精度が良かった.
外れ値なしの場合や X 方向に外れ値がある場合とは逆に,ロバスト・リッジ回帰推定量のほうがリッジ回帰推定量より精度の良い結果となった. わずかな差ではあるが, ロバスト・リッジ回帰推定量の中でも特に τ 推定量が一番精度の良い結果となっている. どの分析手法でも, 各変数は推定量が 0 のあたりに収束しているが, X 方向の外れ値の変数のみ, 図の上のほうにある直線のように 0 ではな
く 1 に収束していた. また, この変数は k=0 のときの値からあまり変化がなかった. これより, X 方向と Y 方向の両方に外れ値がある場合は, ロバスト・リッジ回帰推定量, 特に τ 推定量が一番良いということがわかった.

外れ値なしの場合と X 方向へ外れ値がある場合は LS推定量, X と y 方向の両方に外れ値がある場合は τ 推定量が一番良いということがわかった. 各ロバスト・リッジ
推定は, X と y 方向の両方に外れ値があると, LS 推定より精度が良くなっている. これより, 外れ値が両方向へあるような複雑なデータに対しては, ロバスト・リッジ回帰が有効であり, ロバスト・リッジ推定量でも特に τ 推定量が, 多重共線性と外れ値が混在しているデータに対して安定した分析ができることがわかった. k の値が大きくなるにつれて各変数は 0 付近に収束しており, k の値が 0.01前後で安定することが多かった.

7

おわりに外れ値と多重共線性が混在するデータに対して, τ 推定量に基づくロバスト・リッジ回帰はそれらに影響されないことがわかった. 今回は N (0, 1) に従う独立な説明変数を用いた回帰式に, 多重共線性と外れ値の存在する説明変数を用いた回帰式を当てはめたが, k の値をさらに細かく設定して分析したり, 説明変数をさらに増やすなど研究の余地はまだまだあると感じる. また, 試行回数が少なかったり, τ 推定量以外の他の推定量との比較は行っていないため, 他の推定量に対する τ 推定量の優位性についての分析など, τ 推定を含めたロバスト・リッジ回帰分析のさらなる研究が必要である.

参考文献[1] 阿部智成, 暮石一樹, 木村美善 (2013). ロバスト・リッジ回帰推定量について，ACADEMIA Information
Sciences and Engineering Vol.13.
[2] 金鉉彬・田中豊 (1993).”多重共線性を持つ人工データの作成法の一提案”, 日本計算機統計学.
[3] Matias,Salibian.Barrera,Gert,Willems,andRuben,Zamar.(2008),The Fast-τ Estimator forRegression,Journal of Computational and GraphicalStatistics,659-682[4] 武山嵩弘 (2008). ロバスト・リッジ回帰推定量の研究，
南山大学数理情報研究科修士論文.