- rreg - Stataでロバスト回帰する
回帰分析したいけれど、外れ値をどうにかしてやりたい。良くある事だと思います（当社比）。
Stataではいくつかの方法がありますが、そのうちの1つであるrregについて、使用方法を学んだので、備忘録的にまとめておきたいと思います。
rregコマンドについては、下記の様に説明されています。なるほど、わからん。
rreg first performs an initial screening based on Cook’s distance> 1 to eliminate gross outliers before calculating starting values and then performs Huber iterations followed by biweight iterations, as suggested by Li(1985).
(Stata help fileから引用）
rregは、開始値を計算する前に、最初にCookの距離> 1に基づく初期スクリーニングを行い、総外れ値を除去し、Li(1985)によって提案されたように、Huber反復処理に続いてバイウェイト反復処理を行います。
(DeepL Proによる翻訳を少し編集)
Stataのヘルプファイルに準拠して、私が触った記録です。なお、使用したStataのバージョンは16です。
autoデータで回帰してみる（Example 1）
* データ読み込み
sysuse auto, clear
* 普通に回帰分析
regress mpg weight foreign
* 今回のテーマ：ロバスト回帰分析
rreg mpg weight foreignまず、最初の例では、Stataユーザが親の顔よりも見ているautoデータを用います。車の重さ・値段・燃費とかがまとめられているサンプルデータセットです。
回帰分析では、mpg（燃費、Mile Per Gallonの略）をweight（車体重量、ポンド）とforeign（米国車=0、外国車=1）で回帰しようとしています。
結果がこちらになります。
[result]
* 普通の回帰分析の結果
------------------------------------------------------------------------------
mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
weight | -.0065879 .0006371 -10.34 0.000 -.0078583 -.0053175
foreign | -1.650029 1.075994 -1.53 0.130 -3.7955 .4954422
_cons | 41.6797 2.165547 19.25 0.000 37.36172 45.99768
------------------------------------------------------------------------------[result]
* 今回のテーマ：ロバスト回帰分析
------------------------------------------------------------------------------
mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
weight | -.0063976 .0003718 -17.21 0.000 -.007139 -.0056562
foreign | -3.182639 .627964 -5.07 0.000 -4.434763 -1.930514
_cons | 40.64022 1.263841 32.16 0.000 38.1202 43.16025
------------------------------------------------------------------------------weightの回帰係数は、2つの分析で近い値になっていますが、foreignの回帰係数は大きく異なっています。何があったのでしょうか？
rregでは、解析上の重み（外れ値は重みが軽い）を計算し、その重みに従った解析を行っています。重みが軽い（極論、重み＝0もある）のはどのような車だったのかを見ることで何か分かるかも知れません。
* 推定された重みを新しい変数wに格納する。
rreg mpg weight foreign, genwt(w)
* 推定された重みを表示する。
summarize w, detail表示した結果がこちらです。
[results]
Robust Regression Weight
-------------------------------------------------------------
Percentiles Smallest
1% 0 0
5% .0442957 0
10% .4674935 0 Obs 74
25% .8894815 .0442957 Sum of Wgt. 74
50% .9690193 Mean .8509966
Largest Std. Dev. .2746451
75% .9949395 .9996715
90% .9989245 .9996953 Variance .0754299
95% .9996715 .9997343 Skewness -2.287952
99% .9998585 .9998585 Kurtosis 6.874605Smallestのところを見て下さい。3つのデータで重み=0になっています。つまり、解析から除外されたということです。どんな車が解析から除外されたのか見てみようと思います。
* 算出された重みで並び替え
sort w
* 算出された重みが小さいトップ5をリストアップ
list make mpg weight w in 1/5, sep(0)結果がこちら。
[results]
+-------------------------------------------+
| make mpg weight w |
|-------------------------------------------|
1. | VW Diesel 41 2,040 0 |
2. | Subaru 35 2,050 0 |
3. | Datsun 210 35 2,020 0 |
4. | Plym. Arrow 28 3,260 .04429567 |
5. | Cad. Seville 21 4,290 .08241943 |
+-------------------------------------------+車に詳しい人（？）なら、次の2つのことに気づくと思います。
一番上に上がっている「VW Diesel」は、このデータ中で唯一のディーゼル車です（外れ値）。また、Plymouth Arrowの車体重量は正確ではありません。Google先生によるとPlymouth Arrowの車体重量は2,200ポンドですが、ここでは3260ポンドと記載されています。
まぁ、なにかしら外れ値になる要因があるのでしょう、きっと（匙を投げた）。
重みをつけてregressしたらどうなる？
rregで重みが算出されたので、この重みを使ってregressで回帰分析して見ようと思います。
* pweightとして重みを付けたregress
regress mpg weight foreign [pweight=w]
* aweightとして重みをつけたregress
regress mpg weight foreign [aweight=w]pweight（サンプリングされる確率の逆数で重みとして指定）の場合と、aweight（分散に反比例する重みとして指定）の場合でやっています。
[results]
* pweightの場合
------------------------------------------------------------------------------
| Robust
mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
weight | -.0063976 .0003491 -18.33 0.000 -.0070942 -.005701
foreign | -3.182639 .612648 -5.19 0.000 -4.405159 -1.960119
_cons | 40.64022 1.213473 33.49 0.000 38.21878 43.06167
------------------------------------------------------------------------------
* aweightの場合
------------------------------------------------------------------------------
mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
weight | -.0063976 .0003317 -19.29 0.000 -.0070596 -.0057356
foreign | -3.182639 .5671798 -5.61 0.000 -4.314428 -2.050849
_cons | 40.64022 1.123826 36.16 0.000 38.39766 42.88278
------------------------------------------------------------------------------出力結果を並べて記載しました。両方の場合で回帰係数の点推定値はrregで行った場合と一致しています。しかし、標準誤差はrregよりも小さくなっていますので、注意が必要です。
説明変数を入れない場合はどうなる？（Example 2）
次は、rregに説明変数を入れずに実行します。
* 説明変数なしのモデル
rreg mpg
[results]
------------------------------------------------------------------------------
mpg | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_cons | 20.68825 .641813 32.23 0.000 19.40912 21.96738
------------------------------------------------------------------------------この場合rregは、平均値のロバスト推定を行います。
mpgの平均値のロバスト推定は20.69（95％信頼区間 19.41, 21.97）でした。この結果を非ロバスト推定と比較します。
* 非ロバスト推定で平均値を出す。
ci means mpg
[results]
Variable | Obs Mean Std. Err. [95% Conf. Interval]
-------------+---------------------------------------------------------------
mpg | 74 21.2973 .6725511 19.9569 22.63769非ロバスト推定では、21.30（95%信頼区間 19.96, 22.64）でした。
ロバスト推定の平均値がやや低い値になったのは、燃費が良すぎる外れ値で軽い重みが付与されたからだと考えられます。
沼へようこそ
Stata ユーザーコミュニティには、ロバスト統計の分野があり、活発な議論が行われています。沼にようこそ。
利益相反（COI）の開示
金銭・経済的なCOIもありません。ただし、金銭を頂くことを拒否している訳ではありません。何か贈りたい方は是非お願いします（ダイマ）。