Rに出会ってはや二年。まっとうないちリウマチ科医としてすくすくと成長していた(?)はずの私が、臨床研究をきっかけに出会ったRから始まり、生物統計、疫学とはまって、数年前は考えたこともなかった公衆衛生修士の勉強(など)のために留学中という現状。それほどまでに人生を狂わせるRとはいったい何なのか。。
という話はさておいて、R Advent Calendar 2012( http://atnd.org/events/31973 )の一環として、臨床研究に役立ちそうなRのパッケージを緩めに厳選して紹介したい(と思ったが一個で力尽きた、、)。
ちなみに趣味でやっているR勉強会の資料はこちら
もはや向かう所敵なしのRの敵と言えば、メモリ上に乗り切らない大型のデータ(ff, bigmemoryなどのプロジェクトが進行中)とスタンダードなGUIが無くて統計解析環境としては敷居が高いことである。
臨床研究においてはTable 1というとExposed vs Unexposedの固定した二群間で、年齢だとか性別だとか、いろいろな変数をを比較するのが標準的。ただこれをRでやろうとすると割と面倒。そこで便利なのがDeducerというGUI ( http://www.slideshare.net/kaz_yos/introduction-to-deducer ) だが、ここではあえてDeducerをコマンドラインから使ってみる。
ちなみに、GUIでより高度な解析(生存分析やメタアナリシスなど)を行いたい場合は自治医大さいたま医療センターの神田先生が神懸かり的な努力で作成しているEZRがおすすめである。 http://www.jichi.ac.jp/saitama-sct/SaitamaHP.files/statmed.html 解説の本もある。英文の紹介はBone Marrow Transplantation誌に掲載されている http://www.nature.com/bmt/journal/vaop/ncurrent/pdf/bmt2012244a.pdf
Deducerをコマンドラインから使う場合は先に設定ファイルの.Rprofileに
Sys.setenv(NOAWT = 1)
という謎の呪文(Javaの何からしい)を書いておくのが望ましい。あるいは下記のごとくDeduderの読み込み前に実行する。
Sys.setenv(NOAWT = 1)
library(Deducer)
MASSパッケージのPima.tr(ピマインディアンの糖尿病のデータ)データを使ってみる。これはcross sectional studyのようなのでExposureとかOutcomeとかは言いにくいが、糖尿病の有無で二群に分けて、血圧、BMI、年齢をt-testで比較してみる。前半の文法は下記を意味する。
d(比較する変数1,変数2,変数3) ~ d(群分け変数)
library(MASS)
data(Pima.tr)
head(Pima.tr)
npreg glu bp skin bmi ped age type
1 5 86 68 28 30.2 0.364 24 No
2 7 195 70 33 25.1 0.163 55 Yes
3 5 77 82 41 35.8 0.156 35 No
4 0 165 76 43 47.9 0.259 26 No
5 0 107 60 25 26.4 0.133 23 No
6 5 97 76 27 35.6 0.378 52 Yes
## F-test 分散の違いの検定
two.sample.test(d(bp, bmi, age) ~ d(type), data = Pima.tr, test = var.test)
F test to compare two variances
ratio of variances 95% CI Lower 95% CI Upper F (num df,denom df) p-value
bp 0.91548 0.59341 1.3724 0.91548 (131,67) 0.660331
bmi 1.75945 1.14047 2.6376 1.75945 (131,67) 0.011153
age 0.69107 0.44795 1.0360 0.69107 (131,67) 0.073706
HA: two.sided
H0: ratio of variances = 1
## 等分散を仮定したt-test
two.sample.test(d(bp, bmi, age) ~ d(type), data = Pima.tr, test = t.test, var.equal = TRUE)
Two Sample t-test
mean of No mean of Yes Difference 95% CI Lower 95% CI Upper t df p-value
bp 69.545 74.588 -5.0428 -8.3559 -1.7296 -3.0015 198 0.003032359186
bmi 31.074 34.709 -3.6346 -5.3705 -1.8987 -4.1290 198 0.000053681285
age 29.235 37.691 -8.4563 -11.4706 -5.4420 -5.5323 198 0.000000099265
HA: two.sided
H0: difference in means = 0
## 等分散を仮定しないt-test
two.sample.test(d(bp, bmi, age) ~ d(type), data = Pima.tr, test = t.test, var.equal = FALSE)
Welch Two Sample t-test
mean of No mean of Yes Difference 95% CI Lower 95% CI Upper t df p-value
bp 69.545 74.588 -5.0428 -8.4141 -1.6715 -2.9592 130.28 0.0036648837
bmi 31.074 34.709 -3.6346 -5.2246 -2.0445 -4.5120 171.46 0.0000118817
age 29.235 37.691 -8.4563 -11.6674 -5.2453 -5.2162 115.70 0.0000008106
HA: two.sided
H0: difference in means = 0
Deducerを使うと変数ごとにコマンドを書かなくてすむ分だけ楽になる。BMIのみ等分散が棄却されているので、等分散を仮定しないt-testの結果を読むのが良いだろう。いずれの項目も糖尿病患者で有意に高いという結果であった。ちなみにこの解析はDeducerのGUIからやればなお簡単である( http://www.deducer.org/pmwiki/index.php?n=Main.TwoSampleTest )。
vcdパッケージのArthritis(だいぶ昔のRAの治験データらしい)データを使ってみる。実薬群とプラセボ群で分けて性別と治療効果(Marked, Some, Noneという今こんな論文書いたら1時間以内にrejectされそうな雑な分類。。)をFisher exact testで比較してみる。
library(vcd)
data(Arthritis)
head(Arthritis)
ID Treatment Sex Age Improved
1 57 Treated Male 27 Some
2 46 Treated Male 29 None
3 77 Treated Male 30 None
4 17 Treated Male 32 Marked
5 36 Treated Male 46 Marked
6 23 Treated Male 58 Marked
## クロス集計表を作成
res.tab <- contingency.tables(row.vars = d(Sex, Improved),
col.vars = d(Treatment), data = Arthritis)
## Fisher exact testの結果を追加
res.tab <- add.fishers.exact(res.tab)
## prop.cがここで必要な縦方向の割合になる。
print(res.tab, digits = 1, prop.r = F, prop.c = T, prop.t = F)
====================================================================================================================
======================================================================================
========== Table: Sex by Treatment ==========
| Treatment
Sex | Placebo | Treated | Row Total |
-----------------------|-----------|-----------|-----------|
Female Count | 32 | 27 | 59 |
Column % | 74.4% | 65.9% | |
-----------------------|-----------|-----------|-----------|
Male Count | 11 | 14 | 25 |
Column % | 25.6% | 34.1% | |
-----------------------|-----------|-----------|-----------|
Column Total | 43 | 41 | 84 |
Column % | 51.2% | 48.8% | |
| Exact
Test | Statistic DF p-value
Fishers Exact | 0.476
======================================================================================
========== Table: Improved by Treatment ==========
| Treatment
Improved | Placebo | Treated | Row Total |
-----------------------|-----------|-----------|-----------|
None Count | 29 | 13 | 42 |
Column % | 67.4% | 31.7% | |
-----------------------|-----------|-----------|-----------|
Some Count | 7 | 7 | 14 |
Column % | 16.3% | 17.1% | |
-----------------------|-----------|-----------|-----------|
Marked Count | 7 | 21 | 28 |
Column % | 16.3% | 51.2% | |
-----------------------|-----------|-----------|-----------|
Column Total | 43 | 41 | 84 |
Column % | 51.2% | 48.8% | |
| Exact
Test | Statistic DF p-value
Fishers Exact | 0.001
====================================================================================================================
無効がプラセボ群で67.4%/実薬群で31.7%、著効がプラセボ群で16.3%/実薬群で51.2%、p-value 0.001、ということで治療効果は有意に差があったようだ。
これまたDeducer GUI ( http://www.deducer.org/pmwiki/index.php?n=Main.ContingencyTables )から行った方が簡単だが、コマンドで残しておくと再現性がある(繰り返しができる)点は優れているだろう。
ということで、格調高いプログラミングの話( http://d.hatena.ne.jp/isobe1978/20121210/1355065348 )の直後に、Rをただの統計ソフトとして使おうというゆるい話を書かせていただきました。個人的には、Rをプログラミング言語とは思っていない(あるいは気づいていない?)“統計ソフトRのユーザ"がわりとたくさんいるのが、R界隈の面白さではないかと思います。
本当はこのあとに疫学の基本2x2表の話を書きたかったのですが、途中まで書いて内容がおかしいことに気づいたので、いったんなし。訂正できたら、後日公開します。
@kaz_yos