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(3,1))
NPfit <- coxph(Surv(start,stop,event)~1, data=bladder.data, subset=(rx==1),method="breslow")
KM.1 <- survfit(NPfit, conf.int=.95, type="aalen")
plot(KM.1)
NPfit <- coxph(Surv(start,stop,event)~rx, data=bladder.data, subset=(rx==0),method="breslow")
KM.2 <- survfit(NPfit, conf.int=.95, type="aalen")
lines(KM.2, lty=3)
plot(KM.2)
lines(KM.1, lty=3)

NA.MF <- data.frame(time=c(0,KM$time), na=-log(c(1,KM$surv)))
plot(NA.MF, type="l", lty=1)

次に、上記と同等だが、それぞれの強度関数を推定してみる。 ここでは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)1:number+factor(rx)1:size+number:size+factor(rx)1: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で得たモデルと、ベースラインハザード関数、回帰係数、残差プロットなど、比較してみよ。

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