library(survival)

d <- read.csv("http://p.tl/ZDLT") # 白血病の生存時間データ

res <- coxph(Surv(weeks,status)~group,d)

# Coxモデルの適合性を評価するためのCox-Snell残差

# 原点を通る傾き1の直線に沿えば適合性は高い

# クライン&メシュベルガー, 生存時間解析, 11.2節

r.cs <- d$status - residuals(res,type="martingale") #Cox-Snell残差

r.surv <- survfit(Surv(r.cs,status)~1,type="fleming-harrington",d)

plot(r.surv$time,-log(r.surv$surv),xlab="Cox-Snell残差", ylab="累積ハザード",type="S")

title(main="Coxモデルの適合性を評価するためのCox-Snell残差")

abline(0,1,col="red")

# 共変量の関数形の決定:マルチンゲール残差

# クライン&メシュベルガー, 生存時間解析, 11.3節

# 新たな変数とマルチンゲール残差のトレンドから関数形を見る

r.m  <- residuals(res,type="martingale") # マルチンゲール残差

plot(d$logwbc,r.m)

lines(lowess(x=d$logwbc,y=r.m), col="red")

# モデルにlogwbcを線形で追加

res <- coxph(Surv(weeks,status)~group+logwbc,d)

# 最大対数尤度比残差

# クライン&メシュベルガー, 生存時間解析, 11.5節

r.d  <- residuals(res,type="deviance") #

lp <- predict(res,type="lp")

plot(lp,r.d)

# 各観測値の影響:b - b[j]

# b : 全データでの推定値

# b[j] : j番目の観測値を除いたデータでの推定値

par(mfrow=c(2,2))

r.b  <- residuals(res,type="dfbeta") # 各観測値の影響度

plot(1:nrow(d),r.b[,1],main="group"); abline(h=0,lty="dotdash")

plot(1:nrow(d),r.b[,2],main="logwbc"); abline(h=0,lty="dotdash")

r.bs  <- residuals(res,type="dfbetas") # 各観測値の影響度(標準化)

plot(1:nrow(d),r.bs[,1],main="group"); abline(h=0,lty="dotdash")

plot(1:nrow(d),r.bs[,2],main="logwbc"); abline(h=0,lty="dotdash")