読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

Huzurbazar (2003) 打ち切りデータのヒストグラム(R による拡張版)

R 生存時間分析

アルゴリズム(algorithm)

打ち切りの種類

以前のこれ,
R による打ち切りデータのヒストグラム - 廿TT
は右打ち切りのみを想定していたが, 拡張した → plot a censored data histgram

生存時間分析の主な打ち切りには,

  • 右打ち切り
  • 左打ち切り
  • 2重打ち切り(両側打ち切り)
  • 区間打ち切り

がある.

打ち切られた観測値を z, 真の値を x として, 右打ち切りは z \le x, 左打ち切りは z \ge x, 2重打ち切り(両側打ち切り)は左打ち切りと右打ち切りが混在しているデータセットを指す.

区間打ち切りは  a< x \le b の範囲に x があることがわかっており, a, b のみが観測される状況を指す.

R のパッケージ survival では'right', 'left', 'interval', 'counting', 'interval2', 'mstate' を survival オブジェクトとして扱える.

'right', 'left', 'interval' はともかく, 'counting', 'mstate' がどういう状況で使われるのかは知らない.

Huzurbazar (2003) の打ち切りデータのヒストグラム

Huzurbazar (2003) のヒストグラムは右打ち切りのデータセットのみを想定している.
An Error Occurred Setting Your User Cookie

しかし, 生存率さえ求まればその差分を取り, ビンの幅で割ればヒストグラムが描ける.

下スライドの 6 ページから 7 ページを参照.

R が内側でどういう計算をしているのかはブラックボックスにしている.

引数(argument)

  • sf: survival パッケージの関数 survfit で生成したオブジェクトを入れる.
  • strata: stratum(層)の複数形. survfit オブジェクト内の strata を参照.
  • bw:ビンの幅.
  • digits: 小数点第何位まで表示するか.
  • col: 色.
  • main: グラフのメインタイトル.

例(example)

区間打ち切り(interval censored)

heart データは区間打ち切りされている.

> head(heart)
  start stop event        age      year surgery transplant id
1     0   50     1 -17.155373 0.1232033       0          0  1
2     0    6     1   3.835729 0.2546201       0          0  2
3     0    1     0   6.297057 0.2655715       0          0  3
4     1   16     1   6.297057 0.2655715       0          1  3
5     0   36     0  -7.737166 0.4900753       0          0  4
6    36   39     1  -7.737166 0.4900753       0          1  4
> head(Surv(heart$start, heart$stop, heart$event))
[1] ( 0,50]  ( 0, 6]  ( 0, 1+] ( 1,16]  ( 0,36+] (36,39] 

生存率を求める.

sf <- survfit(Surv(start, stop, event)~surgery, data=heart)
plot(sf,lty=2:1)
legend("topright",legend=c("surgery=0", "surgery=1"), lty=2:1)

f:id:abrahamcow:20141226114451p:plain

ヒストグラム描画.

par(mfrow=c(1,2))
a <- cdh(sf, strata=1, main="surgery=0")
b <- cdh(sf, strata=2, main="surgery=1")

f:id:abrahamcow:20141226114539p:plain

cdh 関数は total probability(生存率の死亡時刻ごとの差分)を返す.

> a
[1] 0.69013858 0.04041671 0.07654679 0.03214965 0.03214965 0.00000000 0.00000000
> b
[1] 0.3942308 0.2942308 0.0000000
> 

右打ち切り(right censored)

Gehan データは右打ち切りされている.

> data(gehan, package="MASS")
> head(gehan)
  pair time cens   treat
1    1    1    1 control
2    1   10    1    6-MP
3    2   22    1 control
4    2    7    1    6-MP
5    3    3    1 control
6    3   32    0    6-MP

生存率を求める.
おまけに survreg 関数でワイブル分布を仮定したパラメトリックなモデルにも当てはめる.

sf2 <- survfit(Surv(time,cens)~treat, data=gehan)
sr <- survreg(Surv(time,cens)~treat, data=gehan) #デフォルトが Weibull
plot(sf2, lty=1:2)
curve(1-pweibull(x,shape=1/sr$scale, scale=exp(coef(sr)[1])) ,add=TRUE, lty=1)
curve(1-pweibull(x,shape=1/sr$scale, scale=exp(coef(sr)[1] + coef(sr)[2])) ,add=TRUE, lty=2)
legend("topright",legend=c("6-MP", "cotrol"), lty=1:2)

f:id:abrahamcow:20141226115140p:plain

ヒストグラムを描いて, 密度関数の曲線を重ねてみる.

par(mfrow=c(1,2))
cdh(sf2, strata=1, main="6-MP")
#
xv <- seq(0,40,by=0.01)
points(xv, dweibull(xv,shape=1/sr$scale, scale=exp(coef(sr)[1] + coef(sr)[2])),
       type="l", lwd=2)
#
cdh(sf2, strata=2, main="control")
points(xv, dweibull(xv,shape=1/sr$scale, scale=exp(coef(sr)[1])), type="l", lwd=2)

f:id:abrahamcow:20141226115400p:plain

めちゃくちゃ当てはまり悪いな, これ. なにか間違えたかな?

今後の課題

気になるところはいろいろある.

いろいろ考えだすときりがないので, いったんここで一段落とした.

シミュレーションスタディもやってみたい.

参考文献(本文中に言及のないもの)

ああああ: R: survreg の分布一覧

生存時間解析

生存時間解析