ロバスト回帰を使用した外れ値の影響の低減

ロバスト線形回帰を使用して、線形回帰モデルの外れ値の影響を減らすことができます。このトピックでは、ロバスト回帰について定義し、それを使用して線形モデルを当てはめる方法を示し、結果を標準近似と比較します。fitlm を名前と値のペアの引数 'RobustOpts' と共に使用して、ロバスト回帰モデルを当てはめることができます。または、robustfit を使用して、単にロバスト回帰係数パラメーターを計算できます。

ロバスト回帰を使用する理由

ロバスト線形回帰は、標準の線形回帰よりも外れ値の影響を受けにくくなります。標準の線形回帰では、通常の最小二乗近似を使用して、応答データを 1 つ以上の係数をもつ予測子データに関連付けるモデルパラメーターを計算します (詳細については、多変量回帰モデルの推定を参照してください)。残差を二乗するとこれらの極値データ点の影響を増幅することになるため、結果として、外れ値が近似に大きく影響します。線形回帰モデルとはで説明した標準線形回帰を使用するモデルは、観測された応答において誤差が正規分布するというような、ある仮定に基づきます。誤差の分布が非対称または外れ値になる傾向がある場合、モデルの仮定は無効になり、パラメーターの推定、信頼区間、その他の計算された統計量が信頼できなくなります。

ロバスト回帰では、反復的に再重み付けした最小二乗と呼ばれる方法を使用して、各データ点に重みを割り当てます。この方法では、狭い部分のデータでの大きな変化に対して感度が低くなります。結果として、ロバスト線形回帰は、標準の線形回帰よりも外れ値の影響を受けにくくなります。

反復的に再重み付けした最小二乗

重み付き最小二乗では、近似プロセスに追加のスケール係数としての重みを含めて、近似を向上させます。重みは、各応答値が最終的なパラメーター推定にどの程度影響するかを決定します。品質が低いデータ点 (たとえば、外れ値) は、近似への影響が少なくなければなりません。重み w_i を計算するには、テューキーの二重平方関数などの定義済みの重み関数を使用できます (その他のオプションについては、fitlm の名前と値のペアの引数 'RobustOpts' を参照)。

"反復的に再重み付けした最小二乗" のアルゴリズムは自動的かつ反復的に重みを計算します。アルゴリズムは初期化時に各データ点に等しい重みを割り当て、通常の最小二乗を使用してモデル係数を推定します。アルゴリズムは各反復で、前の反復においてモデル予測から遠かった点にはより小さな重みが与えられるように、重み w_i を計算します。次に、アルゴリズムは重み付き最小二乗を使用してモデル係数 b を計算します。係数推定の値が指定した許容誤差内に収束すると、反復が停止します。このアルゴリズムは、最小二乗法を使用してデータの大部分に当てはまる曲線を求めると同時に、外れ値の影響を最小限に抑えようとします。

詳細は、反復的に再重み付けした最小二乗の手順を参照してください。

標準およびロバスト最小二乗近似の結果の比較

ライブスクリプトを開く

この例では、関数 fitlm でロバスト回帰を使用する方法を示し、ロバスト近似と標準の最小二乗近似の結果を比較します。

moore データを読み込みます。先頭の 5 列に予測子データがあり、6 列目に応答データが含まれています。

load moore
X = moore(:,1:5);
y = moore(:,6);

最小二乗線形モデルをデータに当てはめます。

mdl = fitlm(X,y)

mdl = 
Linear regression model:
    y ~ 1 + x1 + x2 + x3 + x4 + x5

Estimated Coefficients:
                    Estimate          SE          tStat       pValue 
                   ___________    __________    _________    ________

    (Intercept)        -2.1561       0.91349      -2.3603      0.0333
    x1             -9.0116e-06    0.00051835    -0.017385     0.98637
    x2               0.0013159     0.0012635       1.0415     0.31531
    x3               0.0001278    7.6902e-05       1.6618     0.11876
    x4               0.0078989         0.014      0.56421     0.58154
    x5              0.00014165    7.3749e-05       1.9208    0.075365


Number of observations: 20, Error degrees of freedom: 14
Root Mean Squared Error: 0.262
R-squared: 0.811,  Adjusted R-Squared: 0.743
F-statistic vs. constant model: 12, p-value = 0.000118

名前と値のペアの引数 'RobustOps' を使用して、ロバスト線形モデルをデータに当てはめます。

mdlr = fitlm(X,y,'RobustOpts','on')

mdlr = 
Linear regression model (robust fit):
    y ~ 1 + x1 + x2 + x3 + x4 + x5

Estimated Coefficients:
                    Estimate         SE         tStat       pValue 
                   __________    __________    ________    ________

    (Intercept)       -1.7516       0.86953     -2.0144    0.063595
    x1             1.7006e-05    0.00049341    0.034467     0.97299
    x2             0.00088843     0.0012027      0.7387     0.47229
    x3             0.00015729    7.3202e-05      2.1487    0.049639
    x4              0.0060468      0.013326     0.45375     0.65696
    x5             6.8807e-05    7.0201e-05     0.98015     0.34365


Number of observations: 20, Error degrees of freedom: 14
Root Mean Squared Error: 0.249
R-squared: 0.775,  Adjusted R-Squared: 0.694
F-statistic vs. constant model: 9.64, p-value = 0.000376

2 つのモデルの残差を視覚的に調べます。

tiledlayout(1,2)
nexttile
plotResiduals(mdl,'probability')
title('Linear Fit')
nexttile
plotResiduals(mdlr,'probability')
title('Robust Fit')

ロバスト近似の残差 (プロットの右半分) は直線に近く、顕著な外れ値が 1 つのみあります。

外れ値のインデックスを見つけます。

outlier = find(isoutlier(mdlr.Residuals.Raw))

outlier = 1

ロバスト近似の観測値の重みをプロットします。

figure
b = bar(mdlr.Robust.Weights);
b.FaceColor = 'flat';
b.CData(outlier,:) = [.5 0 .5];
xticks(1:length(mdlr.Residuals.Raw))
xlabel('Observations')
ylabel('Weights')
title('Robust Fit Weights')

Figure contains an axes object. The axes object with title Robust Fit Weights, xlabel Observations, ylabel Weights contains an object of type bar.

ロバスト近似の外れ値の重み (紫色のバー) は、他の観測値の重みよりもはるかに小さい値です。

反復的に再重み付けした最小二乗の手順

反復的に重み付けした最小二乗アルゴリズムは、次の手順に従います。

重みの初期推定から始めて、重み付き最小二乗でモデルを当てはめます。
調整済み残差を計算します。調整済み残差は次によって与えられます。
$r_{adj} = \frac{r_{i}}{\sqrt{1 - h_{i}}}$
ここで、r_i は通常の最小二乗残差、h_i は最小二乗近似のてこ比値です。てこ比による残差の調整では、最小二乗近似に大きな影響を与える、高いてこ比のデータ点の重みを小さくします (ハット行列とてこ比を参照)。
残差を標準化します。標準化された調整済み残差は次によって与えられます。
$u = \frac{r_{adj}}{K s} = \frac{r_{i}}{K s \sqrt{1 - h_{i}}}$
ここで、K は調整定数、s は s = MAD/0.6745 によって与えられる誤差項の標準偏差の推定値です。
MAD は、残差の中央値に対する残差の中央絶対偏差です。定数 0.6745 は、推定を正規分布に対して不偏にします。予測子データ行列 X に p 個の列がある場合、中央値を計算するときに、最小のものから p 個の絶対偏差が除外されます。
u の関数としてロバストな重み w_i を計算します。たとえば、二重平方重みは次によって与えられます。
$w_{i} = {\begin{array}{l} {(1 - u_{i}^{2})}^{2} & , | u_{i} | < 1 \\ 0 & , | u_{i} | \geq 1 \end{array}$
ロバスト回帰係数 b を推定します。重みによって、パラメーター推定 b の式は次のように変更されます。
$b = \hat{β} = {(X^{T} W T)}^{- 1} X^{T} W y$
ここで、W は対角重み付け行列、X は予測子データ行列、y は応答ベクトルです。
次の重み付き最小二乗誤差を推定します。
$e = \sum_{1}^{n} w_{i} {(y_{i} - {\hat{y}}_{i})}^{2} = \sum_{1}^{n} w_{i} r_{i}^{2}$
ここで、w_i は重み、y_i は観測応答、ŷ_i は近似応答、r_i は残差です。
近似が収束するか、最大反復回数に達すると、反復が停止します。そうでない場合は、2 番目の手順に戻って、最小二乗近似の次の反復を実行します。

参考

fitlm | robustfit | LinearModel | plotResiduals