robustbaseパッケージではMM推定、S推定、Koller & Stahel (2017)、Koller (2012)による補正設定を用いたMM推定が可能である。

MM推定

MM推定はYohai(1987)によって提案された方法である。MM推定はM推定で得られた残差標準偏差を最小化するS推定を用いて回帰パラメータを推定した上で、M推定を行うものである。MM推定の目的は、ブレークダウン値が高く、より効率的な推定値を得ることである。

Yohai, V.J. (1987) High breakdown-point and high efficiency estimates for regression. The Annals of Statistics 15, 642–65.

library(robustbase)
mod.robustbase.MM <- lmrob(wage ~ gender + age + education, 
                        method="MM", data = CPS1985)
summary(mod.robustbase.MM)

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -4.20343    1.38053  -3.045  0.00244 ** 
genderfemale -2.22148    0.34372  -6.463 2.33e-10 ***
age           0.11623    0.01634   7.114 3.68e-12 ***
education     0.72453    0.09022   8.030 6.32e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Robust residual standard error: 3.439 
Multiple R-squared:  0.294, Adjusted R-squared:   0.29

S推定

S推定は、RousseeuwとYohai(1984)によって提案された方法。M推定の残差標準偏差(residual scale)を元にして計算する。M推定の弱点は、重み付け値として中央値のみを使用するため、データ分布が考慮されておらず、データ全体の関数ではないことである。S推定は、中央値の弱点を克服するために残差標準偏差を用いる。

1. 1. Rousseeuw and V. J. Yohai, Robust Regression by Mean of SEstimators, Robust and Nonlinear Time Series Analysis, New York, 1984, 256-274, doi: 10.1007/978-1-4615-7821-5-15.

library(robustbase)
mod.robustbase.S <- lmrob(wage ~ gender + age + education, 
                        method="S", data = CPS1985)
summary(mod.robustbase.S)

結果

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -4.20547    1.43947  -2.922  0.00363 ** 
genderfemale -1.72424    0.43121  -3.999 7.28e-05 ***
age           0.12670    0.01895   6.688 5.77e-11 ***
education     0.62001    0.08981   6.903 1.46e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Robust residual standard error: 3.439 
Multiple R-squared:  0.5189,    Adjusted R-squared:  0.5162

set.seed(33)
x1 <- sort(rnorm(30)); x2 <- sort(rnorm(30)); x3 <- sort(rnorm(30))
X. <- cbind(x1, x2, x3)
y <-  10 + X. %*% (10*(2:4)) + rnorm(30)/10
y[1] <- 500   # a moderate outlier
X.[2,1] <- 20 # an X outlier
X1  <- cbind(1, X.)

Koller & Stahel (2017)、Koller (2012)による補正設定を用いたMM設定

library(robustbase)
mod.robustbase.KS2014 <- lmrob(wage ~ gender + age + education, 
                        method="MM", data = CPS1985, setting = "KS2014")
summary(mod.robustbase.KS2014)

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -4.18167    1.04744  -3.992 7.47e-05 ***
genderfemale -2.25123    0.32137  -7.005 7.52e-12 ***
age           0.11518    0.01389   8.291 9.33e-16 ***
education     0.72988    0.06323  11.543  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Robust residual standard error: 3.462 
Multiple R-squared:   0.29, Adjusted R-squared:  0.2859 
Convergence in 17 IRWLS iterations

コードにsetting = "KS2014"を入れる。
KS2014はKoller & Stahel (2017), Koller (2012)によって提案された設定である。

具体的には、KS2014設定では、サンプル再抽出の数（nResample）が1000に増やされ、最良の候補数（best.r.s）を20に設定し、ローカル改善ステップ数（k.fast.s）を2に設定する。これらの変更は、特にファクター変数を含む設計でより安定した推定を行うために設計されている。
また、この設定では、MM推定器（M-estimator with high breakdown point）のために使用される、標準化関数（\psi function）とチューニング定数（tuning.psi、tuning.chi）も特定の値に設定されています。これにより、標準誤差のロバスト性を高めつつ、データの異常値に対しても強い耐性を持たせることができる。この設定を利用することで、ロバスト回帰モデルのパフォーマンスを向上させ、異常値やモデルの仮定違反が結果に与える影響を小さく抑えることが可能となる。データ分析において一般的な仮定違反に対処する際に非常に有用とされている。

KS2014の元になったKS2011の論文はこちら。 - Koller, M., & Stahel, W. A. (2017). Nonsingular subsampling for regression S estimators with categorical predictors. Computational Statistics, 32(2), 631–646. https://doi.org/10.1007/s00180-016-0679-x

KS2014はarXivにあるままである。
- Koller, M. (2012). Nonsingular subsampling for S-estimators with categorical predictors. https://doi.org/10.48550/ARXIV.1208.5595

モデル比較

残差標準誤差の比較

各モデルの比較の一つの方法は各モデルの残差標準誤差を計算することである。残差標準誤差（residual standard error: RSE）は、回帰モデルの残差の標準偏差を測定する方法であり、RSE の値が低いほど、モデルがデータに適合的である。

https://www.rdocumentation.org/packages/robustbase/versions/0.99-2/topics/sigma

## 通常のOLS
sigma.lm <- sigma(mod) 
sigma.lm

## MM推定
sigma.MM <- sigma(mod.robustbase.MM) 
sigma.MM 

## S推定
sigma.S <- sigma(mod.robustbase.S)
sigma.S

## KS2014設定
sigma.KS2014 <- sigma(mod.robustbase.KS2014) 
sigma.KS2014

結果：

## 通常のOLS
4.454083

## MM推定
3.439039

## S推定
3.439039

## KS2014設定
3.461537

R二乗

各モデルの出力にある。

Multiple R-squared:  0.253,  Adjusted R-squared:  0.2488 
Multiple R-squared:  0.294, Adjusted R-squared:   0.29 
Multiple R-squared:  0.294, Adjusted R-squared:   0.29 
Multiple R-squared:   0.29, Adjusted R-squared:  0.2859