読者です 読者をやめる 読者になる 読者になる

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

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

R for beginners vol.4 「データの要約とビジュアライゼーション」

今回はデータの整頓から少し分析に近いことを始めます。その中でデータを要約し、図示することが解析の一歩かと思い今回の内容にしています。

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

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

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

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

1-1.【復習】基本的な統計量の算出方法(Class1で実行済みなのでスキップする)

x <- c(100, 210, 114, 340, 45, 60, 85, 229, 91, 156)
mean(x)
## [1] 143
median(x)
## [1] 107
min(x)
## [1] 45
max(x)
## [1] 340
sd(x)
## [1] 91.96255
sum(x)
## [1] 1430
1-2. 実際のデータセットで実行してみる
sample <- read.csv("sample.csv", sep=",", header=T)
1-3. 要約統計量を算出するのに有効なsummary関数(但し、sdは表示されないことに注意)
#sample内の全ての変数について要約統計量を求める
summary(sample)
##       ID              age            sex         alcohol       exercise       smoking         height          weight           BMI          bodyfat           sbp       
##  Min.   :   1.0   Min.   :30.00   Female:1790   high  : 699   high  : 297   current: 248   Min.   :132.3   Min.   :23.70   Min.   :10.0   Min.   :10.40   Min.   : 52.0  
##  1st Qu.: 751.5   1st Qu.:42.00   Male  :1193   low   :1292   low   :1294   ever   : 946   1st Qu.:151.0   1st Qu.:50.80   1st Qu.:20.8   1st Qu.:23.30   1st Qu.:117.0  
##  Median :1508.0   Median :50.00                 middle: 992   middle:1392   never  :1789   Median :157.0   Median :59.30   Median :23.7   Median :27.20   Median :130.0  
##  Mean   :1503.6   Mean   :49.96                                                            Mean   :157.7   Mean   :58.95   Mean   :23.7   Mean   :27.98   Mean   :130.6  
##  3rd Qu.:2254.5   3rd Qu.:58.00                                                            3rd Qu.:164.0   3rd Qu.:67.10   3rd Qu.:26.5   3rd Qu.:32.40   3rd Qu.:145.0  
##  Max.   :3000.0   Max.   :68.00                                                            Max.   :184.3   Max.   :98.00   Max.   :41.1   Max.   :51.40   Max.   :198.0  
##       dbp               bs               TC            TG             HDLC             LDLC            cre              alb            BMI01            BMI012      
##  Min.   : 31.00   Min.   : 32.00   Min.   : 84   Min.   : 17.0   Min.   : 14.00   Min.   :  8.0   Min.   :0.1800   Min.   :3.500   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.: 67.00   1st Qu.: 79.00   1st Qu.:189   1st Qu.: 62.0   1st Qu.: 50.00   1st Qu.:104.0   1st Qu.:0.6300   1st Qu.:4.100   1st Qu.:0.0000   1st Qu.:0.0000  
##  Median : 76.00   Median : 89.00   Median :213   Median : 89.0   Median : 59.00   Median :125.0   Median :0.7200   Median :4.200   Median :0.0000   Median :1.0000  
##  Mean   : 75.81   Mean   : 89.02   Mean   :212   Mean   :104.8   Mean   : 59.33   Mean   :125.1   Mean   :0.7521   Mean   :4.201   Mean   :0.3667   Mean   :0.9893  
##  3rd Qu.: 85.00   3rd Qu.: 99.00   3rd Qu.:236   3rd Qu.:129.0   3rd Qu.: 68.00   3rd Qu.:146.0   3rd Qu.:0.8500   3rd Qu.:4.300   3rd Qu.:1.0000   3rd Qu.:2.0000  
##  Max.   :125.00   Max.   :144.00   Max.   :324   Max.   :694.0   Max.   :108.00   Max.   :221.0   Max.   :1.4900   Max.   :4.800   Max.   :1.0000   Max.   :2.0000 

#sample全部ではなく、幾つかの列を指定する場合もある
summary(sample[, c(2,3,7)])
##        age            sex           height     
##  Min.   :30.00   Female:1790   Min.   :132.3  
##  1st Qu.:42.00   Male  :1193   1st Qu.:151.0  
##  Median :50.00                 Median :157.0  
##  Mean   :49.96                 Mean   :157.7  
##  3rd Qu.:58.00                 3rd Qu.:164.0  
##  Max.   :68.00                 Max.   :184.3  

#上のコマンドと同義だが、列名を指定しても良い
summary(sample[,c("age", "sex", "height")])
##        age            sex           height     
##  Min.   :30.00   Female:1790   Min.   :132.3  
##  1st Qu.:42.00   Male  :1193   1st Qu.:151.0  
##  Median :50.00                 Median :157.0  
##  Mean   :49.96                 Mean   :157.7  
##  3rd Qu.:58.00                 3rd Qu.:164.0  
##  Max.   :68.00                 Max.   :184.3  
あれ、BMI01とBMI012が数値データになってない?
class(sample$BMI01)
## [1] "integer"
class(sample$BMI012)
## [1] "integer"

#変更して、もう一度確認しよう。
sample$BMI01 <- factor(sample$BMI01)
sample$BMI012 <- factor(sample$BMI012)
summary(sample[,c("BMI01","BMI012")])
## BMI01    BMI012  
## 0:1889   0:1007  
## 1:1094   1:1001  
##          2: 975  

2. Apply系

非常に便利なコマンドであり、それぞれの用途にあったapply系の関数を使用するとデータの要約が上手くいく。(Apply系の概念図は参考資料[1]をご参照ください)

2-1. apply(x)

具体例で説明すると、下の例では、apply関数をsampleの7,8,9,10列目に適応して、2で列ごとに計算するように指定する。そのあとにmeanやsdなどを指定する。(ちなみに1を指定すると、行ごとの値をまとめる)(na.rm=TRUEでnaは除いて計算するように指定する)

apply(sample[,c(7,8,9,10)], 2, mean, na.rm=TRUE)
##    height    weight       BMI   bodyfat 
## 157.68042  58.94841  23.70104  27.97875 
apply(sample[,c(7,8,9,10)], 2, sd, na.rm=TRUE)
##    height    weight       BMI   bodyfat 
##  8.699511 11.736069  4.359196  6.475639 
2-2. tapply(x)

これはおそらく使用頻度の高い関数です。ある値(BMI)をあるグループ(alcohol)ごとに統計量(mean, sd etc...)を表示したい。一元配置分散分析の時に平均値や中央値を比較するのに便利です。

mean <- tapply(sample$BMI, sample$alcohol, mean, na.rm=TRUE)
sd <- tapply(sample$BMI, sample$alcohol, sd, na.rm=TRUE)
mean
##     high      low   middle 
## 24.33219 23.31161 23.76351
sd
##     high      low   middle 
## 3.781853 4.513857 4.481195 
rbind(mean, sd)
##           high       low    middle
## mean 24.332189 23.311610 23.763508
## sd    3.781853  4.513857  4.481195

二つのカテゴリ(アルコールの摂取頻度と性別)で分類し、BMIの平均値を表示したい場合(二元配置分散分析)はこのようにします。

tapply(sample$BMI, list(sample$alcohol, sample$sex), mean, na.rm=TRUE)
##          Female     Male
## high   23.45400 24.68417
## low    22.84431 24.86355
## middle 23.01742 24.89114
2-3. sapply(x)
TG <- c(100.9, 98.0, 58.2, 62.9)
logTG <- sapply(TG, function(x) {return (log(x))})
logTG
## [1] 4.614130 4.584967 4.063885 4.141546
2-4. lapply(x)
TG <- c(100.9, 98.0, 58.2, 62.9)
logTG <- lapply(TG, function(x) {return (log(x))})
logTG
## [[1]]
## [1] 4.61413
## 
## [[2]]
## [1] 4.584967
## 
## [[3]]
## [1] 4.063885
## 
## [[4]]
## [1] 4.141546

基本的なグラフについて(デフォルトの図作成について)

基本的なグラフは新たなパッケージをインストールせずに描く方法が簡便かつ実行速度が速いため、非常に便利です。
3-1. ヒストグラム(hist(x))

デフォルトで描画する

hist(sample$age)

f:id:ryosukefujii0320:20160607001514p:plain

分割数(breaks)を20に指定し、x軸を30から70、y軸を0から250にする。軸のラベルをxlabとylabで指定し、図のタイトルをmainで記述している。

hist(sample$age, breaks=20, xlim=c(30, 70), ylim=c(0, 250), xlab="Age", ylab="Freq", main="The histgram of age")

f:id:ryosukefujii0320:20160607001521p:plain

3-2. 密度関数(plot(density(x)))
様々な設定をして、描画する
plot(density(sample$age), xlim=c(20, 80), xlab="Age", ylab="Prob", 
main="Density plot of age")

f:id:ryosukefujii0320:20160607001528p:plain

3-3. 散布図(plot(x, y))
plot(sample$age, sample$BMI, xlab="Age", ylab="BMI", main="The relationship between age and BMI")

f:id:ryosukefujii0320:20160607001534p:plain

3-3.(番外編)

自動的に散布図に直線を当てはめるscatter.smoothを使用する。
pchは点の種類、colは点の色、lwdは線の太さを指定するのに用いるが、今回は点の大きさを指定している。

scatter.smooth(sample$BMI ~ sample$age, xlab="Age", ylab="BMI", main="The relationship between age and BMI", pch=19, col="gray", lwd=3)

f:id:ryosukefujii0320:20160607001542p:plain

3-4-1. 棒グラフ(barplot(x))(シンプル)
mean <- tapply(sample$BMI, sample$alcohol, mean, na.rm=TRUE)
barplot(mean, ylim=c(0, 30), main="BMI by alcohol intake", 
        xlab="habitual alcohol intake", col=c("orange", "orangered", "darkred"))

f:id:ryosukefujii0320:20160607001553p:plain

#3-4-2. 棒グラフ(barplot(x))(縦積み)

counts <- table(sample$BMI01, sample$age)
barplot(counts, main="Obesity by Age", xlab="Age", col=c("darkblue","red"),
        legend = rownames(counts), ylim=c(0,150))

f:id:ryosukefujii0320:20160607001601p:plain

#3-4-3. 棒グラフ(横)

barplot(counts, main="Obesity by Age", xlab="Age", col=c("darkblue","red"),
        legend = rownames(counts), ylim=c(0,100), beside=TRUE)

f:id:ryosukefujii0320:20160607001608p:plain

#3-5. 箱ひげ図(boxplot(y ~ x))

boxplot(sample$BMI ~ sample$sex, xlab="Sex", ylab="BMI", main="The Boxplot of BMI by sex")

f:id:ryosukefujii0320:20160607001620p:plain

#4. ggplot2(特殊な描画パッケージ)を使用してのVisualization

library(ggplot2)

#4-1. ヒストグラム(geom_histogram)

#sampleを使用して、aesの中にx,yに何を表示するか記述する。
g <- ggplot(sample, aes(x = age)) 
#2歳づつに階級を分けてヒストグラムを描く
g <- g + geom_histogram(binwidth=2)
plot(g)

f:id:ryosukefujii0320:20160607001626p:plain

#4-2. 密度関数(geom_density)

#sampleを使用して、aesの中にx,yに何を表示するか記述する。
g <- ggplot(sample, aes(x = age))
#密度関数を描く
g <- g + geom_density()
plot(g)

f:id:ryosukefujii0320:20160607001632p:plain

#4-3. 散布図(geom_point)

g <- ggplot(sample, aes(x=age, y=BMI))
g <- g + geom_point()
plot(g)

f:id:ryosukefujii0320:20160607001639p:plain


#4-3-1. 線形回帰を当てはめる(geom_smoothで"lm"を指定している)
###95%信頼区間あり

#ここは4-3-2.と同じ
g <- ggplot(sample, aes(x = age, y = BMI))
#点のカタチ、サイズを指定し、NAは除く設定をしている
g <- g +  geom_point(shape = 20, size = 0.8, na.rm = TRUE)
#さらに線形回帰の直線を当てはめる。(glmなども可能)
g <- g + geom_smooth(method = "lm")
#ここからはラベルやタイトルの設定
g <- g + xlab("Age")
g <- g + ylab("BMI")
g <- g + ggtitle("The relationship between age and BMI (95%CI)") 
plot(g)

f:id:ryosukefujii0320:20160607001644p:plain

#4-4-1. 棒グラフ(geom_bar)

#分割したいグループの水準をgruopに格納する
group <- levels(sample$alcohol)
group
#tapplyで水準ごとのBMIの平均値を求める
BMImean <- tapply(sample$BMI, sample$alcohol, mean, na.rm=TRUE)
BMImean
#groupとBMImeanをまとめ、dtbarを作る
dtbar <- data.frame(level=group, mean=BMImean)
dtbar


###シンプル

#まとめたdtbarを用いる。xには水準、yには平均値を取る
g <- ggplot(dtbar, aes(x = level, y = mean))
g <- g + geom_bar(width = 0.8, stat = "identity")
g <- g + xlab("habitual alcohol intake")
g <- g + ylab("BMI")
g <- g + ggtitle("BMI by alcohol intake")
plot(g)

f:id:ryosukefujii0320:20160607001650p:plain

###棒の縁取りを黒、中をオレンジにする

g <- ggplot(dtbar, aes(x = level, y = mean))
g <- g + geom_bar(width = 0.8, stat = "identity", fill="orange", colour="black")
g <- g + xlab("habitual alcohol intake")
g <- g + ylab("BMI")
g <- g + ggtitle("BMI by alcohol intake")
plot(g)

f:id:ryosukefujii0320:20160607001655p:plain

###棒の色をそれぞれのレベルに対応する

g <- ggplot(dtbar, aes(x = level, y = mean, fill=level))
g <- g + geom_bar(width = 0.8, stat = "identity")
g <- g + xlab("habitual alcohol intake")
g <- g + ylab("BMI")
g <- g + ggtitle("BMI by alcohol intake")
plot(g)

f:id:ryosukefujii0320:20160607001701p:plain

#4-6. 箱ひげ図(geom_boxplot)
###シンプル

g <- ggplot(sample, aes(x = sex, y = BMI))
g <- g + geom_boxplot()
g <- g + ggtitle("The Boxplot of BMI by sex") 
g <- g + xlab("Sex") 
g <- g + ylab("BMI")
plot(g)

f:id:ryosukefujii0320:20160607001708p:plain

###色付き

g <- ggplot(sample, aes(x = sex, y = BMI, fill=sex))
g <- g + geom_boxplot()
g <- g + ggtitle("The Boxplot of BMI by sex") 
g <- g + xlab("Sex") 
g <- g + ylab("BMI")
plot(g)

f:id:ryosukefujii0320:20160607001714p:plain

##参考資料
###[1] apply系の関数の使い方(http://takenaka-akio.org/doc/r_auto/chapter_07_apply.html
###[2] 図の作成(http://stat.biopapyrus.net/
###[3] Winston Chang著『Rグラフィックスクックブック-ggplot2によるグラフ作成のレシピ-』日経印刷株式会社
###[4] The R Graph Gallery (http://www.r-graph-gallery.com/
###[5] 山本義郎・飯塚誠也・藤野友和著『統計データの視覚化-Rで学ぶデータサイエンス12-』共立出版

20160607
RF