目次

・線形モデルの前提の再考

・外れ値の存在

・重み付き最小二乗法 (WLS)

・補足：重みの決め方は上記でいいのか？　ほかに方法はないのか？

線形モデルの前提の再考

　これまで様々なパターンの線形モデルを取り扱ってきた。今ここで、検定の時と同じように、線形モデルを使える前提について考えてみよう。線形モデルは、（1）残差の分布が正規分布に従う、（2）残差の分散は説明変数の値によらず一定、（3）説明変数と被説明変数が直線関係（線形）であることを仮定している。これらの仮定が重要であることは、線形回帰にて最小二乗法を使って係数の推定を行っていれば、おのずと理解できる。ゆえに本来は、これらのいずれかの仮定が満たされないときに単純な線形モデルに当てはめて解析するのは問題がある。（3）に関しては以前、非線形最小二乗法でも取り扱った。あるいは、今後、一般化線形モデルでも解決する術を紹介することになるだろう。今回は、（1）や（2）も絡んでくる状況を考えよう。

外れ値の存在

　例えば、下記のようなデータを解析することを考える。

------------------------------------------------------

library(plotn)

x1 <- c(10.1, 9.43, 2.34, 5.28, 17.32, 15.13, 12.68, 5.12, 

        28.99, 18.44, 19.11, 28.44, 18.43, 16.87, 23.41, 

        2.61, 13.98, 6.4, 8.7, 19.95, 10.08, 26.24, 22.27,

        0.44, 20.61, 10.97, 15.66, 17.43, 14.97, 30)

y <- c(11.51, 10.26, 4.85, 5.88, 17.62, 16.34, 13.03, 4.4,

       29.52, 20.36, 17.83, 30.28, 20.39, 18.41, 23.31, 4.38, 

       13.77, 5.98, 9.85, 22.79, 11.05, 27.29, 22.61, 1.63, 

       21.96, 15.15, 16.64, 19.84, 13.62, 120)

d <- data.frame(x1 = x1, y = y)

head(d)

##      x1     y

## 1 10.10 11.51

## 2  9.43 10.26

## 3  2.34  4.85

## 4  5.28  5.88

## 5 17.32 17.62

## 6 15.13 16.34

plotn(y ~ x1, data = d)#図1の描画

------------------------------------------------------

ほとんどのデータは分散も小さく、1直線上に乗っているように見えるが、1点のデータのみ、明らかに外れている。このような、予測から"大きく"外れているような値のことを外れ値outlierと呼んだりする。外れ値に具体的な定義が必ずしもあるわけではないので、カッコ付きで"大きく"と表現した。こんなとき、通常の線形回帰を行い、回帰直線を描くと以下の通りになる。

------------------------------------------------------

fit1 <- lm(y ~ x1, data = d)

summary(fit1)

## 

## Call:

## lm(formula = y ~ x1, data = d)

## 

## Residuals:

##     Min      1Q  Median      3Q     Max 

## -13.168  -5.266  -3.322   1.632  75.597 

## 

## Coefficients:

##             Estimate Std. Error t value Pr(>|t|)    

## (Intercept)  -6.5247     6.0789  -1.073    0.292    

## x1            1.6976     0.3572   4.752 5.46e-05 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 15.56 on 28 degrees of freedom

## Multiple R-squared:  0.4465, Adjusted R-squared:  0.4267 

## F-statistic: 22.59 on 1 and 28 DF,  p-value: 5.457e-05

plotn(y ~ x1, data = d) #図2の描画

xx <- seq(min(x1), max(x1), length = 200)

yy1 <- coef(fit1)[1] + coef(fit1)[2] * xx

overdraw(points(xx, yy1, type = "l", col = "green"))

------------------------------------------------------

30点あるデータのうち29のデータは1直線状なのに、ある1点の外れ値に引っ張られて実測値と回帰直線の間に乖離が生まれている。

外れ値の定義は恣意的になりがちであり、この外れ値をどう取り扱うかは諸説ある。一切の恣意性を排除し、上記の回帰で終わらせるのも手であろう。しかし、外れ値以外の多くのデータに関する予測が、外れ値の存在で大きく外れるのを避けたい場面もあるだろう。その時は、外れ値を除外、という手もあり得る。この時は、どこまでを外れ値とするかをよく考えなければならない。

今回紹介するのは、除外ではなく、外れ値も使いつつ予測を行う、という手法である。その手法の一つが、重み付き最小二乗法weighted ordinary least square (WLS) method である。

重み付き最小二乗法weighted ordinary least square (WLS) method

なぜ、1/30の外れ値に強く影響を受けるかといえば、通常の最小二乗法 (OLS) は、すべてのデータの「重要度」を同じと捉えていることが原因である。これを式で表したものが、線形回帰の仮定の2つ目、残差の分散が説明変数によらず一定である、という仮定である。先取りになるが、OLSの妥当性は、より一般の回帰分析に拡張された手法である最尤法からも導き出される。残差の分散が一定である場合、一つの説明変数だけを持つモデル（推定する係数は切片β0と傾きβ1）は最尤法によって以下の式を解けばいいことになる（ここでは結論だけ示す。詳しくは該当ページ）。

ここで、σが残差の分散を示すのであるが、分散が一定であるという仮定の下では、σは定数だから、結局、分子の最小化と同じ問題となる。

では、説明変数に依存して、分散が変わる場合は、上記のように式は簡単にならない。以下のように、分母に変数となった分散が残ったうえで、最小化の問題を解く必要がある。

この式は、分散一定の時の残差平方和を各データ点における分散で「調整した」平方和を最小化する問題であることを意味する。後述するが、この調整が「重み」と呼ばれる、重要度を示す値である。

　この式の心をもう少し、掘り下げる。もし、あるデータ点における分散が大きいならば、残差平方は小さくなるのだから、そのデータ点の残差平方和への影響は小さくなる。一方で、分散が小さければ、そのデータ点の残差平方和への影響は大きくなる。図で示すなら以下のようなイメージだ。

そこで、より一般に以下のような式を最小化することを考える。ここでwiで示す値が「重み」である。この方法が、重み付き最小二乗法である。まさに、名前通りである。

重みは自由に決定することができるが、特に妥当と考えられる調整方法は、上記の議論からwi = (1/σi)^2とする場合だろう。

　しかしながら、実際には各データ点の分散は未知である。ゆえに、真の分散の代替となる値が必要になる。これに関しては、明確な答えがあるわけではない。以下には、通常のOLSをおこない、そこから求まる予測値と観察値の差の2乗の逆数を重みとすることを考えることにする。つまり、図2の回帰直線（緑）と観察値の差の2乗の逆数を重みとする。おおむね、この重みは、真のデータの分散に依存した数値となっていると予測できるだろう。

R上では、lm関数の引数weightに上記の重みを入力することでWLSを実行できる。以下のとおりである。重みは自由に設定できると述べたが、予測値と観察値の差の2乗の逆数を重みと類似した効果を持つ重みとしては、予測値と観察値の差の絶対値の逆数もありえる。ここでは、その調整の結果も示している。

------------------------------------------------------

w <- abs(y - (coef(fit1)[1] + coef(fit1)[2] * x1)) #通常のOLSの予測値と観察値の差の絶対値

fit2 <- lm(y ~ x1, data = d, weights = 1/w)#絶対値の逆数を重みとする

summary(fit2)

## 

## Call:

## lm(formula = y ~ x1, data = d, weights = 1/w)

## 

## Weighted Residuals:

##     Min      1Q  Median      3Q     Max 

## -1.2079 -0.6160 -0.0884  0.3803 10.0098 

## 

## Coefficients:

##             Estimate Std. Error t value Pr(>|t|)    

## (Intercept)  -0.1172     1.4801  -0.079    0.937    

## x1            1.1028     0.1158   9.523  2.8e-10 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 2.006 on 28 degrees of freedom

## Multiple R-squared:  0.7641, Adjusted R-squared:  0.7557 

## F-statistic: 90.68 on 1 and 28 DF,  p-value: 2.801e-10

yy2 <- coef(fit2)[1] + coef(fit2)[2] * xx

plotn(y ~ x1, data = d)#図3の描画

overdraw(points(xx, yy1, type = "l", col = "green"),

         points(xx, yy2, type = "l", col = "blue"))

fit3 <- lm(y ~ x1, data = d, weights = 1/w^2)#差の2乗の逆数を重みとする

summary(fit3)

## 

## Call:

## lm(formula = y ~ x1, data = d, weights = 1/w^2)

## 

## Weighted Residuals:

##      Min       1Q   Median       3Q      Max 

## -0.64397 -0.21262 -0.05898  0.18448  1.15558 

## 

## Coefficients:

##             Estimate Std. Error t value Pr(>|t|)    

## (Intercept)  0.09229    0.53976   0.171    0.865    

## x1           1.08499    0.05122  21.183   <2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 0.4101 on 28 degrees of freedom

## Multiple R-squared:  0.9413, Adjusted R-squared:  0.9392 

## F-statistic: 448.7 on 1 and 28 DF,  p-value: < 2.2e-16

yy3 <- coef(fit3)[1] + coef(fit3)[2] * xx

overdraw(points(xx, yy3, type = "l", col = "orange"))

------------------------------------------------------

すると、WLSでは外れ値の影響が低減し、29/30のデータ側に沿うような回帰直線となった。このように、WLSを用いることで、必ずしもデータの分散が一定ではない場合、特に外れ値が存在するときに、うまく回帰を実行することが可能となる。

補足：重みの決め方は上記でいいのか？　ほかに方法はないのか？

さて、重みとして、分散の推定値を、通常のOLSから得られる予測値と観測値の差の2乗として計算を行ってきた。重みの決め方は本当にこれでいいのだろうか。確かに分散と関連している値ではあるはずだけど、外れ値に対する重みは過大評価されているように見える。真の重みとするべきなのは、例えば、外れ値を除外したデータで線形回帰し、その予測値と観測値の差にするべきではないか。しかし、これでは同じ問題の繰り返しである。つまり、外れ値とは何か、を明確に定義しないと恣意的なデータ除外になりかねない。けれども、気になるのが人間の性だ。

ところで、上記の重みの設定のもと、WLSをおこなうと、外れ値が「無視」されたかのような回帰直線になっている。実際には考慮されているが、重要度が低減しているので、外れ値以外におおむね合わせたような回帰直線になっているということだ。では、新たにこのWLSで得られた予測値をもとに、観察値との差を重みとすれば、お望みの真の分散に近い重みを得られるのではなかろうか。このような発想の下、繰り返し、WLSをおこなってパラメータの更新を行う方法を反復再重み付き最小二乗法Iteratively reweighted least squares (IRLS) methodと呼ぶ。

このように実行される回帰分析は外れ値の影響を受けにくいため、頑健な（ロバストrobust）回帰分析の一種として取り扱われる。ただし実際には、上記のようなシンプルな重みの設定ではなく、もっと効率的な重みが設定されているようだ（TukeyのbiweightやHuberの方法）が、発想は上記のもので問題ない。

　ほかの方法として、一般化最小二乗法Genenalized least square (GLS) method（とそれを実際に計算可能にしている方法Feasible generalized least squares（fGLS））という手法があるようだが、私の理解はまだ進んでいないので名称の紹介にとどめる。