Help us understand the problem. What is going on with this article?

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

More than 1 year has passed since last update.

この記事はR Advent Calendar 2018 15日目の記事です。
今回は、平成も終わるのでロバスト統計の話をしますね。(異論は受け付けない)

外れ値

何かしらの実験やデータを取った際に、分布を確認すると思います。
そんな時に、いくつかおかしなところにplotされてしまう何かがあるケース、よく出くわしますよね。
他のデータよりも明らかにおかしな値を取っているものを外れ値(Outlier)と言います。
具体的には下記のようなイメージです。
スクリーンショット 2018-12-06 14.24.10.png

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

ロバスト統計とは

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

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

ロバストは日本語で訳すと頑健と言われます。
ロバストネスに関しては、ブラックスワンで有名なナシーム・ニコラス・タレブが著書の「On Robustness And Fragility」(邦訳:強さと脆さ)で取り上げたことでも有名です。

中央値

では、どう対処したら良いのでしょうか。結論から言うと、中央値を用いるのがベターです。
代表値の3つ。最頻値(mode)、平均値(mean)、中央値(median)の違いはWikipediaの下記の図がわかりやすいです。

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

M推定

ロバスト推定で最もポピュラーな手法です。
例えば、データX={x1,...xn}がある母集団からの標本だとします。
この母集団の分布を確率密度関数f(x:θ)とすると、最尤推定量はスクリーンショット 2018-12-06 15.03.51.png
これを次のように書き換えます。
スクリーンショット 2018-12-06 15.06.24.png
スクリーンショット 2018-12-06 15.08.33.png

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

Rで使うには

{robustbase}パッケージにはかなり豊富にロバスト統計を行える関数が用意されています。
執筆時点でのCRANの最新情報だとアップデートも今年の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)

スクリーンショット 2018-12-06 15.16.17.png

ggpairs(carrots)

スクリーンショット 2018-12-07 15.22.43.png

ggpairs(milk)

スクリーンショット 2018-12-06 15.17.58.png

ggpairs(wood)

スクリーンショット 2018-12-06 15.19.35.png
あー。いいですねぇ。ワクワクしますねぇ。

ロバスト回帰モデル

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

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

## robustbase
library(robustbase)

data(carrots)
Cfit1 <- glm(cbind(success, total-success) ~ logdose + block,
             data = carrots, family = binomial)
summary(Cfit1)

まずは普通にglmでやってみます。

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
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)
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"
> 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さんです。楽しみですね。

classi
学校の先生・生徒・保護者向けのB2B2Cの学習支援Webサービス「Classi(クラッシー)」 を開発・運営している会社です。
https://classi.jp/
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away