疫学と医療統計学と遺伝学と時々、大学院生活

疫学を専門とする大学院生の研究に関する備忘録的ページ。

R for beginners vol.3 「データのクリーニングと編集」

今回は先週のデータ入力や出力に続いて、実際の解析を行うにあたって必要なNAの除外や列の追加などのコマンドを学習する。
(今回で3回目であるが、この5回のコースは初心者Rユーザーのためのものであり、基本的な内容で構成されていることを再度確認しておきます。)

R for beginners vol.1 「Rの紹介と基本的なコマンド」
jojoshin.hatenablog.com

R for beginners vol.2 「データの入出力とサブセット」
jojoshin.hatenablog.com


【今回の狙い】Rでのデータインプット後のクリーニングやデータの追加、並び替え方法などの編集が実行できるようになる。

1. NAの除外をする。(is.naを利用する。)

データを収集する上で検出力の低下などの問題もあるので、NAはなるべくないように努力されるものであるが、何らかの原因でNAは含まれることがしばしばある。こうした場合にどう対処するかはそれぞれの判断による。まずはNAを見つけるコマンドを学習し、それを利用して様々な処理方法を考える。

x <- c(20, 1, NA, NA, 5, 40, NA, 67)

1-1. is.naを理解する。

#is.naは「NAなのはどれですか」という質問をした答えを返す。
is.na(x)
## [1] FALSE FALSE  TRUE  TRUE FALSE FALSE  TRUE FALSE

#!is.naは「NAではないのはどれですか」という質問をした答えを返します。
!is.na(x)
## [1]  TRUE  TRUE FALSE FALSE  TRUE  TRUE FALSE  TRUE

1-2. is.naを理解したら、それを利用してNAをどのように処理するかを学ぶ

#NAを0に置換する。replaceを使用する
x1 <- replace(x, is.na(x), 0)
x1
## [1] 20  1  0  0  5 40  0 67

#!is.naを利用して、xのうちNAでない要素のみを選択して、x2という新しいベクトルを作る
x2 <- x[!is.na(x)]
x2
## [1] 20  1  5 40 67

#以下の二つのコマンドも同義である
x3 <- x[which(is.na(x)==FALSE)]
x3
## [1] 20  1  5 40 67

x4 <- x[-which(is.na(x)==TRUE)]
x4
## [1] 20  1  5 40 67

1-3. 実際のデータセットでも使用してみる。

#NAを含んだdtNAを読み込む
dtNA <- read.csv("dtNA.csv", header=T, sep=",")

#dtNAは1000行8列のデータ
dim(dtNA)
## [1] 1000    8

#NAはあるか調べる
anyNA(dtNA)
## [1] TRUE

#欠損値(NA)の除外はna.omitでできる
dtNAomit <- na.omit(dtNA)
dim(dtNAomit)
## [1] 995   8

1'.番外編(NAを持つ人の特性をチェックしたい場合に有効。また、NA以外の値で抽出することも可能)

#どこにNAがあるか探し、その行番号を返す。
na <- which(is.na(dtNA$age)==TRUE | is.na(dtNA$gender)==TRUE | 
             is.na(dtNA$height)==TRUE | is.na(dtNA$weight)==TRUE | 
             is.na(dtNA$BMI)==TRUE | is.na(dtNA$alcohol)==TRUE | 
             is.na(dtNA$BMI01)==TRUE)
na
## [1]   5 661 664 948 997

length(na) #何個NAがあるかを意味する
## [1] 5

#どんな形でNAが入力されているかは確認する必要がある。
dtNA[which(is.na(dtNA$age)==TRUE | is.na(dtNA$gender)==TRUE | 
             is.na(dtNA$height)==TRUE | is.na(dtNA$weight)==TRUE | 
             is.na(dtNA$BMI)==TRUE | is.na(dtNA$alcohol)==TRUE | 
             is.na(dtNA$BMI01)==TRUE), ]
##     age gender   height    weight       BMI alcohol BMI01  ID
## 5    NA Female 152.3321  39.55690  17.04665     low     0   5
## 661  59   <NA> 137.8559  39.92813  21.01009  middle     0 661
## 664  54 Female 999.0000  50.53014 999.00000    high    NA 664
## 948  50 Female 157.7176 999.00000 999.00000     low    NA 948
## 997  44 Female 160.4450  49.39872  19.18948    <NA>     0 997

#naの5人を除外することにすると、、、
dtNAomit <- dtNA[-na,]
dim(dtNAomit)
## [1] 995   8
2. 特定の値や表記を見つけ、置換する。

一般的なデータセットの中には"M"、"men"、"男性"、"Man"などのように一つの性質を表すのに、いくつもの表記方法で記されているものがある。これを解析時には一つに置換したいと考えるのは自然であり、今回はその方法(grep)を説明する。

#性別の列にMale以外にもMenを含んだdtNAを読み込む
dtFIND <- read.csv("dtFIND.csv", header=T, sep=",")

#性別に関して、どんなカテゴリとその人数があるか見てみる
table(dtFIND$gender)

## Female   Male    Men 
##    580    411      9 

#Menの列を特定する
dtFIND[dtFIND$gender=="Men",]
##     age gender   height   weight      BMI alcohol BMI01  ID
## 22   53    Men 163.8752 65.11143 24.24548     low     0  22
## 40   42    Men 168.0086 67.38513 23.87268  middle     0  40
## 147  53    Men 164.5728 77.41877 28.58448     low     1 147
## 185  60    Men 165.4457 78.03305 28.50806  middle     1 185
## 341  47    Men 161.8126 48.53597 18.53699     low     0 341
## 402  61    Men 164.3409 90.20530 33.39958    high     1 402
## 564  39    Men 162.7939 79.00164 29.80985     low     1 564
## 805  42    Men 169.8075 68.57868 23.78348  middle     0 805
## 910  36    Men 174.3097 75.26493 24.77135  middle     0 910

#which関数やgrep関数も同義のコマンド
dtFIND[which(dtFIND$gender=="Men"),]
dtFIND[grep("Men", dtFIND$gender),]

#replaceを用いて、MenをMaleに置き換える
dtFIND$gender <- replace(dtFIND$gender, dtFIND$gender=="Men", "Male")

#再度、性別のカテゴリとその人数を見る
table(dtFIND$gender)

## Female   Male    Men 
##    580    420      0 
3. データの並び替え(sortとorder)
x <- c(0, 20, 40, 60, 50, 10, 30, 100)

3-1. sortはデフォルトでは昇順で「値」を並び替える。(降順の場合はdecreasing=Tを指定)

sort(x)
## [1]   0  10  20  30  40  50  60 100

sort(x, decreasing=T)
## [1] 100  60  50  40  30  20  10   0

3-2. orderはデフォルトでは昇順で「順序」を並び替える。(降順の場合はdecreasing=Tを指定)

order(x)
## [1] 1 6 2 7 3 5 4 8

order(x, decreasing=T)
## [1] 8 4 5 3 7 2 6 1

3-3. 実際のデータでも使用してみる。
例えば、dt.csvBMIの多い順(降順)にデータを並び替えて見る。

dt <- read.csv("dt.csv", header=T, sep=",")
dtBMI <- dt[order(dt$BMI, decreasing=T),]
head(dtBMI)
##     age gender   height    weight      BMI alcohol BMI01  ID
## 985  39   Male 165.7863 101.69633 37.00055    high     1 985
## 546  40   Male 160.3888  92.16509 35.82764  middle     1 546
## 566  33   Male 152.9183  82.26145 35.17851  middle     1 566
## 915  54   Male 160.8557  88.93664 34.37222    high     1 915
## 79   53   Male 154.8380  81.43232 33.96580    high     1  79
## 322  63   Male 160.4150  86.80829 33.73427  middle     1 322
4. データの結合(mergeとrbindとcbind)

4-1. mergeは指定のデータフレーム(2つ)を結合するコマンドである。
mergeで結合を宣言し、どのデータフレームとどのデータフレーム名を記述し、さらに結合のために使用する列を"by"で記述する。

#sample1では6人の3つの変数(ID, Age, Sex)を持つデータを生成する
ID <- c(1, 2, 3, 4, 5, 6)
Age <- c(34, 56, 50, 67, 40, 65)
Sex <- c("M", "M", "F", "F", "F", "M")
Sample1 <- data.frame(ID, Age, Sex)
Sample1
##   ID Age Sex
## 1  1  34   M
## 2  2  56   M
## 3  3  50   F
## 4  4  67   F
## 5  5  40   F
## 6  6  65   M

#sample2では7人の3つの変数(ID, BMI, TG)を持つデータを生成する
ID <- c(1, 2, 3, 4, 5, 6, 7)
BMI <- c(22.5, 18.5, 32, 20.9, 19.0, 25, 30)
TG <- c(100, 89, 109, 190, 123, 100, 86)
Sample2 <- data.frame(ID, BMI, TG)
Sample2
##   ID  BMI  TG
## 1  1 22.5 100
## 2  2 18.5  89
## 3  3 32.0 109
## 4  4 20.9 190
## 5  5 19.0 123
## 6  6 25.0 100
## 7  7 30.0  86

#byでマージする時の基準となる列を指定する
#パターンA(共通する列のみ残す場合)
SampleA <- merge(Sample1, Sample2, by="ID")
SampleA
##   ID Age Sex  BMI  TG
## 1  1  34   M 22.5 100
## 2  2  56   M 18.5  89
## 3  3  50   F 32.0 109
## 4  4  67   F 20.9 190
## 5  5  40   F 19.0 123
## 6  6  65   M 25.0 100

#パターンB(全ての列を残す場合はall=Tを指定する)
SampleB <- merge(Sample1, Sample2, by="ID", all=T)
SampleB
##   ID Age  Sex  BMI  TG
## 1  1  34    M 22.5 100
## 2  2  56    M 18.5  89
## 3  3  50    F 32.0 109
## 4  4  67    F 20.9 190
## 5  5  40    F 19.0 123
## 6  6  65    M 25.0 100
## 7  7  NA <NA> 30.0  86

4-1. rbindは「行」を追加する

a <- c(1, 2, 3, 4, 5)
b <- c(6, 7, 8, 9, 10)
rbind(a, b)
##   [,1] [,2] [,3] [,4] [,5]
## a    1    2    3    4    5
## b    6    7    8    9   10

rbindは3つでも良い。

a <- c(1, 2, 3, 4, 5)
b <- c(6, 7, 8, 9, 10)
c <- c(11, 12, 13, 14, 15)
rbind(a, b, c)
##   [,1] [,2] [,3] [,4] [,5]
## a    1    2    3    4    5
## b    6    7    8    9   10
## c   11   12   13   14   15

4-2. cbindは「列」を追加する

a <- c(1, 2, 3, 4, 5)
b <- c(6, 7, 8, 9, 10)
cbind(a, b)
##      a  b
## [1,] 1  6
## [2,] 2  7
## [3,] 3  8
## [4,] 4  9
## [5,] 5 10

cbindも3つでも良い。

a <- c(1, 2, 3, 4, 5)
b <- c(6, 7, 8, 9, 10)
c <- c(11, 12, 13, 14, 15)
cbind(a, b, c)
     a  b  c
## [1,] 1  6 11
## [2,] 2  7 12
## [3,] 3  8 13
## [4,] 4  9 14
## [5,] 5 10 15

4-3. 実際のデータでもrbindを使用してみる。(例えば、1人のデータをさらに結合したい場合)

d <- c(34, "Male", 157.980, 45.998, 45.998/1.57980/1.57980, "low", 0, 1001)
dt1 <- rbind(dt, d)
tail(dt1)
##      age gender           height           weight              BMI alcohol
## 996   47   Male  162.07143760552 73.0088145070084 27.7947215591229  middle
## 997   44 Female 160.445023328181 49.3987240121513 19.1894810303735     low
## 998   41 Female 150.707270295575 44.6205968014313  19.645675547013     low
## 999   45 Female 155.537203329995 48.3571241587341 19.9890379975865     low
## 1000  39 Female 155.848962643407 44.0017588730915 18.1159955371478  middle
## 1001  34   Male           157.98           45.998 18.4303986840059     low
##      BMI01   ID
## 996      1  996
## 997      0  997
## 998      0  998
## 999      0  999
## 1000     0 1000
## 1001     0 1001

4-4. 実際のデータでもcbindを使用してみる。(例えば、新しい検査項目のデータを結合したい場合)

MBP <- rnorm(1000, mean=120, sd=12)
dt2 <- cbind(dt, MBP)
head(dt2)
##   age gender   height   weight      BMI alcohol BMI01 ID      MBP
## 1  35   Male 170.5088 73.51593 25.28644  middle     1  1 108.6553
## 2  45 Female 154.8704 41.49381 17.30001     low     0  2 126.3991
## 3  41 Female 152.4852 43.05607 18.51736     low     0  3 131.1846
## 4  44 Female 148.0986 40.32888 18.38713  middle     0  4 122.3526
## 5  68 Female 152.3321 39.55690 17.04665     low     0  5 120.6014
## 6  42   Male 174.9575 76.66370 25.04522  middle     1  6 127.2838

4-5. 新しい列を作る。年齢が50歳以上の人を1として、それ以外を0とする新しい列「age01」を作成する。

#シンプルな方法
#dt$age > 50でTRUEorFALSEの列を作り、それを数値化することで01データになる。それをdt$age01と新しく名付けることでデータフレームに新たな列を追加できる
dt$age01 <- as.numeric(dt$age > 50)
head(dt$age01, 30)
## [1] 0 0 0 0 1 0 1 1 0 1 0 0 0 0 0 1 0 1 0 0 1 1 1 0 0 1 1 1 0 1

#複雑な方法(無視してもらって構わない)
dt$age01 <- NA
for (i in 1:1000)
{
  if (dt$age[i] < 50){
  dt$age01[i] <- 0
} else {
  dt$age01[i] <- 1
}
}
5. データのカットとラベル化
#1〜100までの整数を生成する
x <- 1:100

#xを80でカットし、cxとする。この場合は0<x<=80、80<x<=100という二つのグループに分割している
cx <- cut(x, breaks=c(0, 80, 100))
table(cx)
## cx
##   (0,80] (80,100] 
##       80       20 

#同様にxを80でカットし、cx1とする。cx1のラベルを数字にする。
cx1 <- cut(x, breaks=c(0, 80, 100), labels=FALSE)
table(cx1)
## cx1
##  1  2 
## 80 20 
 
#同様にを80でカットし、cx2とし、それぞれのグループにラベル(miss, pass)をつける
cx2 <- cut(x, breaks=c(0, 80, 100), labels=c("miss", "pass"))
table(cx2)
## cx2
## miss pass 
##   80   20 

 #同様にを80でカットし、cx3とするが、今回は0<=x<80、80<=x<100という二つのグループに分割している。範囲の右側は含まれないという意味で「right=F」のコマンドが入っている。つまり、グループ1には1〜79の値をとる79人、グループ2には100を除く80〜99の値をとる20人が含まれる。
cx3 <- cut(x, breaks=c(0, 80, 100), right=F, labels=FALSE)
table(cx3)
## cx3
##  1  2 
## 79 20 
 
#cx3の問題としては100の値を持つ人が漏れていたことであるが、これは「include.lowest=T」で解決できる
cx4 <- cut(x, breaks=c(0, 80, 100), right=F, labels=FALSE, include.lowest=T)
table(cx4)
## cx4
##  1  2 
## 79 21 
 
#単純に、xを3分割して(breaks=3と指定)、cx5とする(これは分割したい変数が等しく分布している場合に有効)
cx5 <- cut(x, breaks=3, labels=c(1, 2, 3))
table(cx5)
## cx5
##  1  2  3 
## 34 33 33 

5-1. 三分位を指定して、分割する(dtのBMIを三分割して、BMI012というカテゴリを作る)

#BMIの0%(min)と33%, 67%, 100%(max)を表示する
c(quantile(dt$BMI, .0), quantile(dt$BMI,.333), quantile(dt$BMI,.667), quantile(dt$BMI, 1.0))

##       0%    33.3%    66.7%     100% 
## 12.06728 18.33456 22.40179 37.00055 

#上記のポイント設定し、BMIを3つに分割する
dt$BMI012 <- cut(dt$BMI, 
                 breaks=c(0, quantile(dt$BMI,.333), quantile(dt$BMI,.667), max(dt$BMI)),
                 labels=c(0, 1, 2))

#結果を表にする
table(dt$BMI012)
## 
##   0   1   2 
## 333 334 333 

#生成した結果と元のデータを比較する
head(dt[,c("BMI", "BMI012")])
##        BMI BMI012
## 1 25.28644      2
## 2 17.30001      0
## 3 18.51736      1
## 4 18.38713      1
## 5 17.04665      0
## 6 25.04522      2
練習問題(9問)

まずはじめに、配布したsampleA.csv(調査票データ)とsampleB.csv(身体計測および生化学データ)を読み込む。(それぞれR内では"sampleA"と"sampleB"と名付ける)

1. それぞれのサンプルデータにNAがあるか確認し、ある場合はそれらを除外し、"sampleA1"と"sampleB1"と名付けよ.また、それぞれのデータの対象人数を答えよ.
2. 1で作成したそれぞれのデータを"ID"という列名を使用して、結合し、"sample"という新しいデータを作成せよ、またこの段階での対象者の数を答えよ.
3. 2で作成したデータに下記に示す新しいアルブミンのデータ(Alb)を追加せよ.
alb <- round(rnorm(nrow(sample), 4.2, 0.2), 1)
4. BMI > 25の人を1とする(それ以外は0となる)新しい列'BMI01'を作成し、それぞれの階級に何人いるか示せ.(dt$BMI01のように列を作成する)
5. BMIを3等分した新しい列'BMI012'を作成し、それぞれの階級に何人いるか示せ.
6. ここまでの過程で完成したデータ"sample"を"sample.csv"として書き出せ.
7. 男女が何人いるか示せ.(復習)
8. 男性でBMIが22以上の人は何人いるか示せ.(復習)
9. 最高血圧が130mmHg以上かつ/もしくは最低血圧が85mmHg以上の人は何人いるか示せ.(復習)

今回は以上である。
次回は6/2までに更新する(次回は簡単な統計検定量および図示についてにする。)

20160526
RF