Rのマハラノビス()関数は、多次元データの外れ値を検出する簡単な手段を提供する。
例えば、身長と体重のデータフレームがあるとする。
hw <- data.frame(Height.cm=c(164, 167, 168, 169, 169, 170, 170, 170, 171, 172, 172, 173, 173, 175, 176, 178),
Weight.kg=c( 54, 57, 58, 60, 61, 60, 61, 62, 62, 64, 62, 62, 64, 56, 66, 70))
これらのデータ（この例ではインタラクティブプロットを使って生成）をプロットするとき、例えば、平均身長や平均体重から2（サンプル）以上の標準偏差がある点を外れ値としてマークすることができる。
is.height.outlier <- abs(scale(hw$Height.cm)) > 2 is.weight.outlier <- abs(scale(hw$Weight.kg)) > 2 pch <- (is.height.outlier | is.weight.outlier) * 16
グラフの描画。
plot(hw, pch=pch)
なお、身長が175cmの点（グラフ右下）は、身長と体重の平均値から2標準偏差未満であるため、異常値としてマークされていない。
しかし、この点は、このデータで見られる身長と体重の直線的な関係に最も明確に従わない点である。これは間違いなく、本当の外れ値である。
ところで、上のグラフのスケールの選択は、やや誤解を招きやすい。少なくともyスケールをゼロからスタートさせれば、身長の体重への影響をより明確に把握することができたはずである。
しかし、このデータは単に外れ値検出の説明のために使っているだけなので、この悪い習慣は見逃す。
次に、mahalanobis()を使ってみよう。
n.outliers <- 2 # 最も極端な2点を外れ値としてマークする m.dist.order <- order(mahalanobis(hw, colMeans(hw), cov(hw)), decreasing=TRUE) is.outlier <- rep(FALSE, nrow(hw)) is.outlier[m.dist.order[1:n.outliers]] <- TRUE pch <- is.outlier * 16
上記のコードでは、マハラノビス距離（一般化二乗距離とも呼ばれる）に応じて、最も極端な2点を外れ値としてマークしている。
これは、非常に大雑把に言うと、データフレームが構成するデータの中心からの各点（データフレームの行）の距離を、各変数（データフレームの列）の標準偏差で正規化し、それらの変数の共分散で調整したものである。
t 統計量は、(x-x) / s で定義され、x は標本の平均値、s は標本の標準偏差であり、多次元統計量のアナログと考えることができる。
見ての通り、今回はグラフの右下隅の点が捕捉されている。
plot(hw, pch=pch)
そして、このテクニックは高次元でも有効である。このコードは3次元のスピン可能な散布図を作成する。
library(rgl)
hwa <- data.frame(Height.cm=c(164, 167, 168, 168, 169, 169, 169, 170, 172, 173, 175, 176, 178),
Weight.kg=c( 55, 57, 58, 56, 57, 61, 61, 61, 64, 62, 56, 66, 70),
Age =c( 13, 12, 14, 17, 15, 14, 16, 16, 13, 15, 16, 14, 16))
m.dist.order <- order(mahalanobis(hwa, colMeans(hwa), cov(hwa)), decreasing=TRUE)
is.outlier <- rep(FALSE, nrow(hwa))
is.outlier[m.dist.order[1:2]] <- TRUE # Mark as outliers the 2 most extreme points
col <- is.outlier + 1
散布図は以下のような感じである。赤い点が外れ値だ。
plot3d(hwa$Height.cm, hwa$Weight.kg, hwa$Age, type="s", col=col)
上記のコードからわかるように、mahalanobis()関数は、与えられた平均のベクトルと与えられた共分散行列を使って、データフレームのマハラノビス距離を計算する。
一般的には、上記の例のように、対象となるデータフレームから計算された平均のベクトルと共分散行列を渡して使用したいと思うだろう。例えば、
# m.dist <- mahalanobis(my.dataframe, colMeans(my.dataframe), cov(my.dataframe))
得られた距離のベクトルは、データフレームから最も極端な行を除外するために使用することができる。例えば、最も極端な5%のポイントを削除したい場合がある。
# percentage.to.remove <- 5 # 5%のポイントを削除 # number.to.remove <- trunc(nrow(my.dataframe) * percentage.to.remove / 100) # m.dist <- mahalanobis(my.dataframe, colMeans(my.dataframe), cov(my.dataframe)) # m.dist.order <- order(m.dist, decreasing=TRUE) # rows.to.keep.index <- m.dist.order[(number.to.remove+1):nrow(my.dataframe)] # my.dataframe <- my.dataframe[rows.to.keep.index,]
これは、実行中の分析が極端な点の影響を過剰に受けていないかどうかを素早くチェックしたい場合によく使われる。
まず、全データセットで分析を実行し、次に上記の手法で最も極端な点を削除する。
結果に大きな差がある場合は、外れ値に対してよりロバストな分析結果を使用することを検討するとよいだろう。
注意事項
マハラノビス距離の式は線形関係しか考慮しないので、データが非線形関係を示す場合は mahalanobis()に注意する必要がある。例えば、以下のようなコードを実行してみよう。
my.dataframe <- data.frame(x=c(4, 8, 10, 16, 17, 22, 27, 33, 38, 40, 47, 48, 53, 55, 63, 71, 76, 85, 85, 92, 96),
y=c(6, 22, 32, 34, 42, 51, 59, 63, 64, 69, 70, 20, 70, 63, 63, 55, 46, 41, 33, 19, 6))
percentage.of.outliers <- 10 # Mark 10% of points as outliers
number.of.outliers <- trunc(nrow(my.dataframe) * percentage.of.outliers / 100)
m.dist <- mahalanobis(my.dataframe, colMeans(my.dataframe), cov(my.dataframe))
m.dist.order <- order(m.dist, decreasing=TRUE)
rows.not.outliers <- m.dist.order[(number.of.outliers+1):nrow(my.dataframe)]
is.outlier <- rep(TRUE, nrow(my.dataframe))
is.outlier[rows.not.outliers] <- FALSE
pch <- is.outlier * 16
グラフ描画。
plot(x=my.dataframe$x, y=my.dataframe$y, pch=pch)
なお、検討中のデータセットの変数間の関係は非線形であるため、最も明白な異常値は検出されていない。
また、マハラノビス距離を用いて多次元データの極値を求めると、t統計量を用いて一次元データの極値を求めるときと同様の問題が生じる可能性があることに注意する必要がある。
1次元の場合、極値が平均値や標準偏差に強い影響を与えるため、t統計量は中心からの距離の合理的な指標にはならなくなるのです。(そのため、中央値や中央値絶対偏差など、外れ値に強い中心値や距離の指標を使った方がよい場合が多いのである)。
マハラノビス距離もt統計量のように、標準偏差と平均からの距離から計算されるので、多次元の場合にも同じことが起こりえる。
