平成も終わるので、ロバスト統計の話をしよう

この記事は[R Advent Calendar 2018](https://qiita.com/advent-calendar/2018/rlang) 15日目の記事です。
今回は、平成も終わるのでロバスト統計の話をしますね。(異論は受け付けない)

## 外れ値
何かしらの実験やデータを取った際に、分布を確認すると思います。
そんな時に、いくつかおかしなところにplotされてしまう何かがあるケース、よく出くわしますよね。
他のデータよりも明らかにおかしな値を取っているものを<b>外れ値(Outlier)</b>と言います。
具体的には下記のようなイメージです。
<img width="500" alt="スクリーンショット 2018-12-06 14.24.10.png" src="https://qiita-image-store.s3.amazonaws.com/0/38233/1754a28f-bc46-d6dd-3df3-5773594652c3.png">

これを元に要約統計量などを算出しようとすると、外れ値の持つ値に引きずられて平均値などの代表値は信頼性を担保できなくなる可能性があります。

## ロバスト統計とは
先ほど説明した外れ値の混入は未然に防ぐことができれば良いですね。しかし、データ取得の前に防げば良いのですが、なかなか思うようにうまくはいきません。サービスをリリースしたりすると、想定していなかった何かが起こることは日常茶飯事です。
そういう時には、外れ値などの混入を起こりうるものとして、対処する必要があります。
このように、外れ値が混入した場合でもメインとなる代表性を担保できるような統計モデルを組み立てたり、推定や検定を行うことを<b>ロバスト統計</b>と言います。


何か外れ値が混入した時に、それを検知するモデルを作るようなアプローチをとると<b>異常検知(Outlier Detection)</b>となります。今回はこちらの話ではありません。

ロバストは日本語で訳すと<b>頑健</b>と言われます。
ロバストネスに関しては、ブラックスワンで有名なナシーム・ニコラス・タレブが著書の「On Robustness And Fragility」(邦訳：強さと脆さ)で取り上げたことでも有名です。
<img src="https://images-na.ssl-images-amazon.com/images/I/41EoMTx-h9L._SX348_BO1,204,203,200_.jpg" width=180>

## 中央値
では、どう対処したら良いのでしょうか。結論から言うと、<b>中央値</b>を用いるのがベターです。
代表値の3つ。最頻値(mode)、平均値(mean)、中央値(median)の違いはWikipediaの下記の図がわかりやすいです。
<img src="https://upload.wikimedia.org/wikipedia/commons/thumb/3/33/Visualisation_mode_median_mean.svg/400px-Visualisation_mode_median_mean.svg.png" width=280>

このように平均では、分布が右に裾野を引いていると、その影響を受けてしまいやすく、中央値の場合は累積50%の値になるため、1つや2つの異常値が混入したとしても、それに引きづられることはありません。

## M推定
ロバスト推定で最もポピュラーな手法です。
例えば、データX={x1,...xn}がある母集団からの標本だとします。
この母集団の分布を確率密度関数f(x:θ)とすると、最尤推定量は<img width="300" alt="スクリーンショット 2018-12-06 15.03.51.png" src="https://qiita-image-store.s3.amazonaws.com/0/38233/b7e561f1-bda7-1dd5-8ba8-ae92174cc387.png">
これを次のように書き換えます。
<img width="300" alt="スクリーンショット 2018-12-06 15.06.24.png" src="https://qiita-image-store.s3.amazonaws.com/0/38233/ae598261-53f2-6a09-da83-157152aacb71.png">
<img width="300" alt="スクリーンショット 2018-12-06 15.08.33.png" src="https://qiita-image-store.s3.amazonaws.com/0/38233/3b0f67d5-6165-8b87-a7e3-a801321290cc.png">

このようにすると、ロス関数を最小化する推定を行うことになります。これがM推定です。
一言で言えば、誤差があるデータに対してその誤差の影響を最小にするものです。

## Rで使うには
{robustbase}パッケージにはかなり豊富にロバスト統計を行える関数が用意されています。
執筆時点での[CRANの最新情報](https://cran.r-project.org/web/packages/robustbase/index.html)だとアップデートも今年の9月に行われており、メンテナンスもされていそうです。

個人的にはこのパッケージ内に用意されているデータセットが豊富なことがテンションが上がるポイントでした。

## データセット一覧
| Data | Description |
| --- | --- |
| Animals2 | Brain and Body Weights for 65 Species of Land Animals |
| CrohnD | Crohn's Disease Adverse Events Data |
| NOxEmissions | NOx Air Pollution Data |
| SiegelsEx | Siegel's Exact Fit Example Data |
| aircraft | Aircraft Data |
| airmay | Air Quality Data |
| alcohol | Alcohol Solubility in Water Data |
| ambientNOxCH | Daily Means of NOx (mono-nitrogen oxides) in air |
| biomassTill | Biomass Tillage Data |
| bushfire | Campbell Bushfire Data |
| carrots | Insect Damages on Carrots |
| cloud | Cloud point of a Liquid |
| coleman | Coleman Data Set |
| condroz | Condroz Data |
| cushny | Cushny and Peebles Prolongation of Sleep Data |
| delivery | Delivery Time Data |
| education | Education Expenditure Data |
| epilepsy | Epilepsy Attacks Data Set |
| exAM | Example Data of Antille and May - for Simple Regression |
| foodstamp | Food Stamp Program Participation |
| hbk | Hawkins Bradu Kass's Artificial Data |
| heart | Heart Catherization Data |
| kootenay | Waterflow Measurements of Kootenay River in Libby and Newgate |
| lactic | Lactic Acid Concentration Measurement Data |
| los | Length of Stay Data |
| milk | Daudin's Milk Composition Data |
| pension | Pension Funds Data |
| phosphor | Phosphorus Content Data |
| pilot | Pilot-Plant Data |
| possum.mat (possumDiv) | Possum Diversity Data |
| possumDiv | Possum Diversity Data |
| pulpfiber | Pulp Fiber and Paper Data |
| radarImage | Satellite Radar Image Data from near Munich |
| salinity | Salinity Data |
| starsCYG | Hertzsprung-Russell Diagram Data of Star Cluster CYG OB1 |
| steamUse | Steam Usage Data (Excerpt) |
| telef | Number of International Calls from Belgium |
| toxicity | Toxicity of Carboxylic Acids Data |
| vaso | Vaso Constriction Skin Data Set |
| wagnerGrowth | Wagner's Hannover Employment Growth Data |
| wood | Modified Data on Wood Specific Gravity |

## データを眺める
こうなると、用意されているデータがどんな按配か確認したくなるのがRユーザーの性ですよね。
普通にpairs関数で散布図行列を描画しても良いのですが、せっかくなので`GGally::ggpairs`で綺麗なplotを描いてみましょう

## ggpairs(aircraft)
<img width="800" alt="スクリーンショット 2018-12-06 15.16.17.png" src="https://qiita-image-store.s3.amazonaws.com/0/38233/8c2bf977-cd66-2493-91be-a1937198152f.png">

## ggpairs(carrots)
<img width="800" alt="スクリーンショット 2018-12-07 15.22.43.png" src="https://qiita-image-store.s3.amazonaws.com/0/38233/0227ee39-770c-4cc0-eecb-5bf30abf121c.png">

## ggpairs(milk)
<img width="800" alt="スクリーンショット 2018-12-06 15.17.58.png" src="https://qiita-image-store.s3.amazonaws.com/0/38233/ac76d463-ada9-48bc-98d4-cd8abfaef2f2.png">

## ggpairs(wood)
<img width="800" alt="スクリーンショット 2018-12-06 15.19.35.png" src="https://qiita-image-store.s3.amazonaws.com/0/38233/723263cf-c91f-18d7-a25c-737377cbb58b.png">
あー。いいですねぇ。ワクワクしますねぇ。

## ロバスト回帰モデル
それでは、実際にパッケージに用意されているDatasetを使ってロバスト回帰を行ってみましょう。
一般的なロバスト回帰である`robustbase::lmrob`もありますが、なんとなく気分が乗らなかったので、一般化線形モデルベースの`robustbase::glmrob`を使ってみます。
`robustbase::glmrob`は柔軟にモデリングができ、分布の指定は`family = binomial, = poisson, = Gamma
and = gaussian`と4種類の選択が可能です。今回は`binomial`を使います。

対象データはにんじんの損害データが記録された`robustbase::data(carrots)`を使います。

```R
## robustbase
library(robustbase)

data(carrots)
Cfit1 <- glm(cbind(success, total-success) ~ logdose + block,
             data = carrots, family = binomial)
summary(Cfit1)
```
まずは普通にglmでやってみます。

```R
Call:
glm(formula = cbind(success, total - success) ~ logdose + block, 
    family = binomial, data = carrots)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9200  -1.0215  -0.3239   1.0602   3.4324  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   2.0226     0.6501   3.111  0.00186 ** 
logdose      -1.8174     0.3439  -5.285 1.26e-07 ***
blockB2       0.3009     0.1991   1.511  0.13073    
blockB3      -0.5424     0.2318  -2.340  0.01929 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 83.344  on 23  degrees of freedom
Residual deviance: 39.976  on 20  degrees of freedom
AIC: 128.61

Number of Fisher Scoring iterations: 4
```
```R
Rfit1 <- glmrob(cbind(success, total-success) ~ logdose + block,
                family = binomial, data = carrots, method= "Mqle", #MrleはHuber Loss
                control= glmrobMqle.control(tcc=1.2))
summary(Rfit1)
```

```R
Call:  glmrob(formula = cbind(success, total - success) ~ logdose +      block, family = binomial, data = carrots, method = "Mqle",      control = glmrobMqle.control(tcc = 1.2)) 


Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   2.3883     0.6923   3.450 0.000561 ***
logdose      -2.0491     0.3685  -5.561 2.68e-08 ***
blockB2       0.2351     0.2122   1.108 0.267828    
blockB3      -0.4496     0.2409  -1.866 0.061989 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Robustness weights w.r * w.x: 
 15 weights are ~= 1. The remaining 9 ones are
     2      5      6      7     13     14     21     22     23 
0.7756 0.7026 0.6751 0.9295 0.8536 0.2626 0.8337 0.9051 0.9009 

Number of observations: 24 
Fitted by method ‘Mqle’  (in 9 iterations)

(Dispersion parameter for binomial family taken to be 1)

No deviance values available 
Algorithmic parameters: 
   acc    tcc 
0.0001 1.2000 
maxit 
   50 
test.acc 
  "coef"
```

```R
> coef(Cfit1)
(Intercept)     logdose     blockB2     blockB3 
  2.0226453  -1.8174044   0.3008816  -0.5423898 

> coef(Rfit1)
(Intercept)     logdose     blockB2     blockB3 
  2.3882515  -2.0491078   0.2351038  -0.4496314 

```
偏回帰係数の推定値の結果が異なっていますね。

`(Intercept)` と `logdose`に関しては有意水準もかなり良さそうです。
一方で `blockB3`は精度が落ちてしまいました。

この辺りはきちんとしたデータ量を増やすことで解決できるかもしれませんね。
現場からは以上です。
Enjoy！

明日は Atsushi776さんです。楽しみですね。