タグ:apply

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

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で連絡いただければ幸いです。

前回は今までの総復習を行いました。
【演習】R初心者が統計をかけるための前準備の流れを復習します : 独学で始める統計×データサイエンス



上記の記事を進めると以下のように統計量をまとめて出すことができますが、summary関数では標準偏差が出ません。
スクリーンショット 2019-02-04 21.49.17


今回は標準偏差を出すいろいろな方法を復習も含め紹介します。

繰り返しのfor関数やapply関数、パイプ演算子(%>%)もはじめて出てきますが、使えるようになるとFIMのような項目が多いものでも1度にまとめて解析できるようになり、いよいよExcelではできない体験をすることができるようになります。

まだ今回の内容であればExcelでも可能ですが、1つずつ進めていきましょう。

1.sd関数で1つずつ出す
2.for関数で繰り返す
3.apply関数を使う
4.tableoneパッケージを使う


今回は前回の記事のデータを使います。

data01.xlsx 

ダウンロードした後、プロジェクトの指定フォルダにファイルを移動させておけば、以下のコマンドで上記の画面まで進みます。

library(readxl)
data01 <- read_excel("data01.xlsx", sheet = "日本語")
data01$性別 <- factor(data01$性別, levels = c("男性","女性"))
data01$歩行 <- as.factor(data01$歩行)
summary(data01)

*注意
もしプロジェクトファイルを作成していない場合はgetwd()と打ち込むと作業フォルダが表示されます、作業フォルダにExcelデータを入れると後は同じです。

getwd()




1.sd関数で1つずつ出す

標準偏差を出す関数はsd関数です。

sd(ベクトル)
sd(data01[[3]])
sd(data01[[4]])
sd(data01[[5]])
sd(data01[[6]])
sd(data01[[7]])
sd(data01[[8]])
スクリーンショット 2019-02-04 21.53.36


data.frameからベクトルを表示するには$または[[ ]]を使います。
[[ ]]の中の数字は列番号です。
$でもいいのですが、このようなケースでは列番号のほうが手っ取り早いです。

列番号を確認するにはnames関数t関数(もしくはdata.frame関数)を使うと便利です。

t(names(data01))
スクリーンショット 2019-02-04 23.20.52

data.frame(names(data01))
スクリーンショット 2019-02-04 23.22.06


2.for関数で繰り返す

1個ずつ出すのもいいですが、繰り返しになる作業はfor関数で繰り返すことができます。
for(i in 3:8){
  x <- sd(data01[[i]])
  print(x)
  }

for関数は(◯ in △)と{}の2つに別れています。
スクリーンショット 2019-02-04 23.51.05
スクリーンショット 2019-02-04 23.59.46


スクリーンショット 2019-02-05 0.03.59

ただ結果は出るのですが結果の数字しか出ません。


3-1.apply関数を使う

library(tidyverse)
apply(data01, 2, sd)

スクリーンショット 2019-02-05 0.17.43

スクリーンショット 2019-02-05 0.08.54


下にエラーが出ます。これは「氏名と性別で集計できないからNAにしてますよ!」と怒られています。なのでselect関数で集計する列だけ抜き出します。


library(tidyverse) #既に実行していたらなくても可
data01_select <- select(data01, 3:8) data01_select apply(data01_select, 2, sd)
スクリーンショット 2019-02-05 0.32.27

今回はdata01_selectという変数名にselect関数で3〜8列目だけ抜き出しapply関数を使いました。


スクリーンショット 2019-02-05 0.33.14


ただこの方法はdata01_selectのように間に変数を挟む必要があります。
tidyverseパッケージはこれをスッキリさせるパイプ演算子(%>%)というものがあります。


3−2.パイプ演算子を使ったapply関数

パイプ演算子(%>%)は解析の流れをスムーズにしてくれます。

データ① %>%
    プログラム② %>%
    プログラム③ %>%
    プログラム④

とすると先程のdata01_selectの用な変数を介することもなく、データ①をプログラム②に送り、その結果をプログラム③に送り、その結果をプログラム④に送るといったことが可能になります

library(tidyverse) #既に実行していたらなくても可
data01 %>% #data01に対して
  select(3:8) %>% #3列目(年齢)〜8列目(MMSE)までを選択し
  apply(., 2, sd) #それぞれの列で標準偏差を求めて!

スクリーンショット 2019-02-05 0.55.55
スクリーンショット 2019-02-05 0.57.27


もし四捨五入するならround関数を使います。

library(tidyverse)#既に実行していたらなくても可
data01 %>% #data01に対して
  select(3:8) %>% #3列目(年齢)〜8列目(MMSE)までを選択し
  apply(., 2, sd) %>% #それぞれの列で標準偏差を求めて!
  round(.,2) #小数点2位より下を四捨五入する
スクリーンショット 2019-02-05 1.13.08


4.tableoneパッケージを使う

tableoneパッケージを使うのも1つの方法です。

tableoneパッケージについてはこの記事で紹介しています。
Rで医療統計で必要なtable1を作るtableoneパッケージについて紹介します : 独学で始める統計×データサイエンス

library(tableone)
CreateTableOne(data = data01,
               vars = c("年齢", "身長", "体重", "SIAS", "BBS", "MMSE"))
スクリーンショット 2019-02-05 1.20.05


まとめ

今回は標準偏差を出すいろいろな方法を通じてfor関数やapply関数、パイプ演算子の使い方を紹介しました。

今回は標準偏差(sd)でしたが、平均(mean)や中央値(median)など他の統計量も求めることができます。

for関数やtidyverse(dplyrなど)に関してはまとまった記事もありますが、これからもちょくちょく出てくると思います。その都度復習しながら使い方のイメージを付けていただければいいかと思います。








↑このページのトップヘ