Applied Statistics and Data Mining

大学院講義。

このファイル(昨年の課題)に収められたコードがすべて動くような環境を整えて下さい。

まずはRをダウンロードしてインストールします。CRANは、Rの開発組織が運営しているウェブサイトで、ここからのダウンロードは、とりあえず信用しても良いです。

• WindowsかMac OX Sなら、まずCRANにアクセスして下さい。
  1. Windowsなら、“Windows”の文字をクリックし、次に現れるページで“base”をクリックすると、“Download R X.Y.Z for Windows”という文字が見つかるので、それをクリックしてダウンロードし、実行する。
  2. Mac OS Xなら、“Mac OS X”の文字をクリックし、次に現れるページで “R-X.Y.Z.pkg (latest version)“という文字が見つかるので、それをクリックしてダウンロードし、実行する。
• Linuxなら、次の2行のどちらかが、インストールしてくれます。sudoが必要かも。

  yum install R R-*
  apt-get install r-base

次に、Rを起動して、必要なパッケージを追加します。今回は次の一つだけ。

• survivalパッケージをインストールするには、次の一行を実行します。

  install.packages(pkgs=c("survival"), repos='http://cran.md.tsukuba.ac.jp/')

この後、このファイルにあるコードをすべて実行してみてください。問題無く実行でき、赤い文字などでエラー(Error)やワーニング(Warning)が表示されなければ、成功です。

もしワーニングやエラーが表示されたら、山本まで、表示されたエラーメッセージなど、Rの出力画面の情報と共に、連絡を下さい。(ただし、メールアドレスの_at_は”@“に取り替えて、送信してください。)

この大学には、プライベートIPアドレスを持つPCが多くあります。それらには、NATルータやファイアウォールを介してインターネットに接続されているものと、NATルータなどは介さずにProxyサーバと呼ばれるサーバを経由してインターネットには接続せずに間接的に利用しているものがあります。以下は、この後者で発生する不具合への対応です。

プライベートIPアドレスを持ち、NATルータを介していないコンピュータからは、インターネットに直接アクセスすることはできず、PROXYを通す必要があります。 Windowsの場合は、Rを起動する際に

C:\Program Files\R\R-X.Y.Z\bin\Rgui.exe  (Xはバージョン、Yはメジャーリビジョン、Zはマイナーリビジョン)

となっているスタートメニュー内のショートカットの代わりに

"C:\Program Files\R\R-X.Y.Z\bin\Rgui.exe" --internet2

というショートカットを新たに作ると、

〔コントロールパネル〕→〔ネットワークとインターネット接続〕→〔インターネットオプション〕→〔接続〕→〔LANの設定〕→〔プロキシサーバー〕

で設定したプロキシを使用できます。

ダウンロードしたrats.datは、改行コードがUNIX用になっているため、一度、ブラウザで開いた上で、すべてをコピーし、メモ帳などに貼り付けてから、rats.datとして保存しないと、Rで読み込む際にエラーとなることがあります。

LaTeXで書いたテキストをdvipdfmxでPDFに変換すると、なぜか「~」の記号のフォントが、置き換えられてしまいます。見た目はほとんど変わりませんが、半角の「~」では無くなっています。改行コードが余計につくこともあります。 テキストのコードをそのままRに貼り付けると、プログラムとして認識できない文字が入力されたことになりますので、「~」は半角文字で打ち直すようにして下さい。 問題なく動作するコードを下に記します。

dimnames(rats)[[2]] <- c("id","start","stop","status","enum","trt")
rats$rtrunc <- rep(NA, nrow(rats))
rats$tstatus <- ifelse(rats$start == 0, 1, 2)
rats[1:12, c("id","start","stop","status","rtrunc","tstatus","enum","trt")]

NPfit <- coxph(Surv(start,stop,status)~1, data=rats, subset=(trt==0),method="breslow")
KM <- survfit(NPfit, type="aalen")
NA.MF <- data.frame(time=c(0,KM$time), na=-log(c(1,KM$surv)))
plot(NA.MF, type="l")

NPfit <- coxph(Surv(start,stop,status)~1, data=rats, subset=(trt==1),method="breslow")
KM <- survfit(NPfit, type="aalen")
NA.MF <- data.frame(time=c(0,KM$time), na=-log(c(1,KM$surv)))
lines(NA.MF, type="l")

NPfit <- coxph(Surv(start, stop, status)~trt, data=rats, method="breslow")
summary(NPfit)
KM <- survfit(NPfit, type="aalen")
NA.MF <- data.frame(time=c(0,KM$time), na=-log(c(1,KM$surv)))
plot(NA.MF)

NPfit <- coxph(Surv(start, stop, status)~trt+cluster(id), data=rats, method="breslow")
summary(NPfit)
KM <- survfit(NPfit, type="aalen")
NA.MF <- data.frame(time=c(0,KM$time), na=-log(c(1,KM$surv)))

library(survival)

2011.06.16に、今日より少し詳しい説明をしますので、コードを実行した出力だけは手元に持っていてもらえると助かります。

今回は Bladder Cancer Recurrences dataと呼ばれる、個人の再来院データを用いる。 観測されているのは、当初時点の、腫瘍の数(number)、最も大きい腫瘍の大きさ(size)、治療の種類(rx)、 それから履歴として、再来院(event)、及び間隔(start, stop)、である。

オリジナルデータには、再来院時の腫瘍の数や、大きさなど、時間依存の変数も含まれているが、今回は上のデータのみから、 病気の進行ではなく、再来院の頻度についての分析を行う。

Rで下記のコードを実行し、データ(bladder2)を少し変換したデータ(bladder.data)を作成する。

library(survival)
data(bladder)
bladder.data <- bladder2
bladder.data$rx[bladder.data$rx==1] <- 0
bladder.data$rx[bladder.data$rx==2] <- 1

最後に使うので、MASSパッケージも読み込んでおく。

library(MASS)

まず治療群(Thiotepaを用いている患者群)の計数過程の再来院確率関数を推定し、グラフに描く。 そしてコントロール群(プラセボを用いている患者群)の計数過程の生存関数を推定し、先のグラフに追記する。 ついでに両側95%信頼区間も描いた。 ここではKaplan-Meier推定量という方法を用いているが、詳細は省略する。

par(mfrow=c(2,1))
NPfit.1 <- coxph(Surv(start,stop,event)~1, data=bladder.data, subset=(rx==1),method="breslow")
KM.1 <- survfit(NPfit.1, conf.int=.95, type="aalen")
plot(KM.1)
NPfit.2 <- coxph(Surv(start,stop,event)~rx, data=bladder.data, subset=(rx==0),method="breslow")
KM.2 <- survfit(NPfit.2, conf.int=.95, type="aalen")
lines(KM.2, lty=3)
plot(KM.2)
lines(KM.1, lty=3)

95%信頼区間を重ね書きできないので、Treatment群の生存関数と信頼区間を描いてPlacebo群の生存関数を上書きしたものと、 Placebo群の生存関数と信頼区間を描いてTreatment群の生存関数を上書きしたもの、の2枚を用意した。

次に、上記と同等だが、それぞれの累積強度関数を推定してみる。 ここではNelson-Aalen推定量という推定方法を用いているが、詳細は省略する。

par(mfrow=c(1,1))
NPfit <- coxph(Surv(start,stop,event)~1, data=bladder.data, subset=(rx==1),method="breslow")
KM <- survfit(NPfit, conf.int=.95, type="aalen")
NA.MF <- data.frame(time=c(0,KM$time), na=-log(c(1,KM$surv)))
plot(NA.MF, type="l", lty=1,ylim=c(0,2.4), xlim=c(0,50))

NPfit <- coxph(Surv(start,stop,event)~rx, data=bladder.data, subset=(rx==0),method="breslow")
KM <- survfit(NPfit, conf.int=.95, type="aalen")
NA.MF <- data.frame(time=c(0,KM$time), na=-log(c(1,KM$surv)))
lines(NA.MF, type="l", lty=2)

これら2枚のグラフから、何が読み取れるか？

上と同じ事を、比例ハザードモデルで推定する。

NPfit <- coxph(Surv(start, stop, event)~factor(rx), data=bladder.data, method="breslow")
summary(NPfit)
KM <- survfit(NPfit, type="aalen")
NA.MF <- data.frame(time=c(0,KM$time), na=-log(c(1,KM$surv)))
lines(NA.MF, type="l", lty=4)

出力は次のようになる。

Call:
coxph(formula = Surv(start, stop, event) ~ factor(rx), data = bladder.data, method = "breslow")

  n= 178 

                  coef       exp(coef) se(coef)          z Pr(>|z|)  
factor(rx)1 -0.3655    0.6939   0.1976 -1.85   0.0644 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

            exp(coef) exp(-coef) lower .95 upper .95
factor(rx)1    0.6939      1.441    0.4711     1.022

Rsquare= 0.02   (max possible= 0.994 )
Likelihood ratio test= 3.52  on 1 df,   p=0.0605
Wald test            = 3.42  on 1 df,   p=0.06438
Score (logrank) test = 3.46  on 1 df,   p=0.06293

それぞれ、直訳すると、次の通り。

Call:
coxph(formula = Surv(start, stop, event) ~ factor(rx), data = bladder.data, method = "breslow")

  n= 178 

            　　推定値  exp(推定値)   標準誤差           t検定のp値 (5%以下なら有意)  
factor(rx)1 -0.3655    0.6939     0.1976 -1.85   0.0644 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

                 exp(coef) exp(-coef) 下側95%   上側95%
factor(rx)1    0.6939      1.441    0.4711     1.022

重相関係数の平方= 0.02   (max possible= 0.994 )  (今回は無視)
尤度比検定統計量= 3.52  on 1 df,   p=0.0605   (有意になると、モデルを統計的に支持する)
ワルド検定統計量= 3.42  on 1 df,   p=0.06438   (有意になると、モデルを統計的に支持する)
ログランク検定統計量= 3.46  on 1 df,   p=0.06293   (有意になると、モデルを統計的に支持する)

以下のモデルから出発し、t検定のp値を眺めて、5%以上の変数を外したりしながら、各種検定統計量が有意になるようなモデルを探して下さい。

NPfit <- coxph(Surv(start, stop, event)~factor(rx)*number*size+cluster(id), 
               data=bladder.data, method="breslow")
summary(NPfit)
scatter.smooth(residuals(NPfit))
abline(h=0,lty=3,col=2)

表示されるグラフは、残差プロットです。 本来はこのグラフも眺めながら、分析を進めていくものですが、今回はモデルを得ながら、このグラフも眺めてみよ。 注目するのは、点のばらつき具合。

上のモデルはフルモデルと呼ばれ、すべての変数の効果、すべての交互作用の効果を含んでいるため、下記のように書くこともできる。

NPfit <- coxph(Surv(start, stop, event)~factor(rx)+number+size
               +factor(rx):number+factor(rx):size+number:size+factor(rx):number:size+cluster(id), 
               data=bladder.data, method="breslow")
summary(NPfit)
scatter.smooth(residuals(NPfit))
abline(h=0,lty=3,col=2)

こちらから出発して、”+変数”を削る方が作業は簡単でしょう。 なお推定結果にcluster(id)の係数は現れません。 “+cluster(id)“は、各変数の標準誤差を、モデルが正しくはないかもしれない、ということを考慮しながら推定する、サンドイッチ推定量という方法を用いて推定するために必要な項なので、 これだけは削らないこと。

また、上のモデル当てはめと同時に、先の強度関数の推定値に、ベースライン強度関数の推定値を重ねていくと、モデルを変えてもほとんど変化しないことが見て取れる。 これは、回帰係数とベースライン強度関数を別々に推定できている、Cox比例ハザードモデルの特徴でもある。

AICというモデル選択基準は、これを最小にすると、データの予測精度の意味でのモデルの当てはまりが一番良い、という解釈を持つ。 上では、手動でフルモデルから変数を削ってもらったが、このAIC基準を用いて、自動的にモデル選択をさせることもできる。

NPfit <- coxph(Surv(start, stop, event)~factor(rx)*number*size+cluster(id), data=bladder.data, method="breslow")
stepAIC(NPfit)

上で得た、AICで最適なモデルに+cluster(id)を加え直し、問2で得たモデルと、ベースラインハザード関数、回帰係数、残差プロットなど、比較してみよ。

このデータの再来院に関して、以上の分析から、何が言えるか、まとめよ。

0.0 0 0
4.8 0 16012
6.0 0 16012
14.3 7 32027
22.8 7 48042
32.1 7 58854
41.4 7 69669
51.2 11 80483
60.6 12 91295
70.0 13 102110
79.9 15 112925
91.3 20 120367
97.0 21 127812
107.7 22 135257
119.1 28 142702
127.6 40 150147
135.1 44 152806
142.8 46 155464
148.9 48 158123
156.6 52 167081
163.9 52 167704
169.7 59 174626
170.1 59 174626
170.6 59 174626
174.7 63 181548
179.6 68 188473
185.5 71 194626
194.0 88 200782
200.3 93 206937
207.2 97 213093
211.9 98 219248
217.0 105 221355
223.5 113 223462
227.0 113 225568
234.1 122 227675
241.6 129 229784
250.7 141 233557
259.8 155 237330
268.3 166 241103
277.2 178 244879
285.5 186 247946
294.2 190 251016
295.7 190 251016
298.0 190 254086
305.2 195 257155
312.3 201 260225
318.2 209 260705
328.9 224 261188
334.8 231 261669
342.7 243 262889
350.5 252 263629
356.3 259 264367
360.6 271 265107
365.7 277 265845
374.9 282 266585
386.5 290 267325
396.5 300 268607
408.0 310 269891
417.3 312 271175
424.9 321 272457
434.2 326 273741
442.7 339 275025
451.4 346 276556
456.1 347 278087
460.8 351 279618
466.0 356 281149
472.3 359 283592
476.4 362 286036
480.9 367 288480
486.8 374 290923
495.8 376 293367
505.7 380 295811
516.0 392 298254
526.2 399 300698
527.3 401 300698
535.8 405 303142
546.3 415 304063
556.1 425 305009
568.1 440 305956
577.2 457 306902
578.3 457 306902
587.2 467 307849
595.5 473 308795
605.6 480 309742
613.9 491 310688
621.6 496 311635
623.4 496 311635
636.3 502 311750
649.7 517 311866
663.9 527 312467
675.1 540 313069
677.4 543 313069
677.9 544 313069
688.4 553 313671
698.1 561 314273
710.5 573 314783
720.9 581 315294
731.6 584 315805
732.7 585 315805
733.6 585 315805
746.7 586 316316
761.0 598 316827
776.5 612 318476
793.5 621 320125
807.2 636 321774
811.8 639 321774
812.5 639 321774
829.0 648 323423
844.4 658 325072
860.5 666 326179
876.7 674 327286
892.0 679 328393
895.5 686 328393
910.8 690 329500
925.1 701 330608
938.3 710 330435
952.0 720 330263
965.0 729 330091
967.7 729 330091
968.6 731 330091
981.3 740 329919
1013.9 759 330036
1030.1 776 330326
1044.0 781 330616
1047.0 782 330616
1059.7 783 330906
1072.6 787 331196
1085.7 793 331486
1098.4 796 331577
1112.4 797 331669
1113.5 798 331669
1114.1 798 331669
1128.0 802 331760
1139.1 805 331852
1163.2 823 332167
1174.3 827 332391
1184.6 832 332615
1198.3 834 332839
1210.3 836 333053
1221.1 839 333267
1230.5 842 333481
1231.6 842 333481
1240.9 844 333695
1249.5 845 333909
1262.2 849 335920
1271.3 851 337932
1279.8 854 339943
1287.4 855 341955
1295.1 859 341967
1304.8 860 341979
1305.8 865 342073
1313.3 867 342168
1314.4 867 342168
1320.0 867 342262
1325.3 867 342357
1330.6 870 342357
1334.2 870 342358
1336.7 870 342358