R上でロバストな標準化をしたい

データ分析をする際には各データを標準化するのが一般的だが、大抵の場合は特に何も考えずに平均0、分散1になるよう処理するのが基本だと思う。
これはRならscale()で実行できる。
ただし、この方法では外れ値によって大きく影響を受けることがある¹。

ここ最近になってその辺りを解決できる標準化が必要になったので調べていたら、ロバストzスコアというそのものずばりな方法があった。

biolab.sakura.ne.jp

統計ガチ勢の方々には当たり前の話だと思うけど、一般的な標準化はデータXにおける平均値 $\mu$ と標準偏差 $\sigma$ を用いて計算している。

$\displaystyle Z = \frac{X - \mu}{\sigma}$

この計算式における平均値を中央値へ、標準偏差を四分位範囲（IQR）²を標準正規分布に対応させた正規四分位範囲（NIQR）³へそれぞれ置き換えることによって、ロバストzスコアが算出できる。

式として書くと以下の通りである。

$\displaystyle Z = \frac{X - \overline{X}}{NIQR}$

という訳で、とりあえずRで実装してみる。
なお、いつも通り結果の貼り付け方がよく分かっていないので、基本的には手元の環境で試してください。

Rのbase関数にはIQR()が入っているが、NIQRはIQRを定数で割ればいい⁴だけなので関数としては存在しない。
今回は使い勝手を重視して自作関数にしておく。

#IQR()と同様の使い方ができるように変数を入れておく
NIQR <- function(x, na.rm = FALSE, type = 7){
  IQR(x, na.rm = na.rm, type = type) / (qnorm(0.75)-qnorm(0.25))
}

中央値とNIQRを使った標準化をやってみる。
対象データはRユーザーお馴染みのiris。分かりやすいように1行目だけで実施してみる。

dat <- iris[[1]]
(dat - median(dat)) / NIQP(dat)

もっと綺麗にコードを書けないか試していたら、scale()で簡単に実現できた。
scale()のヘルプ内にあるExampleには載っていない⁵が、centerとscaleにはロジカル（TRUE、FALSE）じゃなく整数も渡せる⁶ので、中央値とNIQRを渡せば処理してくれる。

scale(dat, median(dat), NIQP(dat))

こうすると出力される結果にcenterとscaleの情報も保存されるので、後で確認したい場合などには便利かもしれない。

上記では偏差としてNIQRを使用したが、他にも中央絶対偏差（MAD）⁷を用いて標準化する手法も存在する。
また、「MADは非対称性のあるデータにおいては性能が低下する」として、SnやQnというものも提案されている⁸。

調べた感じだとMADはそれなりに見かける印象だが、SnやQnは論文などで少し出てくる程度だった。
最近のコンピューターの性能的に「算出方法が複雑＝計算量が多いから使いにくい」ってことではないと思うのだが、詳細は不明⁹。
計算式はかなりややこしいので、注釈につけた論文を直接見てください。

これも細かい計算云々は抜きにして、試しに計算してみる。
MADの計算はRに関数として元々入っているのでそれを活用。
また、SnとQnについてはrobustbaseパッケージを入れることで簡単に扱えるようになる。

cran.r-project.org

mad(dat)

library(robustbase)
Sn(dat)
Qn(dat)

robustbase::s_*(x, mu.too=T)を使うと、中央値と偏差の両方が同時に取れる。
共分散などの計算式へ値を渡す時に活躍するらしい。

s_mad(dat, mu.too=T)
s_Sn(dat, mu.too=T)
s_Qn(dat, mu.too=T)

MAD、Sn、Qnを用いた標準化は、上でやったのと同様の処理によって実現できる¹⁰。

scale(dat, median(dat), mad(dat))
scale(dat, median(dat), Sn(dat))
scale(dat, median(dat), Qn(dat))

今回調べた結果として標準化や偏差の計算手法は色々と確認できたのだが、それぞれの使い分けとかについてはほとんど情報が得られなかった。
とりあえずIQRかMADを用いておけば問題はなさそうだけど、またどこかのタイミングで改めて確認したい。

Enjoy!!

例えばc(2, 3, 5, 6, 9)のデータに1000とかを混ぜて標準化すると悲惨なことになる（参考元：https://en.wikipedia.org/wiki/Robust_statistics#Examples）。↩
昇順で並び替えたデータの25%～75%の範囲のこと。Rではsummary(x)やquantile(x)で確認することができる。↩
IQRを『対応する範囲の標準正規分布の確率変数』（qnorm(0.75) - qnorm(0.25) ≒ 1.3489）で割ったもの。↩
計算的には$IQR \times 0.7413$（$IQR_{N(0,1)}$の逆数）として実施することが多いみたいですが、今回は定義に沿って実装しています。↩
一応Argumentsにはそれっぽいことが書いてある。ヘルプはちゃんと読もう（戒め）。↩
整数どころかベクトルとかも大丈夫なので、行列xのそれぞれの行で別々の値を元に標準化したい場合はcenterとscaleにベクトル化して渡すとその通りに計算してくれる。たぶんtibbleデータに対してscale()を使った時と似たような動作じゃないかと思われる。あくまでも予想だけど。↩
Median Absolute Deviationの略なのだが、平均値を使って同様の処理をする平均絶対偏差もMean Absolute DeviationでMADと、かなりややこしい状況だったりする。なお、媒体によっては混同を避けるためにAAD（Average Absolute Deviation）としているが、Wikipediaの英語版には「AADは両方のMADを含みます」と書いてあったりして本当にややこしい。この辺りの話は後日書きます。↩
Rousseeuw,P.J.,andCroux,C.(1993).Alternativestothemedianabsolutedeviation,J.Amer.Statist.Assoc.88,1273-1283.↩
計算式が複雑だから言語やソフトによっては実装が難しいとか……？いやでもやれないことないでしょ……？みたいな感じで絶賛迷宮入り。誰か詳しい人いないですかね……。↩
MADを用いた標準化がこの方法でやっていて、かつ論文読んだ感じだとMADの代替（というか発展形）としてSnやQnが使えるとのことだったので、たぶん問題はない（はず）。間違ってたらすみません。↩

戯言日記

Rの話だと思ったら唐突にサバゲーが混じってくる何か。

R上でロバストな標準化をしたい