#======================================================
#	11 生存時間解析
#======================================================
#パッケージのロード
library(survival) #生存時間分析のため
library(MASS)     #stepAICのため
library(KMsurv)   #データhodgのため

#======================================================
#	12.1 基本的な概念
#======================================================
#---------------<生存関数>---------------
#主要な関数のプロット（ワイブル分布の例）

time <- 1:420    #35年間×12か月

#図11.1
windows(10, 7)
op <- par(mfrow=c(2, 2), mar=c(4.1, 4.1, 3.1, 1.1))
S.t <- 1 - pweibull(time, shape=2.5, scale=250)
plot(time, S.t, type="l", ylab="S(t)", main="生存関数")
f.t <- dweibull(time, shape=2.5, scale=exp(5.5))
plot(time, f.t, type="l", ylab="f(t)", main="確率密度関数")
h.t <- f.t/S.t
plot(time, h.t, type="l", ylab="h(t)", main="ハザード関数")
H.t <- -log(S.t)
plot(time, H.t, type="l", ylab="H(t)", main="累積ハザード関数")
par(op)

#---------------<打ち切り>---------------

#生存時間データの形式
head(aml)
Surv(aml$time, aml$status)

#時間依存共変量がある場合
head(heart)
Surv(heart$start, heart$stop, heart$event)

#======================================================
#	12.2 ノンパラメトリック法
#======================================================
#---------------<Kaplan-Meier推定量>---------------

#推定(卵巣がん)
head(ovarian)
ofit<-survfit(Surv(futime, fustat) ~ rx, data=ovarian, conf.type="plain")
summary(ofit)

#プロット(卵巣がん)
#図11.2
plot(ofit, xscale=365.25, lty=1:2,
       xlab="年", ylab="生存曲線")
legend(x=2, y=1, legend=c("rx = 1", "rx = 2"), lty = 1:2, bty="n")

#肺がん
lung.fit <- survfit(Surv(time, status) ~ sex, data=lung)
#図11.3
plot(lung.fit, col=c("blue", "red"), lty=1:2, xlab="経過日数", 
       ylab="生存曲線", main="進行した肺がん患者の性別生存曲線")
legend(800, 1, legend=c("男性", "女性"), col=c("blue", "red"), 
       lty=1:2, bty="n")

#---------------<2以上の生存曲線の比較>---------------

#卵巣がん
survdiff(Surv(futime, fustat)~rx, data=ovarian)

#肺がん
survdiff(Surv(time, status) ~ sex, data=lung)

#Peto検定
survdiff(Surv(futime, fustat)~rx, data=ovarian, rho=1)

#======================================================
#	12.3 パラメトリック法
#======================================================
#---------------<ワイブル分布>---------------

#lungデータの概観

summary(lung)  #データの要約
idx <- !is.na(lung$ph.ecog) #欠損値でないデータ

#survregの当てはめ

lung.reg <- survreg(Surv(time, status) ~ sex + factor(ph.ecog), 
        data=lung[idx, ])
summary(lung.reg)

#モデル行列と回帰係数による線形予測子の計算

X <- model.matrix(lung.reg)            #モデル行列
head(X)                                #最初のいくつか
(beta <- coef(lung.reg))               #回帰係数
X.beta <- X %*% beta                   #線形予測子の計算
head(X.beta)                           #最初のいくつか

#対数生存時間を予測する線形予測子の確認

pred.lp <- predict(lung.reg, type="lp") #線形予測子（対数生存時間）
head(pred.lp)                           #最初のいくつか（上記と同じ）
head(exp(pred.lp))                      #指数変換（原スケール）

#生存時間の予測

pred.res <- predict(lung.reg)           #応答（生存時間）：既定値
head(pred.res)                          #同じ

#生存曲線の予測

#Kaplan-Meier推定量の当てはめとプロット

lung.fit <- survfit(Surv(time, status) ~ sex + ph.ecog, 
        data=lung[idx, ], subset=(sex==1 & ph.ecog==1))
#図11.4
plot(lung.fit, conf.int=F, col="steelblue", xlab="経過日数", 
       ylab="生存曲線", main="進行した肺がん患者の生存曲線（男性、ECOG=1）")

#ワイブル分布による生存曲線の当てはめとプロット

pct <- 1:98 / 100                                          #生存率
pred.quant <- predict(lung.reg, type="quantile", 
        newdata=data.frame(sex=1, ph.ecog=1), p=pct, se=T) #生存曲線の予測
fit <- pred.quant$fit                      #推定値
upper <- fit + 1.96 * pred.quant$se.fit    #95％上限
lower <- fit - 1.96 * pred.quant$se.fit    #95％下限

lines(fit, 1-pct, col="tomato", lwd=3)     #推定値を描く
lines(upper, 1-pct, lty=2)                 #95％上限を描く
lines(lower, 1-pct, lty=2)                 #95％下限を描く
legend(650, 1, 
       legend=c("Kaplan-Meier推定量", "ワイブル当てはめ", "95％信頼区間"),
	bty="n", lty=c(1, 1, 2), col=c("steelblue", "tomato", "black"),
       lwd=c(1, 3, 1))

#関数pweibullを使用して確認

days <- seq(0, 1000, length=98)         #経過日数   
alpha <- 1 / lung.reg$scale             #shapeパラメータ
x.0 <- c(1, 1, 1 , 0, 0)                #共変量（男性, ECOG=1）
beta <- coef(lung.reg)                  #回帰係数
lp.0 <- beta %*% x.0                    #線形予測子
scale.0 <- exp(lp.0)                    #線形予測子の指数変換(= scale)
S.0 <- 1 - pweibull(days, shape = alpha, scale = scale.0) #生存曲線
lines(days, S.0, col="skyblue", lwd=3)  #生存曲線の上書き

#ハザード関数

f.0 <- dweibull(days, shape = alpha, scale = scale.0) #確率密度関数
h.0 <- f.0 / S.0                                      #ハザード関数

#図11.5
plot(days, h.0, type="l",xlab="経過日数", ylab="ハザード率", 
       main="進行した肺がん患者のハザード関数（男性、ECOG=1）",
	col="steelblue", lwd=3) 

#---------------<対数ロジスティック分布>---------------

#survregの当てはめ

lung.llogis <- survreg(Surv(time, status) ~ sex + factor(ph.ecog), 
        dist="loglogistic", data=lung[idx, ])
summary(lung.llogis)

#生存曲線の予測

#Kaplan-Meier推定量
#図11.6
plot(lung.fit, conf.int=F, col="steelblue", xlab="経過日数", 
       ylab="生存曲線", main="進行した肺がん患者の生存曲線（男性、ECOG=1）")

#対数ロジスティック分布による生存曲線の当てはめとプロット

pct <- 1:98 / 100                                          #生存率
pred.quant <- predict(lung.llogis, type="quantile", 
        newdata=data.frame(sex=1, ph.ecog=1), p=pct, se=T) #生存曲線の予測
fit <- pred.quant$fit                      #推定値
upper <- fit + 1.96 * pred.quant$se.fit    #95％上限
lower <- fit - 1.96 * pred.quant$se.fit    #95％下限

lines(fit, 1-pct, col="tomato", lwd=3)     #推定値を描く
lines(upper, 1-pct, lty=2)                 #95％上限を描く
lines(lower, 1-pct, lty=2)                 #95％下限を描く
legend(650, 1, 
       legend=c("Kaplan-Meier推定量", "対数ロジスティック", "95％信頼区間"),
	bty="n", lty=c(1, 1, 2), col=c("steelblue", "tomato", "black"),
       lwd=c(1, 3, 1))

#生存関数S(t)の定義

sllogis <- function(t, alpha, lambda) {
        ans <- 1 / (1 + (lambda * t)^alpha)	
}

#関数sllogisを使用して確認

days <- seq(0, 1000, length=98)         #経過日数
alpha <- 1 / lung.llogis$scale          #alphaパラメータ
x.0 <- c(1, 1, 1 , 0, 0)                #共変量（男性, ECOG=1）
beta <- coef(lung.llogis)               #回帰係数
lp.0 <- beta %*% x.0                    #線形予測子
lambda <- exp(-lp.0)                    #lambdaパラメータ	
S.0 <- sllogis(days, alpha, lambda)     #生存曲線
lines(days, S.0, col="skyblue", lwd=3)  #生存曲線の上書き

#ハザード関数h(t)の定義

hllogis <- function(t, alpha, lambda) {
        S <- sllogis(t, alpha, lambda)
        num <- lambda * alpha * (lambda * t)^(alpha - 1)
        ans <- S * num		
}

#ハザード関数のプロット

h.0 <- hllogis(days, alpha, lambda)      #ハザード関数
#図11.7
plot(days, h.0, type="l",xlab="経過日数", ylab="ハザード率", 
       main="対数ロジスティック分布によるハザード関数（男性、ECOG=1）",
	col="steelblue", lwd=3) 

#---------------<対数正規分布>---------------

#survregの当てはめ
lung.lnorm <- survreg(Surv(time, status) ~ sex + factor(ph.ecog), 
        dist="lognormal", data=lung[idx, ])
summary(lung.lnorm)

#生存曲線の予測
#Kaplan-Meier推定量
#図11.8
plot(lung.fit, conf.int=F, col="steelblue", xlab="経過日数", 
       ylab="生存曲線", main="進行した肺がん患者の生存曲線（男性、ECOG=1）")

#対数正規分布による生存曲線の当てはめとプロット

pct <- 1:98 / 100                                          #生存率
pred.quant <- predict(lung.lnorm, type="quantile", 
        newdata=data.frame(sex=1, ph.ecog=1), p=pct, se=T) #生存曲線の予測
fit <- pred.quant$fit                      #推定値
upper <- fit + 1.96 * pred.quant$se.fit    #95％上限
lower <- fit - 1.96 * pred.quant$se.fit    #95％下限

lines(fit, 1-pct, col="tomato", lwd=3)     #推定値を描く
lines(upper, 1-pct, lty=2)                 #95％上限を描く
lines(lower, 1-pct, lty=2)                 #95％下限を描く
legend(650, 1, 
       legend=c("Kaplan-Meier推定量", "対数正規", "95％信頼区間"),
	bty="n", lty=c(1, 1, 2), col=c("steelblue", "tomato", "black"),
       lwd=c(1, 3, 1))

#関数plnormを使用して確認

days <- seq(0, 1000, length=98)         #経過日数   
alpha <- lung.lnorm$scale               #sdlogパラメータ
x.0 <- c(1, 1, 1 , 0, 0)                #共変量（男性, ECOG=1）
beta <- coef(lung.lnorm)                #回帰係数
lp.0 <- beta %*% x.0                    #線形予測子(= meanlog)
S.0 <- 1 - plnorm(days, meanlog = lp.0, sdlog = alpha) #生存曲線
lines(days, S.0, col="skyblue", lwd=3)  #生存曲線の上書き

#ハザード関数

f.0 <- dlnorm(days, meanlog = lp.0, sdlog = alpha) #確率密度関数
h.0 <- f.0 / S.0                                   #ハザード関数
#図11.9
plot(days, h.0, type="l",xlab="経過日数", ylab="ハザード率", 
       main="対数正規分布によるハザード関数（男性、ECOG=1）",
	col="steelblue", lwd=3) 

#AIC
p <- length(lung.reg$coef)        #パラメータ数
-2 * lung.reg$loglik[2] + 2 * (p + 2)    #ワイブル分布
-2 * lung.llogis$loglik[2] + 2 * (p + 2) #対数ロジスティック分布
-2 * lung.lnorm$loglik[2] + 2 * (p + 2)  #対数正規分布
 
#======================================================
#	11.4 Cox比例ハザードモデル
#======================================================
#---------------<比例ハザード性の仮定>---------------

#coxphの当てはめ
lung.cox <- coxph(Surv(time, status) ~ sex + factor(ph.ecog), 
        data=lung[idx, ])
summary(lung.cox)

#補対数・対数プロット
#図11.10
op <- par(mfrow=c(1, 2), mar=c(4.1, 4.1, 3.1, 1.1))
lung.fit2 <- survfit(Surv(time, status) ~ sex, data=lung[idx, ])
plot(lung.fit2, fun="cloglog", xlab="経過日数", 
    ylab="log(-log(生存曲線))", main="性別", lty=1:2, 
    col=c("steelblue", "tomato"), lwd=2)
legend(5, 1, legend=c("男性", "女性"), lty=1:2, 
        col=c("steelblue", "tomato"), lwd=2, bty="n")

lung.fit3 <- survfit(Surv(time, status) ~ ph.ecog, data=lung[idx, ])
plot(lung.fit3, fun="cloglog", xlab="経過日数",
    ylab="log(-log(生存曲線))", main="ECOGスコア", lty=1:3,
    col=c("steelblue", "tomato", "brown"), lwd=2)
legend(5, 1, legend=c("スコア０", "スコア１", "スコア２"), lty=1:3, 
        col=c("steelblue", "tomato", "brown"), lwd=2, bty="n")
par(op)

#stepAICによるモデル選択

lung.0 <- na.omit(lung)                  #欠損値を除く
lung.0$ph.ecog2 <- lung.0$ph.ecog        #ph.ecogをコピー
lung.0$ph.ecog2[lung.0$ph.ecog2==3] <- 2 #スコア３を２とまとめる。
lung.cox.full <- coxph(Surv(time, status) ~ (age + sex + factor(ph.ecog2)
        + ph.karno + pat.karno + meal.cal + wt.loss)^2, data=lung.0)
lung.cox.best <- stepAIC(lung.cox.full)  #変数選択
summary(lung.cox.best)

#生存曲線の予測
#Kaplan-Meier推定量
#図11.11
plot(lung.fit, conf.int=F, col="steelblue", xlab="経過日数", 
       ylab="生存曲線", main="進行した肺がん患者の生存曲線（男性、ECOG=1）")

#coxphによる生存曲線の当てはめとプロット

idx2 <- lung.0$sex==1 & lung.0$ph.ecog2==1 #該当する患者を示すインデックス
lung.cox.surv <- survfit(lung.cox.best, newdata=data.frame(sex=1, 
        ph.ecog2=1, 
        ph.karno=mean(lung.0$ph.karno[idx2]),
        pat.karno=mean(lung.0$pat.karno[idx2]),
        meal.cal=mean(lung.0$meal.cal[idx2]),
        wt.loss=mean(lung.0$wt.loss[idx2])),
        se.fit=T)
lines(lung.cox.surv, conf.int=T, col=c("tomato", "black", "black"),
	lwd=c(3, 1, 1))
legend(650, 1, 
       legend=c("Kaplan-Meier推定量", "coxph", "95％信頼区間"),
       bty="n", lty=c(1, 1, 2), col=c("steelblue", "tomato", "black"),
       lwd=c(1, 3, 1))

#層化Coxモデル

lung.cox.strata <- coxph(Surv(time, status) ~ (age + factor(ph.ecog2)
        + ph.karno + pat.karno + meal.cal + wt.loss)^2 + strata(sex), 
        data=lung.0)
lung.cox.best2 <- stepAIC(lung.cox.strata) #変数選択
summary(lung.cox.best2)

#層化Coxモデルの生存曲線の予測
#Kaplan-Meier推定量
#図11.12
plot(lung.fit, conf.int=F, col="steelblue", xlab="経過日数", 
       ylab="生存曲線", main="進行した肺がん患者の生存曲線（男性、ECOG=1）")

#層化Coxモデルによる生存曲線の当てはめとプロット

lung.cox.surv2 <- survfit(lung.cox.best2, newdata=data.frame(
        age=mean(lung.0$age[idx2]),
        sex=1, 
        ph.ecog2=1, 
        ph.karno=mean(lung.0$ph.karno[idx2]),
        pat.karno=mean(lung.0$pat.karno[idx2]),
        wt.loss=mean(lung.0$wt.loss[idx2])),
        se.fit=T)
lung.cox.surv2$strata    #層化因子
lines(lung.cox.surv2$time[1:91], lung.cox.surv2$surv[1:91],
        col="tomato", lwd=3)                  #推定値を描く
lines(lung.cox.surv2$time[1:91], 
        lung.cox.surv2$upper[1:91], lty=2)    #95％上限を描く
lines(lung.cox.surv2$time[1:91], 
        lung.cox.surv2$lower[1:91], lty=2)    #95％下限を描く
legend(650, 1, 
       legend=c("Kaplan-Meier推定量", "層化Coxモデル", "95％信頼区間"),
       bty="n", lty=c(1, 1, 2), col=c("steelblue", "tomato", "black"),
       lwd=c(1, 3, 1))

#-----------<残差診断>-----------
#-----------<マルチンゲール残差>-----------
data(hodg)                    #データのロード
hodg.cox <- coxph(Surv(time, delta) ~ (gtype + dtype)^2 + score, 
        data = hodg)          #当てはめ
resid.hodg <- resid(hodg.cox) #マルチンゲール残差

#待ち時間なし
#図11.13
plot(hodg$wtime, resid.hodg, xlab="待ち時間（月）", 
        ylab="マルチンゲール残差")  #マルチンゲール残差のプロット
abline(h=0, lty=2)                  #残差０
lines(loess.smooth(hodg$wtime, resid.hodg), 
        col="steelblue", lwd=3)     #平滑化曲線

#pspline(待ち時間)
hodg.cox2 <- coxph(Surv(time, delta) ~ (gtype + dtype)^2 + score
        + pspline(wtime, df=0), data = hodg) #pspline
resid.hodg2 <- resid(hodg.cox2)     #マルチンゲール残差
#図11.14
plot(hodg$wtime, resid.hodg2, xlab="待ち時間（月）", 
        ylab="マルチンゲール残差")  #マルチンゲール残差のプロット
abline(h=0, lty=2)                  #残差０
lines(loess.smooth(hodg$wtime, resid.hodg2), 
        col="steelblue", lwd=3)     #平滑化曲線

#-----------<スコア残差、dfbeta、dfbetas>-----------

dfbetas.hodg2 <- resid(hodg.cox2, type="dfbetas")      #dfbetas
dimnames(dfbetas.hodg2)[[2]] <- names(coef(hodg.cox2)) #列名を付ける。
head(dfbetas.hodg2)                                    #確認

#図11.15
op <- par(mfrow=c(2, 2), mar=c(4.1, 4.1, 3.1, 1.1))
plot(dfbetas.hodg2[ ,1],　xlab="観測値", type="h",
        ylab="dfbetas", main="gtype", ylim=c(-2, 2))
plot(dfbetas.hodg2[ ,2], xlab="観測値", type="h", 
        ylab="dfbetas", main="dtype", ylim=c(-2, 2))
plot(dfbetas.hodg2[ ,3], xlab="観測値", type="h",
        ylab="dfbetas", main="score", ylim=c(-2, 2))
plot(dfbetas.hodg2[ ,21], xlab="観測値", type="h", 
        ylab="dfbetas", main="gtype:dtype", ylim=c(-2, 2))
par(op)

#-----------<Schoenfeld残差>-----------
lung.zph <- cox.zph(lung.cox.best2)
lung.zph

#図11.16
op <- par(mfrow=c(2, 2), mar=c(4.1, 4.1, 3.1, 1.1))
plot(lung.zph[1], main="age"); abline(h=0, lty=2)
plot(lung.zph[4], main="ph.karno"); abline(h=0, lty=2)
plot(lung.zph[5], main="pat.karno"); abline(h=0, lty=2)
plot(lung.zph[6], main="wt.loss"); abline(h=0, lty=2)
par(op)

#------------<多重イベント>------------
#------------<多重イベントの概念図>------------
#図11.17
op <- par(mfrow=c(2, 2), mar=c(4.1, 4.1, 3.1, 1.1))
#AG
plot(c(0, 10), c(0, 10), type="n", axes=F, xlab="", ylab="", 
        main="AG", cex.main=3)
symbols(c(4, 8), c(5, 5), inches=F, 
        rectangles=matrix(c(2, 2), nrow=2, ncol=2, byrow=T), add=T)
text(c(4 ,8), c(5 ,5), c("0", "1+"), cex=2)
arrows(x0=5, x1=7, y0=5, y1=5, length=0.15)
segments(x0=8, x1=8, y0=6, y1=8)
segments(x0=8, x1=1, y0=8, y1=8)
segments(x0=1, x1=1, y0=8, y1=5)
arrows(x0=1, x1=3, y0=5, y1=5, length=0.15)

#PWP
plot(c(0, 10), c(0, 10), type="n", axes=F, xlab="", ylab="", 
        main="PWP", cex.main=3)
symbols(c(1, 3.5, 6, 8.5), rep(5 ,4), inches=F, 
        rectangles=matrix(c(2, 2), nrow=4, ncol=2, byrow=T), add=T)
arrows(x0=c(2, 4.5, 7), x1=c(2.5, 5, 7.5), y0=5, y1=5, length=0.15)
text(c(1, 3.5, 6, 8.5), rep(5 ,4), c("0", "1", "2", "3+"), cex=2)

#WLW
plot(c(0, 10), c(0, 10), type="n", axes=F, xlab="", ylab="", 
        main="WLW", cex.main=3)
symbols(c(2, 8, 8, 8), c(5, 8, 5, 2), inches=F, 
        rectangles=matrix(c(2, 2), nrow=4, ncol=2, byrow=T), add=T)
arrows(x0=rep(3, 3), x1=rep(7, 3), y0=5, y1=c(8, 5, 2), length=0.15)
text(c(2, 8, 8, 8), c(5, 8, 5, 2), c("0", "1", "2", "3+"), cex=2)

par(op)

#------------<CGD>------------

head(cgd)
summary(cgd)
table(cgd$enum)
cgd$enum2 <- cgd$enum         #コピー
cgd$enum2[cgd$enum > 3] <- 3  #３超を３とする。

#------------<AGモデル>------------

cgd.ag <- coxph(Surv(tstart, tstop, status) ~ treat + age
       + steroids + cluster(id), data=cgd)
summary(cgd.ag)

#------------<PWPモデル>------------

cgd.pwp <- coxph(Surv(tstop - tstart, status) ~ treat + age
       + steroids + cluster(id) + strata(enum2), data=cgd)
summary(cgd.pwp)
#図11.18
pwp.surv <- survfit(cgd.pwp)
plot(pwp.surv, fun="cloglog", lty=1:3)
legend(2, 0, legend=c("enum 1", "enum 2", "enum 3+"), lty=1:3)
 
#======================================================
#	11.5 randomForest
#======================================================
library(randomForestSRC)
help(package=randomForestSRC)

data(wihs)
summary(wihs)

wihs.rf <- rfsrc(Surv(time, status) ~ ., data=wihs, 
    nsplit = 3, ntree = 500, seed=1)
wihs.rf

#図11.19
plot(wihs.rf)

#図11.20
plot.competing.risk(wihs.rf)

#図11.21
cif <- wihs.rf$cif
Time <- wihs.rf$time.interest
idu <- wihs$idu
cif.haart <- cbind(apply(cif[,,1][idu == 0,], 2, mean), 
    apply(cif[,,1][idu == 1,], 2, mean))
cif.aids  <- cbind(apply(cif[,,2][idu == 0,], 2, mean), 
    apply(cif[,,2][idu == 1,], 2, mean))
matplot(Time, cbind(cif.haart, cif.aids), type = "l",
    lty = c(1,2,1,2), col = c(4, 4, 2, 2), lwd = 3,
    ylab = "Cumulative Incidence")
legend("topleft",
    legend = c("HAART (Non-IDU)", "HAART (IDU)", 
    "AIDS (Non-IDU)", "AIDS (IDU)"),
    lty = c(1,2,1,2), col = c(4, 4, 2, 2), lwd = 3, cex=1.35)

#図11.22
plot.variable(wihs.rf)
