データを扱う中で「外れ値を外す方法」に関して検討してみます。

apply関数の記事で以下のコメントがありました。




大変勉強になりました。
列ごとに、ファイルの値(NAがある)が2SDより高い値をNAにしたいです。
Data01 <-read_excel("raw.xlsx",1)
Data02 <-apply(Data01, 2, sd, na.rm = TRUE)
Data01[Data01 > 2*Data02]<- NA
以上のコードをして最後は列ごとにうまく出来なかったです。
ご提言があれば大変嬉しいです。

自分もそうですが、実際にやってみるとうまくいかないことってたくさんあります。
今回は質問を一緒に考えながらどう考えるか?どんな解決方法があるかを考えて、今までの理解を深めるきっかけになればと思います。

1.データの作成

コードを確認してみます。
Data01 <-read_excel("raw.xlsx",1)
これはraw.xlsxの1シート目を読み込むという意味なので、今回は仮のデータを作ります。

set.seed(2020)
a <- c(rnorm(998,10,1),NA,30)
b <- c(rnorm(998,20,1),NA,40)
c <- c(rnorm(998,30,1),NA,50)
Data01 <- data.frame(a,b,c)
head(Data01)
簡単に言うと以下のとおりです。
a : 平均10, 標準偏差1のデータ998個にNAと30を加えた1000個のデータ
b : 平均20, 標準偏差1のデータ998個にNAと40を加えた1000個のデータ
c : 平均30, 標準偏差1のデータ998個にNAと50を加えた1000個のデータ

スクリーンショット 2020-03-28 2.34.45
この赤色が2SDを超えたデータです。やりたいことはこの赤い点のデータを削除することです(これの是非については今回問いません)

2.apply関数のおさらい

次の行を確認します。
Data02 <-apply(Data01, 2, sd, na.rm = TRUE)
apply関数【2-1】Rのfor関数、apply関数を使ってまとめて標準偏差などの統計量を求める方法で紹介しました。下はイメージ図です。
スクリーンショット 2020-03-27 13.46.36
上図ではsd関数(標準偏差)を求めました。
また欠損値(NA)があるのでrm.na = TRUEをつけることでNAの値を外して計算してくれます。

スクリーンショット 2020-03-27 13.54.01
それぞれの列の標準偏差が表示されました。
ただ元データが2標準偏差を超えているかどうかを確認するにはこれだけでは足りません。


3.どう考えるか?
スクリーンショット 2020-03-28 4.00.30

それぞれのデータが2SDを超えているかどうかを判定する必要があります。
具体的には(平均-2SD) < それぞれの値 < (平均 + 2SD)の条件を満たしているかどうかとなります。

そのため以下の2つを行えばなんとかできそうです。
①a,b,cそれぞれの平均±2標準偏差の値を求める
②条件に合えばそのまま、条件にあわなければNAと変換する


 4.平均±2SDを求める

まずa,b,cの平均と標準偏差を求めます。
これには前述のapply関数かsummarize関数を使います。
まずはapply関数でやってみます。

Data01_mean <-apply(Data01, 2, mean, na.rm = TRUE)
Data01_sd <-apply(Data01, 2, sd, na.rm = TRUE)
Data01_mean
Data01_sd

low <- Data01_mean - 2*Data01_sd
high <- Data01_mean + 2*Data01_sd
low
high
スクリーンショット 2020-03-28 3.17.41

ここからlow[2]、もしくはlow["b"]とすればbの平均-2SDが取り出せます。
high[1]、もしくはhigh["a"]であればaの平均+2SDとなります。

5.条件に合わないものをNAにする

次に条件に合わないものを抽出しNAにします。
データを抽出するにはデータ[条件]を使います

(aのデータ < 平均-2sd) または (平均+2SD < aのデータ)だとこうなります
Data01$a[(Data01$a < low["a"]) | (high["a"] < Data01$a)]
スクリーンショット 2020-03-28 4.08.45

条件式の記号については以下をご参照ください。今回はまたはの意味の を使いました。
スクリーンショット 2019-02-11 2.04.40
抽出した条件に <- を使うことでデータを置き換えることができます。
Data02 <- Data01
Data02$a[(Data02$a < low["a"]) | (high["a"] < Data02$a)] <- NA
Data02$b[(Data02$b < low["b"]) | (high["b"] < Data02$b)] <- NA
Data02$c[(Data02$c < low["c"]) | (high["c"] < Data02$c)] <- NA
これで完成です!
ただ他にも方法があります。どれもよく使う方法なので、もし余力があればぜひ比べてみてください。


6.tidyverseパッケージを使う(dplyrパッケージでも同じ)

第2章で紹介しているtidyverseパッケージのsummarize関数とmutate関数を使うこともできます。
tidyverseパッケージを使うときには先にlibrary関数で呼び出します。なければ先にインストールします。
#もし初めての場合はインストール
install.packages("tidyverse")
#インストール済ならlibrary関数で呼び出す
library(tidyverse)
low_high <-
  summarize(Data01,
          a_low = mean(a, na.rm = TRUE)-2*sd(a, na.rm = TRUE),
          a_high = mean(a, na.rm = TRUE)+2*sd(a, na.rm = TRUE),
          b_low = mean(b, na.rm = TRUE)-2*sd(b, na.rm = TRUE),
          b_high = mean(b, na.rm = TRUE)+2*sd(b, na.rm = TRUE),
          c_low = mean(c, na.rm = TRUE)-2*sd(c, na.rm = TRUE),
          c_high = mean(c, na.rm = TRUE)+2*sd(c, na.rm = TRUE)) 
class(low_high) #型を確認するとdata.frameとなっている
low_high <- unlist(low_high)
class(low_high) #型を確認するとnumeric(数値型)となっている
low_high
スクリーンショット 2020-03-28 3.37.43
文字は多いですが、コピペを使えば作業自体は楽になります。
ただsummarize関数の結果は表形式(data.frame型)で、このままでは数値として計算できません。
そのためunlist関数を使って数値型に戻します。
low_high[5]、もしくはlow_high["c_low"]とすればcの平均-2SDとなります。


既に第2章を読まれた方は%>%を使う事もできます。
#%>%使う場合
low_high <- summarize(Data01, a_low = mean(a, na.rm = TRUE)-2*sd(a, na.rm = TRUE), a_high = mean(a, na.rm = TRUE)+2*sd(a, na.rm = TRUE), b_low = mean(b, na.rm = TRUE)-2*sd(b, na.rm = TRUE), b_high = mean(b, na.rm = TRUE)+2*sd(b, na.rm = TRUE), c_low = mean(c, na.rm = TRUE)-2*sd(c, na.rm = TRUE), c_high = mean(c, na.rm = TRUE)+2*sd(c, na.rm = TRUE)) %>%
unlist()
low_high
スクリーンショット 2020-03-28 3.32.12

次に条件に応じて分類するifelse関数を使うこともできます。

ifelse(条件TRUEの場合FALSEの場合)

列の追加や修正を行うmutate関数を使うことでこうなります。
Data02 <- 
  Data01 %>% 
  mutate(a = ifelse((low_high["a_low"] < a) & (a < low_high["a_high"]), a, NA),
         b = ifelse((low_high["b_low"] < b) & (b < low_high["b_high"]), b, NA),
         c = ifelse((low_high["c_low"] < c) & (c < low_high["c_high"]), c, NA))
スクリーンショット 2020-03-28 4.57.19

条件式を少し変えているので確認してみてください。今回は&を使っています。

tidyverseパッケージは現在のRを使う上でとても重要になります。
ifelse,mutate関数に関してはこちらで紹介していますので難しかったという方はご参照ください。




7.関数を自作しapply関数でまとめて処理する
他の方法としては関数を自作するという方法があります。

関数名 <- function(○) {中身のコード}

今回はdrop_2sdという名前の関数を作ります(名前は何でもいいけど文字の最初に数字を使うのは×)。


a = ifelse((aの平均-2SD < a) & (a < aの平均+2SD), a, NA),
b = ifelse((bの平均-2SD < b) & (b < bの平均+2SD), b, NA),
c = ifelse((cの平均-2SD < c) & (c < cの平均+2SD), c, NA)
 ↓
function(x){
 ifelse((xの平均-2SD < x) & (x < xの平均+2SD), x, NA)
}
と共通部分をxにします。ちなみにxでなくてもcolとか別の名前をつけても大丈夫です。

drop_2sd <- function(x){
  ifelse((mean(x, na.rm = TRUE) - 2 * sd(x, na.rm = TRUE) < x) & (x < mean(x, na.rm = TRUE) + 2 * sd(x, na.rm = TRUE)), x, NA)
}
そしてapply関数で各列にdrop_2sd関数を使います。
Data02 <- apply(Data01, 2, drop_2sd)
スクリーンショット 2020-03-28 4.50.32
2行のコードで終わりました。
関数を作るのは最初は慣れないですが、使えるようになると強力です。


8.まとめ

今回は標準偏差±2SDの除去を3つの方法で紹介しました。

本やサイトで独学すると1つの方法でしか紹介してないことが多いです。

本当は自分の知っている知識でもできるはずなのに、知らない方法を偶然読んでしまいドツボにはまることもあります。そういった紆余曲折は勉強する上で必要な苦労でもありますが、いくつかの方法を見比べることでショートカットできることもあるのではと思っています。

また自分が躓いたところや他の記事でもコメントがあれば検討していきますので、ブログのコメントやtwitterで連絡いただければ幸いです。