重み付き回帰分析(WLS)について
重み付き回帰分析(WLS)ってなに？
重み付き回帰分析とは、加重最小二乗法やWeighted least squares
(WLS)などと呼ばれる
回帰分析を行う際の係数の推定方法の一つです。
通常の最小二乗法(Ordinary least
squares⇒OLS)と何が違うかというと、
通常、回帰分析に掛けるデータは等分散性を前提としているのですが、
どうしても観測データに不等分散性が生じる場合もあり、それを加重して補正しましょう、といった場合に使われます。
今回は説明を簡単にするため、単回帰分析で行います。
まず参考として単回帰分析(OLS)の基本モデルを列記します。
\[ Y_i = \beta_0 + \beta_1X_i + \epsilon_i\\ \begin{array}{l} Y_i&:&目的変数(従属変数)\\ X_i&:&説明変数(独立変数)\\ \beta_0,\beta_1&:&回帰係数\\ \epsilon_i&:&誤差項(平均0、分散\sigma^2)\\ \end{array} \]
例えば下記のようなデータ
\[ y = 3x + 2 \]
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 2.041 | 0.03183 | 64.12 | 0 |
| x | 3.088 | 0.03211 | 96.18 | 0 |
| Observations | Residual Std. Error | \(R^2\) | Adjusted \(R^2\) |
|---|---|---|---|
| 1000 | 1.006 | 0.9026 | 0.9025 |
これは一般的な単回帰分析ですが、上の式を元に乱数を発生させた結果です。
回帰分析の結果もそれなりに当てはまりが良く結果を得られています。
次に、同じxとyの関係について、例えばxの値が裾野に広がるほど、
yの値も暴れる(分散が大きくなる)、といったかたちを作ってみます。
yの値を乱数で出す際に\(\sigma^2 =
x^2\)で生成しました。
set.seed(123)
x = rnorm(n)
s = x^2
y <- 2 + 3 * x + rnorm(n, sd = s)
lm <- lm(y ~ x)
par(las =1,ps = 8)
plot(x,y)
abline(lm$coefficients[1],lm$coefficients[2],col = 2)
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 2.028 | 0.05332 | 38.02 | 2.786e-196 |
| x | 3.344 | 0.05379 | 62.16 | 0 |
| Observations | Residual Std. Error | \(R^2\) | Adjusted \(R^2\) |
|---|---|---|---|
| 1000 | 1.686 | 0.7947 | 0.7945 |
先ほどのように係数が3、とはいかず、やや不安定となっています。
例えば比率を用いたデータなどでは、裾野にいくほど数値が不安定となる、
といったデータを取り扱うことがあるかと思います。
この影響を取り除くために重み付き回帰分析(WLS)が使われます。
どうやって影響を取り除くのか
数式的な説明は横断的にサイトを巡っていただくのが良いとは思いますが
単回帰係数の簡単な計算部分だけ下記に記載します。
\[ \hat{\beta}_1 = \frac{\sum_{i=1}^n w_i (X_i – \bar{X}_w)(Y_i – \bar{Y}_w)}{\sum_{i=1}^n w_i (X_i – \bar{X}_w)^2}\\ \hat{\beta}_0 = \bar{Y}_w – \hat{\beta}_1 \bar X_w \\ \\ ここで、\\ \bar{X}_w = \frac{\sum_{i=1}^n w_i X_i}{\sum_{i=1}^n w_i}\\ \bar{Y}_w = \frac{\sum_{i=1}^n w_i Y_i}{\sum_{i=1}^n w_i}\\ \]
今回は重みを下記のように設定しました。
\[ 加重値：w_i = \frac{1}{(x_i^2)^2}\\ \]
weights <- 1 / s^2
lm <- lm(y ~ x, weights = weights)
par(las =1,ps = 8)
plot(x,y)
abline(lm$coefficients[1],lm$coefficients[2],col = 2)
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 2 | 4.389e-06 | 455684 | 0 |
| x | 3.002 | 0.001994 | 1505 | 0 |
| Observations | Residual Std. Error | \(R^2\) | Adjusted \(R^2\) |
|---|---|---|---|
| 1000 | 1.01 | 0.9996 | 0.9996 |
重み付けをすることで精度の高い結果を得ることができました。
かなり無茶な重み付けも補正してくれる
先ほどは裾野に広がるほど数値が暴れる、というパターンを作成しましたが
逆に「裾野に広がるほど数値が大人しくなる、逆に中心付近で暴れる」場合で処理してみます。
set.seed(123)
x = rnorm(n)
s = 1/x^2
y <- 2 + 3 * x + rnorm(n, sd = s)
lm <- lm(y ~ x)
par(las =1,ps = 8)
plot(x,y)
abline(lm$coefficients[1],lm$coefficients[2],col = 2)
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 263.8 | 148.6 | 1.775 | 0.0762 |
| x | -0.3199 | 149.9 | -0.002134 | 0.9983 |
| Observations | Residual Std. Error | \(R^2\) | Adjusted \(R^2\) |
|---|---|---|---|
| 1000 | 4699 | 4.564e-09 | -0.001002 |
yの値が中心で非常に高い値となり、見づらいので拡大します。
このとおり、まったく合わない結果となりました。
この場合でも加重値値を正しく設定すれば、正しく補正してくれます。
weights <- 1 / s^2
lm <- lm(y ~ x, weights = weights)
par(las =1,ps = 8)
plot(x,y,ylim = c(-300,300))
abline(lm$coefficients[1],lm$coefficients[2],col = 2)
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 2.003 | 0.01913 | 104.7 | 0 |
| x | 3.027 | 0.009147 | 330.9 | 0 |
| Observations | Residual Std. Error | \(R^2\) | Adjusted \(R^2\) |
|---|---|---|---|
| 1000 | 1.006 | 0.991 | 0.991 |
重み付けが与える影響について
今回は重み付けとしてこんな式を設定しました。
\[ 加重値：w_i = \frac{1}{(\frac{1}{x_i^2})^2}\\ \]
したがって、中心付近の数値が暴れている界隈では、影響度を小さくし
裾野に広がるほど影響度を高くする、といった補正を掛けて処理をしています。
x,y,加重(weights)の関係性は、下図のようなイメージです。
どんな重み付けをしていいかわからない場合
今回は意図的にxの値を用いてyの値を暴れさせていますが、
実際の観測データでは、どういった重み付けをしていいかわからない、という場合があります。
提案の一つとしては、通常の回帰分析(OLS)の残差を誤差として扱い、計算をする方法があります。
lm0 <- lm(y ~ x)
weights <- 1 / abs(lm0$residuals)
lm <- lm(y ~ x, weights = weights)
par(las =1,ps = 8)
plot(x,y,ylim = c(-300,300))
abline(lm$coefficients[1],lm$coefficients[2],col = 2)
| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 36.99 | 9.851 | 3.755 | 0.0001836 |
| x | 2.215 | 10.49 | 0.2111 | 0.8328 |
| Observations | Residual Std. Error | \(R^2\) | Adjusted \(R^2\) |
|---|---|---|---|
| 1000 | 20.35 | 4.467e-05 | -0.0009573 |
係数の線が先ほどの図よりちょっとだけ近づいたように思いますが、
あまり当てはまりが良いとは言えませんでした。
残差を元に繰り返し計算をするとどうなるのか。
先ほどの処理で、ほんのちょっとだけ、誤差が修正されました。
ということは、残差はやや小さくなり、その残差を加重値に設定することで
より精度が高まることが考えられます。
調べてみると、IWLS(Iteratively weighted least
squares)という考え方があるようです。
直訳すれば反復、繰り返し、ですね。
演算としては
・通常の回帰分析を行う
・残差を元に重み付けをした回帰分析を行う（WLS）
・繰り返す（IWLS）
といった手順です。
処理が多いのでアニメーション化しました。
左がxと加重値の関係性、右が数値が跳ねた場合も見えるよう対数化した図です。
赤い線は、理論値、今回で言えば\(w_i =
\frac{1}{(x_i^2)^2}\)の値に線を入れており
加重の修正が正しくされているかの指標としています。
回帰係数が3となるのはかなり初期の段階で達成されます。
今回は決定係数が高止まりしたら計算終了としましたが
アニメーションのとおり特定の値がどんどん基準から外れていくことが確認されます。
値が大きくなる＝今回は「当てはまりが良すぎる」値になるため
極端に当てはまりが良い値に対して強い加重が掛かる関係で
TSS(全平方和)とRSS(残差平方和)の差が離れていくことが確認できます。
式を補足掲載します。
\[ \begin{array}{l} \text{重み付き全平方和} & TSS &=& \sum_{i=1}^n w_i (Y_i – \bar{Y}_w)^2\\ \text{重み付き平均} & \bar{Y}_w &=& \frac{\sum w_i Y_i}{\sum w_i}\\ \text{重み付き残差平方和} & RSS &=& \sum_{i=1}^n w_i (Y_i – \hat{Y}_i)^2\\ 決定係数 & R^2 &=& 1 – \frac{\text{RSS}}{\text{TSS}}\\ \end{array}\\ \]
通常、残差が大きくなるということは、TSS(全平方和)に対してRSS(残差平方和)も大きくなり
決定係数が高くならないことが一般的ですが
加重を掛けた場合、wの補正がTSS、RSSに掛かり、
また残差の大きい値は無視（補正を小さく）、
残差の小さい値は補正を大きくする関係である段階から決定係数がどんどん高くなっていくようです。
これを解決するために実践の場では様々な手法が開発されているとは思いますが
また参考になる情報を見つけたらご紹介します。
