[R program]生存時間解析

Rのsurvivalパッケージだけでいろいろな生存時間解析ができます。

  • カップランマイヤー（Kaplan-Meier）曲線を描く
  • Cox回帰をする
  • Log-Rank検定をする

この3つが生存時間解析の主軸になりますが、それ以外にも次のようなものができます。

  • 条件付ロジスティック
  • Aalen additive model（ノンパラのモデル）
  • ケースコホート研究のハザードの計算
  • 人年の計算
  • アメリカ国勢調査データから調整した期待死亡率を計算（SMRみたいな感じ？）
  • 周辺モデル
  • Frailtyモデル
  • 比例ハザード性の検定
  • pスプラインによる共変量のスムージング（pspline）
  • リッジ回帰
  • Harrington-Fleming検定
  • 競合リスク（Competing risk）の計算

とりあえず今日実行した中からいくつか結果を紹介します。

  • KM曲線

f:id:isseing333:20100702160758j:image

  • 比例ハザード性のチェック（水平になっていれば大丈夫）

f:id:isseing333:20100702160933j:image

f:id:isseing333:20100702161044j:image


SASではLog-Logプロットとか描けるんだけど、survivalでは見当たらなかったなぁ。

あとKM曲線の信頼区間とか。

自分で計算して出すのか、他のパッケージにあるのかね。

まぁ、predictとかで推定値が出るから簡単に計算できるのだろうけど。


コードは以下のようになります↓

#------生存時間解析
library(survival)

example(aareg)				#Aalen additive model（ノンパラの生存時間）
example(cch)				#ケースコホート研究のハザード（case-cohort）
example(coxph.detail)			#Cox回帰結果の詳細
example(pyears)				#人年の計算
example(survexp)


#---カップランマイヤー曲線
leukemia.surv <- survfit(Surv(time, status) ~ x, data = aml)
plot(leukemia.surv, lty = 2:3)
legend(100, .9, c("Maintenance", "No Maintenance"), lty = 2:3)
title("Kaplan-Meier Curves\nfor AML Maintenance Study")
lsurv2 <- survfit(Surv(time, status) ~ x, aml, type='fleming')
plot(lsurv2, lty=2:3, fun="cumhaz",
xlab="Months", ylab="Cumulative Hazard")

#---Kaplan-Meier曲線の修正
fit <- survfit(Surv(time, status) ~ sex, pbc,subset=1:312)			#胆汁性肝硬変
plot(fit, mark.time=FALSE, xscale=365.25, xlab='Years', ylab='Survival')
lines(fit[1], lwd=2, xscale=365.24)

#ratetableでアメリカ国勢調査データから調整した期待死亡率を計算してくれる
efit <- survexp(~ ratetable(sex=sex,age=age*365.35,year=as.Date('1979/1/1')) + sex, data=pbc, times=(0:24)*182)
temp <- lines(efit, lty=2, xscale=365.24, lwd=2:1)				#試験参加日がないから試験の中央値でまかなう
text(temp, c("Male", "Female"), adj= -.1)					#labels just past the ends
title(main="Primary Biliary Cirrhosis, Observed and Expected")			#アメリカ国勢調査から計算した期待値



#---条件付ロジスティック（層を指定するが層のパラメータは推定しない→層の調整みたいなイメージかな）
clogit(case~spontaneous+induced+strata(stratum),data=infert)
a <-glm(case~spontaneous+induced+factor(stratum), data=infert, family="binomial")	#普通のロジスティック回帰
a						#層の数が多くてまともにやったら推定できない
exp(a$coeff)

#---周辺モデルとFrailtyモデル
marginal.model <- coxph(Surv(time, status) ~ rx + cluster(litter), rats)	#clusterで相関のある対象者を指定
frailty.model <- coxph(Surv(time, status) ~ rx + frailty(litter), rats)		#S Ripatti and J Palmgren, Estimation of multivariate frailty models using penalized partial likelihood, Biometrics, 56:1016-1022, 2000.

#Frailtyモデル
coxph(Surv(time, status) ~ age + frailty(inst, df=4), lung)

#Litter effects for the rats data（Litter効果って何だ？）
rfit2a <- survreg(Surv(time, status) ~ rx + frailty.gaussian(litter, df=13, sparse=FALSE), rats )
rfit2b <- survreg(Surv(time, status) ~ rx + frailty.gaussian(litter, df=13, sparse=TRUE), rats )


#---比例ハザード性のチェック
#P. Grambsch and T. Therneau (1994), Proportional hazards tests and diagnostics based on weighted residuals. Biometrika, 81, 515-26.
fit <- coxph(Surv(futime, fustat) ~ age + ecog.ps, data=ovarian)
temp <- cox.zph(fit)
print(temp)					# display the results、ハザード比が一定→推定バラメータが時間に対して一定（水平）になるはず
plot(temp)					# plot curves　　　　　→有意になると傾きが水平じゃない→比例ハザード性が崩壊

#cox.zphのプロット
vfit <- coxph(Surv(time,status) ~ trt + factor(celltype) + karno + age, data=veteran, x=TRUE)
temp <- cox.zph(vfit)

plot(temp[5])					#tempをみると変数karnoは比例ハザードではない
abline(0, 0, lty=3)
abline(lm(temp$y[,5] ~ temp$x)$coefficients, lty=4, col=3)			#直線回帰を当てはめると確かに上昇している
title(main="VA Lung Study")


#---pスプラインによる共変量スムージング（penalized spline）
lfit6 <- survreg(Surv(time, status)~pspline(age, df=2), cancer)
plot(cancer$age, predict(lfit6), xlab='Age', ylab="Spline prediction")
title("Cancer Data")

(fit0 <- coxph(Surv(time, status) ~ ph.ecog + age, cancer))
(fit1 <- coxph(Surv(time, status) ~ ph.ecog + pspline(age,3), cancer))

#---リッジ回帰（ridge regression）
coxph(Surv(futime, fustat) ~ rx + ridge(age, ecog.ps, theta=1), ovarian)


#---層別Cox回帰、IDでクラスター
bladder1 <- bladder[bladder$enum < 5, ]
coxph(Surv(stop, event) ~ (rx + size + number) * strata(enum) + cluster(id), bladder1)

#---生存時間解析に関連する分布
x <- seq(.1, 3, length=30)
haz <- dsurvreg(x, 2, 3)/ (1-psurvreg(x, 2, 3))
plot(x, haz, log='xy', ylab="Hazard")		#line with slope (1/scale -1)


#---Harrington-Fleming検定
#Harrington, D. P. and Fleming, T. R. (1982). A class of rank test procedures for censored survival data. Biometrika 69, 553-566.
survdiff(Surv(futime, fustat) ~ rx,data=ovarian)
survdiff(Surv(futime, fustat) ~ rx,data=ovarian, rho=1)				#rho=0がデフォルト→log-rank検定、1は修正Gehan-Wilcoxon検定


#---Cox回帰後に分散分析（devianceをみてるらしい、モデル間の比較・検定をやってくれてるのか？）
fit <- coxph(Surv(futime, fustat) ~ resid.ds *rx + ecog.ps, data = ovarian)
anova(fit)
fit2 <- coxph(Surv(futime, fustat) ~ resid.ds +rx + ecog.ps, data=ovarian)
anova(fit2,fit)


#---競合リスク（Competing risk curves (cumulative incidence)）
fit1 <- survfit(Surv(stop, event=='progression') ~1, data=mgus1, subset=(start==0))
fit2 <- survfit(Surv(stop, status) ~1, data=mgus1, subset=(start==0), etype=event)	#competing risks

# CI curves are always plotted from 0 upwards, rather than 1 down
plot(fit2, fun='event', xscale=365.25, xmax=7300, mark.time=FALSE, col=2:3, xlab="Years post diagnosis of MGUS")
lines(fit1, fun='event', xscale=365.25, xmax=7300, mark.time=FALSE, conf.int=FALSE)
text(10, .4, "Competing Risk: death", col=3)
text(16, .15,"Competing Risk: progression", col=2)
text(15, .30,"KM:prog")

ページTOPへ