文書の過去の版を表示しています。
データマイニング
大学院講義。
演習の準備 2011.05.26
このファイル(昨年の課題)に収められたコードがすべて動くような環境を整えて下さい。
まずはRをダウンロードしてインストールします。CRANは、Rの開発組織が運営しているウェブサイトで、ここからのダウンロードは、とりあえず信用しても良いです。
- WindowsかMac OX Sなら、まずCRANにアクセスして下さい。
- Windowsなら、“Windows”の文字をクリックし、次に現れるページで“base”をクリックすると、“Download R X.Y.Z for Windows”という文字が見つかるので、それをクリックしてダウンロードし、実行する。
- 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で読み込む際にエラーとなることがあります。
PDFが原因
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)
演習1
日程
出題 | 2011.06.09 |
〆切 | 2011.06.23 |
2011.06.16に、今日より少し詳しい説明をしますので、コードを実行した出力だけは手元に持っていてもらえると助かります。
データの説明
今回は Bladder Cancer Recurrences dataと呼ばれる、個人の再来院データを用いる。 観測されているのは、当初時点の、腫瘍の数(number)、最も大きい腫瘍の大きさ(size)、治療の種類(rx)、 それから履歴として、再来院(event)、及び間隔(start, stop)、である。
変数 | 説明 | データ |
number | 当初の腫瘍の数 | 整数 |
size | 当初の最も大きい腫瘍の大きさ | 整数 |
rx | 治療の種類 | 整数 (0=プラセボ、1=pyrodixine or thiotepa) |
event | 事象の種類 | 整数 (0=打ち切り、1=来院) |
start | 区間の開始時点 | 整数(単位は月) |
stop | 区間の終了時点 | 整数(単位は月) |
id | 患者番号 | 整数 |
オリジナルデータには、再来院時の腫瘍の数や、大きさなど、時間依存の変数も含まれているが、今回は上のデータのみから、 病気の進行ではなく、再来院の頻度についての分析を行う。
準備
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)
問1
これら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 (有意になると、モデルを統計的に支持する)
問2
以下のモデルから出発し、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)
問3
上で得た、AICで最適なモデルに+cluster(id)を加え直し、問2で得たモデルと、ベースラインハザード関数、回帰係数、残差プロットなど、比較してみよ。
まとめ
このデータの再来院に関して、以上の分析から、何が言えるか、まとめよ。
演習
生データ
t Nt Ct 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