KU-R Session 2
Kazuki Yoshida
2015-03-23
本日の内容
解析に向けてのデータの整理
- データフレーム、ベクター
- カテゴリ変数
- 二値化
- サブセット
基本統計量
- 連続変数のサマリー mean(), median(), quantile(), range(), min(), max()
- カテゴリー変数のサマリー table(), xtabs(), formulaの説明
- summary()
- tableoneパッケージ
- by関数、doByパッケージ、dplyrパッケージによるグループ集計
検定
- 連続変数: t.test(), wilcox.test(), oneway.test(), kruskal.test()
- カテゴリカル変数: chisq.test(), fisher.test()
- 生存時間変数: survfit()とカプランマイヤー曲線, survdiff()
グラフ作成
- 利用可能なグラフシステムの種類
- baseグラフィックによる作図
- ggplot2パッケージによる作図
- 3Dグラフ
- グラフの出力
解析に向けてのデータの整理(データフレーム、ベクター、カテゴリ変数、二値化、サブセット)
データ形式
他の統計解析環境と比べた場合のユニークな点として、データ形式が非常に豊富ということがあげられます。通常のいゆるデータセット(SASやStataのそれと似たもの)はdata frameと呼ばれる形式です。これはデータの最小単位ではなく、data frameの列を形成しているvectorが最小単位になります。
## データ読み込み
nhefs <- XLConnect::readWorksheetFromFile("./nhefs_book.xls", sheet = 1)
## nhefsデータセットのデータ形式を確認
class(nhefs)## [1] "data.frame"## 一部を閲覧
head(nhefs)## seqn qsmk death yrdth sbp dbp sex age race income marital school
## 1 233 0 0 NA 175 96 0 42 1 19 2 7
## 2 235 0 0 NA 123 80 0 36 0 18 2 9
## 3 244 0 0 NA 115 75 1 56 1 15 3 11
## 4 245 0 1 85 148 78 0 68 1 15 3 5
## 5 252 0 0 NA 118 77 0 40 0 18 2 11
## 6 257 0 0 NA 141 83 1 43 1 11 4 9
## ht wt71 wt82 wt82_71 birthplace smokeintensity
## 1 174.1875 79.04 68.94604 -10.093960 47 30
## 2 159.3750 58.63 61.23497 2.604970 42 20
## 3 168.5000 56.81 66.22449 9.414486 51 20
## 4 170.1875 59.42 64.41012 4.990117 37 3
## 5 181.8750 87.09 92.07925 4.989251 42 20
## 6 162.1875 99.00 103.41906 4.419060 34 10
## smkintensity82_71 smokeyrs asthma bronch tb hf hbp pepticulcer colitis
## 1 -10 29 0 0 0 0 1 1 0
## 2 -10 24 0 0 0 0 0 0 0
## 3 -14 26 0 0 0 0 0 0 0
## 4 4 53 0 0 0 0 1 0 0
## 5 0 19 0 0 0 0 0 0 0
## 6 10 21 0 0 0 0 0 0 0
## hepatitis chroniccough hayfever diabetes polio tumor nervousbreak
## 1 0 0 0 1 0 0 0
## 2 0 0 0 0 0 0 0
## 3 0 0 1 0 0 1 0
## 4 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0
## 6 0 0 0 0 0 0 0
## alcoholpy alcoholfreq alcoholtype alcoholhowmuch pica headache otherpain
## 1 1 1 3 7 0 1 0
## 2 1 0 1 4 0 1 0
## 3 1 3 4 NA 0 1 1
## 4 1 2 3 4 0 0 1
## 5 1 2 1 2 0 1 0
## 6 1 3 2 1 0 1 0
## weakheart allergies nerves lackpep hbpmed boweltrouble wtloss infection
## 1 0 0 0 0 1 0 0 0
## 2 0 0 0 0 0 0 0 1
## 3 0 0 1 0 0 0 0 0
## 4 1 0 0 0 0 0 0 0
## 5 0 0 0 0 0 1 0 0
## 6 0 0 0 0 0 0 0 0
## active exercise birthcontrol pregnancies cholesterol hightax82 price71
## 1 0 2 2 NA 197 0 2.183594
## 2 0 0 2 NA 301 0 2.346680
## 3 0 2 0 2 157 0 1.569580
## 4 1 2 2 NA 174 0 1.506592
## 5 1 1 2 NA 216 0 2.346680
## 6 1 1 0 1 212 1 2.209961
## price82 tax71 tax82 price71_82 tax71_82
## 1 1.739990 1.1022949 0.4619751 0.44378662 0.6403809
## 2 1.797363 1.3649902 0.5718994 0.54931641 0.7929688
## 3 1.513428 0.5512695 0.2309875 0.05619812 0.3202515
## 4 1.451904 0.5249023 0.2199707 0.05479431 0.3049927
## 5 1.797363 1.3649902 0.5718994 0.54931641 0.7929688
## 6 2.025879 1.1547852 0.7479248 0.18408203 0.4069824## qsmk変数のvectorをとりだす
## $はこのdata frameの中でqsmkという名前の付いた部分を取り出すという意味です
head(nhefs$qsmk, 100)## [1] 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 1 0 0 1 0 1 0 0 0 0 1 0
## [36] 0 1 0 0 1 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 1 0
## [71] 0 1 0 0 1 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1下記は、ほぼ同じ意味です。
## このdata frameの中でqsmkという名前の付いた部分を取り出す
nhefs["qsmk"]
## qsmkという列を取り出す
nhefs[, "qsmk"]Rに特徴的なものとしてはlistと呼ばれるデータ形式があります。これは任意のデータを束ねておける非常に強力なデータ型式です。Rの解析結果は、ほぼ必ずこのlistと呼ばれるデータ型式で返ってきます。
## 適当なデータをlistに格納
lst1 <- list(data1 = head(nhefs$qsmk), data2 = tail(nhefs[,1:7]), number = 3, name = "test list")
## listを表示
lst1## $data1
## [1] 0 0 0 0 0 0
##
## $data2
## seqn qsmk death yrdth sbp dbp sex
## 1741 25014 0 0 NA 115 66 0
## 1742 25016 0 0 NA 124 80 1
## 1743 25024 0 0 NA NA NA 1
## 1744 25032 0 0 NA 171 77 0
## 1745 25033 0 0 NA 115 90 0
## 1746 25061 1 0 NA 136 90 0
##
## $number
## [1] 3
##
## $name
## [1] "test list"## listの一部だけ取り出す
lst1[["name"]]## [1] "test list"実はdata frameもlistの一種です。同じ長さのvectorが束ねられて表になっているわけです。
## data frameの1:10行目と1:5列目を抜き出し
nhefs[1:10, 1:5]## seqn qsmk death yrdth sbp
## 1 233 0 0 NA 175
## 2 235 0 0 NA 123
## 3 244 0 0 NA 115
## 4 245 0 1 85 148
## 5 252 0 0 NA 118
## 6 257 0 0 NA 141
## 7 262 0 0 NA 132
## 8 266 0 0 NA 100
## 9 419 0 1 84 163
## 10 420 0 1 86 184## デフォルトの表示方法をすると、listとして表示される。
print.default(nhefs[1:10, 1:5])## $seqn
## [1] 233 235 244 245 252 257 262 266 419 420
##
## $qsmk
## [1] 0 0 0 0 0 0 0 0 0 0
##
## $death
## [1] 0 0 0 1 0 0 0 0 1 1
##
## $yrdth
## [1] NA NA NA 85 NA NA NA NA 84 86
##
## $sbp
## [1] 175 123 115 148 118 141 132 100 163 184
##
## attr(,"class")
## [1] "data.frame"データの変換
解析の準備に必要な典型的なデータの変換の例を示します。
## 年齢を二値化して新しい変数としてデータセットに作成
summary(nhefs$age)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 25.00 33.00 43.00 43.64 53.00 74.00## 65歳以上かどうかの0,1変数を作成する。
nhefs$elderly <- as.numeric(nhefs$age >= 65)
table(nhefs$elderly)##
## 0 1
## 1653 93自動的にカテゴリ変数として扱って欲しい場合はfactorというデータ形式を与えます。名前を付けることもできます。
## カテゴリ変数を数値データからfactorに変換
nhefs$diab_cat <- factor(nhefs$diabetes)
## 元の変数は数値としてサマリーされる
summary(nhefs$diabetes)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.9868 2.0000 2.0000## 新しい変数はカテゴリとしてサマリーされる。
summary(nhefs$diab_cat)## 0 1 2
## 877 15 854## カテゴリ変数のレベルに名前を付ける
nhefs$diab_cat <- factor(nhefs$diabetes, levels = c(0,1,2), labels = c("Never","Ever","Missing"))
summary(nhefs$diab_cat)## Never Ever Missing
## 877 15 854欠測値はNAとしてあらわされます。前述の変数はNAとなっていなかったので変換します。
## nhefs$diab_cat変数で"Missing"と一致する場所をとりだす。
head(nhefs$diab_cat[nhefs$diab_cat == "Missing"])## [1] Missing Missing Missing Missing Missing Missing
## Levels: Never Ever Missing## それらをNAで置き換える
nhefs$diab_cat[nhefs$diab_cat == "Missing"] <- NA
## 結果確認
summary(nhefs$diab_cat)## Never Ever Missing NA's
## 877 15 0 854## NAのチェックは等号では行えません。
sum(nhefs$diab_cat == NA)## [1] NA## is.na()を用います。
sum(is.na(nhefs$diab_cat))## [1] 854データセットの一部のみでデータセットを作成する場合はsubsetを利用します。
## 高齢者の行だけでデータセットを作成
nhefs_elderly <- subset(nhefs, elderly == 1)
## 行数を確認
nrow(nhefs_elderly)## [1] 93## 高齢者変数の長さ
length(nhefs_elderly$elderly)## [1] 93## 変数の内容を確認
nhefs_elderly$elderly## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [36] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [71] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1年齢を10歳ごとに区切ってカテゴリ変数を作成します。
nhefs$ageCat <- cut(nhefs$age,
breaks = c(-Inf,30,40,50,60,70,Inf), # カットオフの指定。一番外は無限を使うとよい
labels = c("20s","30s","40s","50s","60s","70s"), # カテゴリ名
include.lowest = TRUE, # 一番最初の値も含める
right = FALSE) # [30, 40) 範囲の左側を閉じる。[以上,未満)にする。
summary(nhefs$ageCat)## 20s 30s 40s 50s 60s 70s
## 268 441 452 393 168 24年齢ときちんと対応しているかチェックします。ソートにはdplyrパッケージのarrange()関数が扱いやすいです。
library(dplyr)##
## Attaching package: 'dplyr'
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionunique(arrange(nhefs[c("age","ageCat")], age))## age ageCat
## 1 25 20s
## 63 26 20s
## 121 27 20s
## 165 28 20s
## 225 29 20s
## 269 30 30s
## 314 31 30s
## 356 32 30s
## 411 33 30s
## 457 34 30s
## 516 35 30s
## 557 36 30s
## 592 37 30s
## 626 38 30s
## 659 39 30s
## 710 40 40s
## 751 41 40s
## 787 42 40s
## 830 43 40s
## 878 44 40s
## 912 45 40s
## 962 46 40s
## 1013 47 40s
## 1065 48 40s
## 1113 49 40s
## 1162 50 50s
## 1221 51 50s
## 1266 52 50s
## 1296 53 50s
## 1340 54 50s
## 1391 55 50s
## 1433 56 50s
## 1476 57 50s
## 1507 58 50s
## 1537 59 50s
## 1555 60 60s
## 1577 61 60s
## 1599 62 60s
## 1624 63 60s
## 1639 64 60s
## 1654 65 60s
## 1672 66 60s
## 1685 67 60s
## 1697 68 60s
## 1711 69 60s
## 1723 70 70s
## 1731 71 70s
## 1733 72 70s
## 1741 74 70s基本統計量
各種の統計量は個別の関数になっていますので、まずはそちらから解説します。まずは、前回同様のデータを読み込みます。
## エクセルファイルを読み込みます
nhefs <- XLConnect::readWorksheetFromFile("./nhefs_book.xls", sheet = 1)連続変数のサマリー mean(), median(), quantile(), range(), min(), max()
連続変数の各種サマリーについて解説します。
## 年齢変数の平均
mean(nhefs$age)## [1] 43.63688## 死亡年齢の平均: 欠測値があると欠測値が返ります。
mean(nhefs$yrdth)## [1] NA## 死亡年齢の平均 欠測値以外を平均するには欠測値削除を指定します。
mean(nhefs$yrdth, na.rm = TRUE)## [1] 87.64545## 中央値
median(nhefs$age)## [1] 43## パーセンタイル
quantile(nhefs$age)## 0% 25% 50% 75% 100%
## 25 33 43 53 74## パーセンタイル 任意にprobオプションで指定できます
quantile(nhefs$age, probs = c(0.025, 0.25, 0.5, 0.75, 0.975))## 2.5% 25% 50% 75% 97.5%
## 25 33 43 53 68## range関数は最小値、最大値を返します。
range(nhefs$age)## [1] 25 74## 下記のように二つの結果をvectorにまとめて返すこともできます。
c(min(nhefs$age), max(nhefs$age))## [1] 25 74## 本当のデータの範囲は差をとれば良いです
diff(range(nhefs$age))## [1] 49## 関数として定義してみます。任意の関数を自作することが出来ます。
true_range <- function(x, na.rm = FALSE) {
max(x, na.rm = na.rm) - min(x, na.rm = na.rm)
}
true_range(nhefs$age)## [1] 49カテゴリー変数のサマリー table(), xtabs(), formulaの説明
カテゴリー変数のサマリーに使える関数とくりかえし出てくるformula指定の方法について説明します。
## カテゴリー化されている収入の変数を見てみます。
table(nhefs$income)##
## 11 12 13 14 15 16 17 18 19 20 21 22
## 29 60 68 70 85 80 67 312 439 245 127 98## NAを含めるには除外を無効にします。
table(nhefs$income, exclude = NULL)##
## 11 12 13 14 15 16 17 18 19 20 21 22 <NA>
## 29 60 68 70 85 80 67 312 439 245 127 98 66## 割合を見るにはテーブルの結果を割合変換の関数に渡します。
## このように関数の結果を他の関数に渡すときは処理が内側から外側に進行します
prop.table(table(nhefs$income, exclude = NULL))##
## 11 12 13 14 15 16
## 0.01660939 0.03436426 0.03894616 0.04009164 0.04868270 0.04581901
## 17 18 19 20 21 22
## 0.03837342 0.17869416 0.25143184 0.14032073 0.07273769 0.05612829
## <NA>
## 0.03780069## 見にくいのでパーセント表示にして、round関数で四捨五入します。
## 一段階目の結果を保存して、そちらに追加処理をしてみます。
propTab <- prop.table(table(nhefs$income, exclude = NULL))
round(propTab * 100, digits = 2)##
## 11 12 13 14 15 16 17 18 19 20 21 22
## 1.66 3.44 3.89 4.01 4.87 4.58 3.84 17.87 25.14 14.03 7.27 5.61
## <NA>
## 3.78## クロス集計ではxtabs関数がわかりやすいです。
## formulaは 左辺 ~ 右辺 の書式で、左辺を右辺で説明するという指定です。
## この場合は左辺は変数ではなく症例数なので空欄にします。
## 症例数をincomeとschoolで説明します。行、列の順番で指定します。
## dataにdata frameを指定します。
xtabs1 <- xtabs(formula = ~ income + school, data = nhefs)
xtabs1## school
## income 0 1 2 3 4 5 6 7 8 9 10 11 12 13 15 16 17
## 11 1 2 1 0 2 0 3 1 2 2 6 1 7 0 0 0 1
## 12 2 3 3 1 3 0 5 4 12 7 8 3 7 1 0 0 1
## 13 1 0 0 0 1 3 9 8 7 3 9 7 16 2 0 2 0
## 14 1 0 1 4 1 0 2 1 13 1 11 6 19 0 2 0 1
## 15 1 0 0 2 4 3 2 9 11 5 9 10 23 2 1 0 1
## 16 1 0 0 3 0 4 1 2 11 8 9 8 25 4 0 1 1
## 17 2 0 0 1 0 0 3 3 11 10 8 4 14 2 1 2 4
## 18 2 0 2 0 3 3 7 14 28 20 37 20 116 13 5 9 5
## 19 1 1 1 2 1 2 4 10 26 15 33 25 221 28 8 24 15
## 20 0 0 0 0 0 0 0 0 14 8 11 12 103 20 6 24 22
## 21 0 0 0 0 0 0 0 0 2 1 7 4 52 8 8 21 11
## 22 0 0 0 0 0 0 0 0 4 1 1 4 33 8 3 14 18## マージンに合計をつけたい場合は下記の様にします。
addmargins(xtabs1)## school
## income 0 1 2 3 4 5 6 7 8 9 10 11 12
## 11 1 2 1 0 2 0 3 1 2 2 6 1 7
## 12 2 3 3 1 3 0 5 4 12 7 8 3 7
## 13 1 0 0 0 1 3 9 8 7 3 9 7 16
## 14 1 0 1 4 1 0 2 1 13 1 11 6 19
## 15 1 0 0 2 4 3 2 9 11 5 9 10 23
## 16 1 0 0 3 0 4 1 2 11 8 9 8 25
## 17 2 0 0 1 0 0 3 3 11 10 8 4 14
## 18 2 0 2 0 3 3 7 14 28 20 37 20 116
## 19 1 1 1 2 1 2 4 10 26 15 33 25 221
## 20 0 0 0 0 0 0 0 0 14 8 11 12 103
## 21 0 0 0 0 0 0 0 0 2 1 7 4 52
## 22 0 0 0 0 0 0 0 0 4 1 1 4 33
## Sum 12 6 8 13 15 15 36 52 141 81 149 104 636
## school
## income 13 15 16 17 Sum
## 11 0 0 0 1 29
## 12 1 0 0 1 60
## 13 2 0 2 0 68
## 14 0 2 0 1 63
## 15 2 1 0 1 83
## 16 4 0 1 1 78
## 17 2 1 2 4 65
## 18 13 5 9 5 284
## 19 28 8 24 15 417
## 20 20 6 24 22 220
## 21 8 8 21 11 114
## 22 8 3 14 18 86
## Sum 88 34 97 80 1567## 二次元の表になっているので一部を切り出すことも出来ます。
xtabs1[1:2, 0:5]## school
## income 0 1 2 3 4
## 11 1 2 1 0 2
## 12 2 3 3 1 3## さらに説明変数を加えれば3次元以上の表になります。
## 行、列、階層化変数(複数可)という形の指定になります。
xtabs2 <- xtabs(formula = ~ income + school + sex, data = nhefs)
xtabs2## , , sex = 0
##
## school
## income 0 1 2 3 4 5 6 7 8 9 10 11 12 13 15 16 17
## 11 1 1 1 0 1 0 2 0 1 1 1 0 4 0 0 0 1
## 12 1 3 2 0 3 0 1 2 4 2 1 2 3 0 0 0 0
## 13 1 0 0 0 1 1 4 5 6 0 2 3 4 0 0 0 0
## 14 0 0 1 3 1 0 2 1 8 1 2 1 8 0 1 0 1
## 15 0 0 0 2 0 3 1 6 6 0 2 3 8 1 0 0 0
## 16 1 0 0 2 0 4 0 2 7 5 5 4 9 2 0 1 0
## 17 1 0 0 1 0 0 2 2 7 5 4 0 9 1 1 0 2
## 18 1 0 2 0 1 2 6 10 13 8 16 8 47 5 3 6 3
## 19 1 1 0 2 0 2 3 6 16 7 18 11 108 16 5 11 13
## 20 0 0 0 0 0 0 0 0 10 4 6 3 41 9 2 12 17
## 21 0 0 0 0 0 0 0 0 2 1 7 3 25 3 4 11 8
## 22 0 0 0 0 0 0 0 0 4 1 1 1 13 2 2 7 15
##
## , , sex = 1
##
## school
## income 0 1 2 3 4 5 6 7 8 9 10 11 12 13 15 16 17
## 11 0 1 0 0 1 0 1 1 1 1 5 1 3 0 0 0 0
## 12 1 0 1 1 0 0 4 2 8 5 7 1 4 1 0 0 1
## 13 0 0 0 0 0 2 5 3 1 3 7 4 12 2 0 2 0
## 14 1 0 0 1 0 0 0 0 5 0 9 5 11 0 1 0 0
## 15 1 0 0 0 4 0 1 3 5 5 7 7 15 1 1 0 1
## 16 0 0 0 1 0 0 1 0 4 3 4 4 16 2 0 0 1
## 17 1 0 0 0 0 0 1 1 4 5 4 4 5 1 0 2 2
## 18 1 0 0 0 2 1 1 4 15 12 21 12 69 8 2 3 2
## 19 0 0 1 0 1 0 1 4 10 8 15 14 113 12 3 13 2
## 20 0 0 0 0 0 0 0 0 4 4 5 9 62 11 4 12 5
## 21 0 0 0 0 0 0 0 0 0 0 0 1 27 5 4 10 3
## 22 0 0 0 0 0 0 0 0 0 0 0 3 20 6 1 7 3## sex = 1の階層だけ取り出してみます。"1"として文字列として指定してください。
xtabs2[,,"1"]## school
## income 0 1 2 3 4 5 6 7 8 9 10 11 12 13 15 16 17
## 11 0 1 0 0 1 0 1 1 1 1 5 1 3 0 0 0 0
## 12 1 0 1 1 0 0 4 2 8 5 7 1 4 1 0 0 1
## 13 0 0 0 0 0 2 5 3 1 3 7 4 12 2 0 2 0
## 14 1 0 0 1 0 0 0 0 5 0 9 5 11 0 1 0 0
## 15 1 0 0 0 4 0 1 3 5 5 7 7 15 1 1 0 1
## 16 0 0 0 1 0 0 1 0 4 3 4 4 16 2 0 0 1
## 17 1 0 0 0 0 0 1 1 4 5 4 4 5 1 0 2 2
## 18 1 0 0 0 2 1 1 4 15 12 21 12 69 8 2 3 2
## 19 0 0 1 0 1 0 1 4 10 8 15 14 113 12 3 13 2
## 20 0 0 0 0 0 0 0 0 4 4 5 9 62 11 4 12 5
## 21 0 0 0 0 0 0 0 0 0 0 0 1 27 5 4 10 3
## 22 0 0 0 0 0 0 0 0 0 0 0 3 20 6 1 7 3summary()
## summary()関数は与えたオブジェクトによって適切なサマリーを返します。
summary(nhefs)## seqn qsmk death yrdth
## Min. : 233 Min. :0.0000 Min. :0.0000 Min. :83.00
## 1st Qu.:10606 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:85.00
## Median :20346 Median :0.0000 Median :0.0000 Median :88.00
## Mean :16560 Mean :0.2658 Mean :0.1919 Mean :87.65
## 3rd Qu.:22763 3rd Qu.:1.0000 3rd Qu.:0.0000 3rd Qu.:90.00
## Max. :25061 Max. :1.0000 Max. :1.0000 Max. :93.00
## NA's :1416
## sbp dbp sex age
## Min. : 87.0 Min. : 47.00 Min. :0.0000 Min. :25.00
## 1st Qu.:115.0 1st Qu.: 70.00 1st Qu.:0.0000 1st Qu.:33.00
## Median :126.0 Median : 77.00 Median :1.0000 Median :43.00
## Mean :128.6 Mean : 77.78 Mean :0.5052 Mean :43.64
## 3rd Qu.:139.0 3rd Qu.: 85.00 3rd Qu.:1.0000 3rd Qu.:53.00
## Max. :229.0 Max. :130.00 Max. :1.0000 Max. :74.00
## NA's :79 NA's :84
## race income marital school
## Min. :0.0000 Min. :11.00 Min. :2.000 Min. : 0.00
## 1st Qu.:0.0000 1st Qu.:17.00 1st Qu.:2.000 1st Qu.:10.00
## Median :0.0000 Median :19.00 Median :2.000 Median :12.00
## Mean :0.1346 Mean :18.02 Mean :2.507 Mean :11.14
## 3rd Qu.:0.0000 3rd Qu.:20.00 3rd Qu.:2.000 3rd Qu.:12.00
## Max. :1.0000 Max. :22.00 Max. :8.000 Max. :17.00
## NA's :66 NA's :117
## ht wt71 wt82 wt82_71
## Min. :142.9 Min. : 36.17 Min. : 35.38 Min. :-41.280
## 1st Qu.:161.9 1st Qu.: 59.65 1st Qu.: 61.69 1st Qu.: -1.252
## Median :168.6 Median : 69.57 Median : 72.57 Median : 2.723
## Mean :168.9 Mean : 71.11 Mean : 73.63 Mean : 2.748
## 3rd Qu.:175.5 3rd Qu.: 80.26 3rd Qu.: 83.46 3rd Qu.: 6.805
## Max. :198.1 Max. :180.53 Max. :136.53 Max. : 48.538
## NA's :67 NA's :67
## birthplace smokeintensity smkintensity82_71 smokeyrs
## Min. : 1.00 Min. : 1.00 Min. :-80.000 Min. : 1.00
## 1st Qu.:21.00 1st Qu.:10.00 1st Qu.:-10.000 1st Qu.:15.00
## Median :34.00 Median :20.00 Median : -1.000 Median :24.00
## Mean :31.43 Mean :20.58 Mean : -4.812 Mean :24.59
## 3rd Qu.:42.00 3rd Qu.:30.00 3rd Qu.: 1.000 3rd Qu.:33.00
## Max. :56.00 Max. :80.00 Max. : 50.000 Max. :64.00
## NA's :97
## asthma bronch tb hf
## Min. :0.00000 Min. :0.00000 Min. :0.00000 Min. :0.000000
## 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:0.000000
## Median :0.00000 Median :0.00000 Median :0.00000 Median :0.000000
## Mean :0.04811 Mean :0.08419 Mean :0.01432 Mean :0.005155
## 3rd Qu.:0.00000 3rd Qu.:0.00000 3rd Qu.:0.00000 3rd Qu.:0.000000
## Max. :1.00000 Max. :1.00000 Max. :1.00000 Max. :1.000000
##
## hbp pepticulcer colitis hepatitis
## Min. :0.000 Min. :0.0000 Min. :0.00000 Min. :0.00000
## 1st Qu.:0.000 1st Qu.:0.0000 1st Qu.:0.00000 1st Qu.:0.00000
## Median :1.000 Median :0.0000 Median :0.00000 Median :0.00000
## Mean :1.056 Mean :0.1014 Mean :0.03322 Mean :0.01718
## 3rd Qu.:2.000 3rd Qu.:0.0000 3rd Qu.:0.00000 3rd Qu.:0.00000
## Max. :2.000 Max. :1.0000 Max. :1.00000 Max. :1.00000
##
## chroniccough hayfever diabetes polio
## Min. :0.00000 Min. :0.00000 Min. :0.0000 Min. :0.00000
## 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.:0.00000
## Median :0.00000 Median :0.00000 Median :0.0000 Median :0.00000
## Mean :0.05212 Mean :0.08877 Mean :0.9868 Mean :0.01489
## 3rd Qu.:0.00000 3rd Qu.:0.00000 3rd Qu.:2.0000 3rd Qu.:0.00000
## Max. :1.00000 Max. :1.00000 Max. :2.0000 Max. :1.00000
##
## tumor nervousbreak alcoholpy alcoholfreq
## Min. :0.00000 Min. :0.00000 Min. :0.000 Min. :0.000
## 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:1.000 1st Qu.:1.000
## Median :0.00000 Median :0.00000 Median :1.000 Median :2.000
## Mean :0.02348 Mean :0.02749 Mean :0.882 Mean :1.897
## 3rd Qu.:0.00000 3rd Qu.:0.00000 3rd Qu.:1.000 3rd Qu.:3.000
## Max. :1.00000 Max. :1.00000 Max. :2.000 Max. :5.000
##
## alcoholtype alcoholhowmuch pica headache
## Min. :1.000 Min. : 1.000 Min. :0.0000 Min. :0.0000
## 1st Qu.:1.000 1st Qu.: 2.000 1st Qu.:0.0000 1st Qu.:0.0000
## Median :3.000 Median : 2.000 Median :0.0000 Median :1.0000
## Mean :2.462 Mean : 3.281 Mean :0.9822 Mean :0.6266
## 3rd Qu.:3.000 3rd Qu.: 4.000 3rd Qu.:2.0000 3rd Qu.:1.0000
## Max. :4.000 Max. :48.000 Max. :2.0000 Max. :1.0000
## NA's :431
## otherpain weakheart allergies nerves
## Min. :0.0000 Min. :0.00000 Min. :0.00000 Min. :0.0000
## 1st Qu.:0.0000 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:0.0000
## Median :0.0000 Median :0.00000 Median :0.00000 Median :0.0000
## Mean :0.2434 Mean :0.02176 Mean :0.06243 Mean :0.1397
## 3rd Qu.:0.0000 3rd Qu.:0.00000 3rd Qu.:0.00000 3rd Qu.:0.0000
## Max. :1.0000 Max. :1.00000 Max. :1.00000 Max. :1.0000
##
## lackpep hbpmed boweltrouble wtloss
## Min. :0.00000 Min. :0.000 Min. :0.000 Min. :0.00000
## 1st Qu.:0.00000 1st Qu.:0.000 1st Qu.:0.000 1st Qu.:0.00000
## Median :0.00000 Median :1.000 Median :1.000 Median :0.00000
## Mean :0.04983 Mean :1.011 Mean :1.041 Mean :0.02577
## 3rd Qu.:0.00000 3rd Qu.:2.000 3rd Qu.:2.000 3rd Qu.:0.00000
## Max. :1.00000 Max. :2.000 Max. :2.000 Max. :1.00000
##
## infection active exercise birthcontrol
## Min. :0.000 Min. :0.0000 Min. :0.000 Min. :0.000
## 1st Qu.:0.000 1st Qu.:0.0000 1st Qu.:1.000 1st Qu.:0.000
## Median :0.000 Median :1.0000 Median :1.000 Median :1.000
## Mean :0.146 Mean :0.6558 Mean :1.189 Mean :1.089
## 3rd Qu.:0.000 3rd Qu.:1.0000 3rd Qu.:2.000 3rd Qu.:2.000
## Max. :1.000 Max. :2.0000 Max. :2.000 Max. :2.000
##
## pregnancies cholesterol hightax82 price71
## Min. : 1.000 Min. : 78.0 Min. :0.0000 Min. :1.507
## 1st Qu.: 2.000 1st Qu.:189.0 1st Qu.:0.0000 1st Qu.:2.037
## Median : 3.000 Median :216.0 Median :0.0000 Median :2.168
## Mean : 3.662 Mean :219.6 Mean :0.1631 Mean :2.139
## 3rd Qu.: 5.000 3rd Qu.:245.0 3rd Qu.:0.0000 3rd Qu.:2.242
## Max. :15.000 Max. :416.0 Max. :1.0000 Max. :2.693
## NA's :977 NA's :20 NA's :97 NA's :97
## price82 tax71 tax82 price71_82
## Min. :1.452 Min. :0.5249 Min. :0.2200 Min. :-0.2027
## 1st Qu.:1.740 1st Qu.:0.9449 1st Qu.:0.4399 1st Qu.: 0.2010
## Median :1.815 Median :1.0498 Median :0.5060 Median : 0.3430
## Mean :1.807 Mean :1.0587 Mean :0.5062 Mean : 0.3324
## 3rd Qu.:1.868 3rd Qu.:1.1548 3rd Qu.:0.5719 3rd Qu.: 0.4438
## Max. :2.103 Max. :1.5225 Max. :0.7479 Max. : 0.6121
## NA's :97 NA's :97 NA's :97 NA's :97
## tax71_82 elderly diab_cat ageCat
## Min. :0.0360 Min. :0.00000 Never :877 20s:268
## 1st Qu.:0.4610 1st Qu.:0.00000 Ever : 15 30s:441
## Median :0.5440 Median :0.00000 Missing: 0 40s:452
## Mean :0.5525 Mean :0.05326 NA's :854 50s:393
## 3rd Qu.:0.6220 3rd Qu.:0.00000 60s:168
## Max. :0.8844 Max. :1.00000 70s: 24
## NA's :97summary(nhefs$age)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 25.00 33.00 43.00 43.64 53.00 74.00summary(nhefs$ageCat)## 20s 30s 40s 50s 60s 70s
## 268 441 452 393 168 24tableoneパッケージ
医療系の研究でよく見られる表1を作成するためのパッケージを書きました。
library(tableone)
## qsmkで群分けした表
tab1a <- CreateTableOne(vars = c("age","death","sbp","dbp","sex"), strata = "qsmk", data = nhefs)## NOTE: no factor/logical/character variables supplied, using CreateContTable()tab1a## Stratified by qsmk
## 0 1 p test
## n 1282 464
## age (mean (sd)) 42.68 (11.80) 46.27 (12.53) <0.001
## death (mean (sd)) 0.18 (0.38) 0.23 (0.42) 0.009
## sbp (mean (sd)) 127.67 (18.84) 131.27 (19.59) 0.001
## dbp (mean (sd)) 77.44 (10.60) 78.76 (10.77) 0.027
## sex (mean (sd)) 0.53 (0.50) 0.45 (0.50) 0.003## deathとsexを表作成時にカテゴリカル化します
tab1b <- CreateTableOne(vars = c("age","death","sbp","dbp","sex"), strata = "qsmk", data = nhefs,
factorVars = c("death","sex"))
tab1b## Stratified by qsmk
## 0 1 p test
## n 1282 464
## age (mean (sd)) 42.68 (11.80) 46.27 (12.53) <0.001
## death = 1 (%) 227 (17.7) 108 (23.3) 0.011
## sbp (mean (sd)) 127.67 (18.84) 131.27 (19.59) 0.001
## dbp (mean (sd)) 77.44 (10.60) 78.76 (10.77) 0.027
## sex = 1 (%) 675 (52.7) 207 (44.6) 0.004## すべての連続変数をnonnormalとして扱う
print(tab1b, nonnormal = TRUE)## Stratified by qsmk
## 0 1
## n 1282 464
## age (median [IQR]) 42.00 [32.00, 51.00] 46.50 [35.00, 56.00]
## death = 1 (%) 227 (17.7) 108 (23.3)
## sbp (median [IQR]) 125.00 [115.00, 138.00] 129.00 [117.00, 142.00]
## dbp (median [IQR]) 77.00 [70.00, 84.00] 78.00 [71.00, 86.00]
## sex = 1 (%) 675 (52.7) 207 (44.6)
## Stratified by qsmk
## p test
## n
## age (median [IQR]) <0.001 nonnorm
## death = 1 (%) 0.011
## sbp (median [IQR]) <0.001 nonnorm
## dbp (median [IQR]) 0.028 nonnorm
## sex = 1 (%) 0.004## ageのみをnonnormalとして扱う
print(tab1b, nonnormal = "age")## Stratified by qsmk
## 0 1 p
## n 1282 464
## age (median [IQR]) 42.00 [32.00, 51.00] 46.50 [35.00, 56.00] <0.001
## death = 1 (%) 227 (17.7) 108 (23.3) 0.011
## sbp (mean (sd)) 127.67 (18.84) 131.27 (19.59) 0.001
## dbp (mean (sd)) 77.44 (10.60) 78.76 (10.77) 0.027
## sex = 1 (%) 675 (52.7) 207 (44.6) 0.004
## Stratified by qsmk
## test
## n
## age (median [IQR]) nonnorm
## death = 1 (%)
## sbp (mean (sd))
## dbp (mean (sd))
## sex = 1 (%)by関数、doByパッケージ、dplyrパッケージによるグループ集計
グループごとの集計はtableoneが簡易的ですが、他にもいくつか方法があります。
## by関数で体重を年齢カテゴリごとにサマリーしてみます。
by(nhefs$wt71, nhefs$ageCat, summary)## nhefs$ageCat: 20s
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 40.82 57.01 67.02 69.98 80.66 169.20
## --------------------------------------------------------
## nhefs$ageCat: 30s
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 42.30 57.83 68.27 70.88 80.29 180.50
## --------------------------------------------------------
## nhefs$ageCat: 40s
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 40.82 60.75 70.53 71.94 81.59 134.80
## --------------------------------------------------------
## nhefs$ageCat: 50s
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 36.17 60.89 70.42 71.55 79.95 135.60
## --------------------------------------------------------
## nhefs$ageCat: 60s
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 39.58 60.44 70.70 70.50 78.70 121.00
## --------------------------------------------------------
## nhefs$ageCat: 70s
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 44.34 64.21 70.93 69.51 76.18 93.67## doByパッケージのsummaryByも便利かもしれません。
## formulaの左側にサマリーされる変数、右側にグループ分けの変数を書きます。
library(doBy)## Loading required package: survivalsummaryBy(wt71 + ht ~ ageCat, data = nhefs, FUN = mean)## ageCat wt71.mean ht.mean
## 1 20s 69.98157 169.8759
## 2 30s 70.87832 169.2659
## 3 40s 71.93912 169.7029
## 4 50s 71.55331 167.7697
## 5 60s 70.50173 167.4388
## 6 70s 69.51333 163.7839## 最近はdplyrパッケージが人気です。
## パイプ %>% と呼ばれる書式が使われます。
nhefs %>%
group_by(ageCat) %>%
summarize(weight_median = median(wt71),
weight_25perc = quantile(wt71, probs = 0.25),
weight_75perc = quantile(wt71, probs = 0.75),
height_median = median(ht))## Source: local data frame [6 x 5]
##
## ageCat weight_median weight_25perc weight_75perc height_median
## 1 20s 67.015 57.0125 80.6575 169.5312
## 2 30s 68.270 57.8300 80.2900 168.5938
## 3 40s 70.530 60.7525 81.5875 169.7344
## 4 50s 70.420 60.8900 79.9500 167.6875
## 5 60s 70.705 60.4400 78.7000 166.1875
## 6 70s 70.930 64.2075 76.1775 165.9375検定
連続変数: t.test(), wilcox.test(), oneway.test(), kruskal.test()
## 二群
## 正規分布
## デフォルトでは両群の分散の一致を仮定しない検定になります。
t.test(age ~ sex, data = nhefs)##
## Welch Two Sample t-test
##
## data: age by sex
## t = 1.7044, df = 1735.968, p-value = 0.08848
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.1487584 2.1225396
## sample estimates:
## mean in group 0 mean in group 1
## 44.13542 43.14853## 通常のt検定はこちら
t.test(age ~ sex, data = nhefs, var.equal = TRUE)##
## Two Sample t-test
##
## data: age by sex
## t = 1.7053, df = 1744, p-value = 0.08833
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.1481994 2.1219806
## sample estimates:
## mean in group 0 mean in group 1
## 44.13542 43.14853## 正規性のチェックはQQ plotが良いでしょうか。
car::qqPlot(nhefs$age)
## 非正規分布
## Wilcoxon test()
wilcox.test(age ~ sex, data = nhefs)##
## Wilcoxon rank sum test with continuity correction
##
## data: age by sex
## W = 397971.5, p-value = 0.1075
## alternative hypothesis: true location shift is not equal to 0## 多群
## ANOVAはoneway.test()で出来ます。
## これもデフォルトでは分散の一致を仮定しないものになります。
oneway.test(age ~ alcoholtype, data = nhefs)##
## One-way analysis of means (not assuming equal variances)
##
## data: age and alcoholtype
## F = 32.5121, num df = 3.000, denom df = 579.648, p-value < 2.2e-16oneway.test(age ~ alcoholtype, data = nhefs, var.equal = TRUE)##
## One-way analysis of means
##
## data: age and alcoholtype
## F = 31.3619, num df = 3, denom df = 1742, p-value < 2.2e-16## 非正規の場合はkruskal.testが使えます
kruskal.test(age ~ alcoholtype, data = nhefs)##
## Kruskal-Wallis rank sum test
##
## data: age by alcoholtype
## Kruskal-Wallis chi-squared = 95.596, df = 3, p-value < 2.2e-16ANOVAに関しては実際は回帰分析の関数を使って行う方が一般的です。平方和のテーブルなどが出ます。
lm1 <- lm(age ~ alcoholtype, data = nhefs)
anova(lm1)## Analysis of Variance Table
##
## Response: age
## Df Sum Sq Mean Sq F value Pr(>F)
## alcoholtype 1 11699 11699.3 83.734 < 2.2e-16 ***
## Residuals 1744 243671 139.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Rでは結果がオブジェクトとして返ってきます。名前を付けて保存して内容をチェックします。
resTTest <- t.test(age ~ sex, data = nhefs)
## 名前だけ評価すると適切な方法で表示されます
resTTest##
## Welch Two Sample t-test
##
## data: age by sex
## t = 1.7044, df = 1735.968, p-value = 0.08848
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.1487584 2.1225396
## sample estimates:
## mean in group 0 mean in group 1
## 44.13542 43.14853## 要素の名前をチェックします。
names(resTTest)## [1] "statistic" "parameter" "p.value" "conf.int" "estimate"
## [6] "null.value" "alternative" "method" "data.name"## p-valueだけ抜きだします。
resTTest$p.value## [1] 0.08848213カテゴリカル変数: chisq.test(), fisher.test()
カテゴリカル変数では前述の関数を使ってクロス集計表を作成して検定します。
## 表作成
tabSexQsmk <- xtabs( ~ sex + qsmk, data = nhefs)
## カイ二乗検定 デフォルトでは連続性補正をします。
chisq.test(tabSexQsmk)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tabSexQsmk
## X-squared = 8.4915, df = 1, p-value = 0.003568## 通常のカイ二乗検定
chisq.test(tabSexQsmk, correct = FALSE)##
## Pearson's Chi-squared test
##
## data: tabSexQsmk
## X-squared = 8.8102, df = 1, p-value = 0.002996## すでに表のデータの場合は下記のように入力しても良いでしょう。
## 4つの値を与えて、2列の行列を作成しています。標準では縦方向に数字が入ります。
## byrowで横向きに変更できます。
mat <- matrix(c(5,1,10,20), ncol = 2, byrow = TRUE)
resFisher <- fisher.test(mat)
resFisher##
## Fisher's Exact Test for Count Data
##
## data: mat
## p-value = 0.06281
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.8811891 494.4191744
## sample estimates:
## odds ratio
## 9.364462## 同じようにオブジェクトが返っています。
names(resFisher)## [1] "p.value" "conf.int" "estimate" "null.value" "alternative"
## [6] "method" "data.name"## 2x2表のオッズ比を抜き出してみます。
resFisher$estimate## odds ratio
## 9.364462生存時間変数: survfit()とカプランマイヤー曲線, survdiff()
デモデータが生存分析データではないので別のデータを使います。
## survivalパッケージから白血病のデータを読み込みます。
library(survival)
data(leukemia)
## チェック
head(leukemia)## time status x
## 1 9 1 Maintained
## 2 13 1 Maintained
## 3 13 0 Maintained
## 4 18 1 Maintained
## 5 23 1 Maintained
## 6 28 0 Maintained## カプランーマイヤーはこのように行います。
## 群分けしない場合はformulaの右には1とします (モデルに説明変数がなく切片のみという意味)
## 信頼区間のタイプはlog-log変換にしておいた方が狭くなります。
surv1 <- survfit(Surv(time, status) ~ 1, data = leukemia, conf.type = "log-log")
## 標準の結果表示では限られた情報のみ
surv1## Call: survfit(formula = Surv(time, status) ~ 1, data = leukemia, conf.type = "log-log")
##
## records n.max n.start events median 0.95LCL 0.95UCL
## 23 23 23 18 27 13 34## summaryするとイベント時間ごとの集計がみられます。
summary(surv1)## Call: survfit(formula = Surv(time, status) ~ 1, data = leukemia, conf.type = "log-log")
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 5 23 2 0.9130 0.0588 0.69495 0.978
## 8 21 2 0.8261 0.0790 0.60061 0.931
## 9 19 1 0.7826 0.0860 0.55421 0.903
## 12 18 1 0.7391 0.0916 0.50921 0.873
## 13 17 1 0.6957 0.0959 0.46564 0.842
## 18 14 1 0.6460 0.1011 0.41395 0.805
## 23 13 2 0.5466 0.1073 0.31925 0.726
## 27 11 1 0.4969 0.1084 0.27557 0.684
## 30 9 1 0.4417 0.1095 0.22738 0.637
## 31 8 1 0.3865 0.1089 0.18284 0.587
## 33 7 1 0.3313 0.1064 0.14183 0.535
## 34 6 1 0.2761 0.1020 0.10444 0.480
## 43 5 1 0.2208 0.0954 0.07100 0.422
## 45 4 1 0.1656 0.0860 0.04211 0.360
## 48 2 1 0.0828 0.0727 0.00696 0.287## 任意の時間を指定するにはtimesを使います。
summary(surv1, times = c(0, 20, 40, 60))## Call: survfit(formula = Surv(time, status) ~ 1, data = leukemia, conf.type = "log-log")
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0 23 0 1.0000 0.0000 1.00000 1.000
## 20 13 8 0.6460 0.1011 0.41395 0.805
## 40 5 7 0.2761 0.1020 0.10444 0.480
## 60 1 3 0.0828 0.0727 0.00696 0.287## KMはこの結果オブジェクトをplotすると描画されます。
plot(surv1)
## 信頼区間が不要な場合
plot(surv1, conf.int = FALSE)
## 群分けする場合はformulaを変えます。
surv2 <- survfit(Surv(time, status) ~ x, data = leukemia, conf.type = "log-log")
surv2## Call: survfit(formula = Surv(time, status) ~ x, data = leukemia, conf.type = "log-log")
##
## records n.max n.start events median 0.95LCL 0.95UCL
## x=Maintained 11 11 11 7 31 13 NA
## x=Nonmaintained 12 12 12 11 23 5 33summary(surv2)## Call: survfit(formula = Surv(time, status) ~ x, data = leukemia, conf.type = "log-log")
##
## x=Maintained
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 9 11 1 0.909 0.0867 0.5081 0.987
## 13 10 1 0.818 0.1163 0.4474 0.951
## 18 8 1 0.716 0.1397 0.3502 0.899
## 23 7 1 0.614 0.1526 0.2658 0.835
## 31 5 1 0.491 0.1642 0.1673 0.753
## 34 4 1 0.368 0.1627 0.0928 0.657
## 48 2 1 0.184 0.1535 0.0117 0.525
##
## x=Nonmaintained
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 5 12 2 0.8333 0.1076 0.48171 0.956
## 8 10 2 0.6667 0.1361 0.33702 0.860
## 12 8 1 0.5833 0.1423 0.27014 0.801
## 23 6 1 0.4861 0.1481 0.19188 0.730
## 27 5 1 0.3889 0.1470 0.12627 0.650
## 30 4 1 0.2917 0.1387 0.07240 0.561
## 33 3 1 0.1944 0.1219 0.03120 0.461
## 43 2 1 0.0972 0.0919 0.00575 0.349
## 45 1 1 0.0000 NaN NA NAsummary(surv2, times = c(0, 20, 40, 60))## Call: survfit(formula = Surv(time, status) ~ x, data = leukemia, conf.type = "log-log")
##
## x=Maintained
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0 11 0 1.000 0.000 1.0000 1.000
## 20 7 3 0.716 0.140 0.3502 0.899
## 40 3 3 0.368 0.163 0.0928 0.657
## 60 1 1 0.184 0.153 0.0117 0.525
##
## x=Nonmaintained
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0 12 0 1.000 0.000 1.0000 1.000
## 20 6 5 0.583 0.142 0.2701 0.801
## 40 2 4 0.194 0.122 0.0312 0.461plot(surv2)
## log-rank検定の書式はsurvfitと同じです。
survdiff(Surv(time, status) ~ x, data = leukemia)## Call:
## survdiff(formula = Surv(time, status) ~ x, data = leukemia)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## x=Maintained 11 7 10.69 1.27 3.4
## x=Nonmaintained 12 11 7.31 1.86 3.4
##
## Chisq= 3.4 on 1 degrees of freedom, p= 0.0653グラフ作成
利用可能なグラフシステムの種類
Rにはグラフ作成のためのシステムがいくつかあります。おおまかにいって、traditional (base)と呼ばれるシステムとgridと呼ばれるシステムがあります。gridの方は最近は、実質的にはggplot2というパッケージが主に使われます。
baseグラフィックによる作図
通常のplot関数はbaseの方になります。formulaを使ったグラフ指定ができます。plot関数はオブジェクトの種類に応じてさまざまなものをplotできます。統計モデルに使用するとモデル評価を表示しますがこれらはモデリングのところで扱います。先程すでに出たカプランマイヤー曲線の描画もbaseを使用しています。
## 連続変数 ~ カテゴリカル変数は自動的にboxplotになります。
plot(time ~ x, data = aml)
## 連続変数同士は散布図になります
plot(wt71 ~ age, nhefs)
## 軸のラベルは下記のようにおこないます。
plot(wt71 ~ age, nhefs,
xlab = "Age in 1971", ylab = "Weight in 1971", main = "Scatterplot", sub = "subtitle")
## 軸の範囲は下記のように指定します。
plot(wt71 ~ age, nhefs,
xlim = c(0,100), ylim = c(0,200))
ggplot2パッケージによる作図
Grammer of Graphics (gg)というグラフをどのように言語として表現すべきかという考えかたがあります。これをRパッケージとして具現化したのがggplot2パッケージです。かなりくせがあるように感じますが、実に理路整然とした形でグラフを言語表現します。簡略記法もあるのですが、逆にggの概念を隠してしまうので、あえて完全記法で記載します。
## ggplot2パッケージを読み込みます。
library(ggplot2)
## デフォルトのデータセットとどの軸(scale)にどの変数を割り振るかを指定します。
## +サインで表現を繋いでいきます。
## レイヤーとして形態(geometry)が散布図(point)のものを加えます。
## baseと違ってグラフをオブジェクトとして保存できます。
plot1 <- ggplot(data = nhefs, mapping = aes(x = age, y = wt71)) +
layer(geom = "point")
## 名前だけ評価すれば表示されます。
plot1
## themeという要素をつかってグラフの仕上げを調節します。
## theme_bwというプリセットで背景を白にできます。
plot1 + theme_bw()
## パネルを変数で割るにはfacetという要素をつかいます。
## 1変数でわるならfacet_wrap
plot1 + facet_wrap( ~ sex)
## 縦横に展開するならfacet_gridを使います
## ラベリングはlabsを使います。
plot1 + facet_grid(sex ~ qsmk) + labs(title = "sex by qsmk", x = "Age in 1971", y = "Weight in 1971")
## 男女で色分けしてみます。カテゴリカルで与える必要があるのでfactorを使用します。
plot2 <- ggplot(data = nhefs, mapping = aes(x = age, y = wt71, color = factor(sex))) +
layer(geom = "point") + labs(color = "sex")
plot2
## パネルを分けた方が見易いですね。
plot2 + facet_wrap( ~ sex)
3Dグラフ
いくつか3Dグラフを描画する方法もあります。
library(scatterplot3d)
## with関数を使うと表記を簡略化できます
with(nhefs,
scatterplot3d(x = age, y = ht, z = wt71))
## あるいは
scatterplot3d(x = nhefs$age, y = nhefs$ht, z = nhefs$wt71)
実用性は微妙ですがインタラクティブな3Dグラフを作成することもできます。
##
library(rgl)
## 男女で青赤に塗り分けています
with(nhefs,
plot3d(x = age, y = ht, z = wt71,
type = "s", size = 1,
col = c("blue","red")[sex + 1]))グラフの出力
グラフの出力はファイル形式に応じた関数を用いて行います。
## PDFの場合は複数ページ対応しています。
pdf(file = "plots.pdf", width = 7, height = 7, family = "sans")
plot(time ~ x, data = aml)
plot(wt71 ~ age, nhefs)
plot1
## ファイルを閉じて完了というシグナルを送る必要があります。
dev.off()
pdf(file = "plot2.pdf", width = 7, height = 7, family = "sans")
## 3Dグラフの見易い角度を探すためにforループで複数角度を変えながら作成してみます
for (i in seq(0, 360, 10)) {
scatterplot3d(x = nhefs$age, y = nhefs$ht, z = nhefs$wt71, angle = i, main = i)
}
dev.off()
## pngの場合(jpegは関数名以外同じ)は一枚づつです。
png(file = "plot1.png", width = 700, height = 700, family = "sans")
plot1
dev.off()
png(file = "plot2.png", width = 700, height = 700, family = "sans")
plot2
dev.off()