KMMC plot

https://www.page-meeting.org/pdf_assets/4280-2012-06%20PAGE%20KMMC.pdf

Survival model の共変量の効果を視覚的に評価する方法。うまい方法です。
シミュレーションが必要になる。R コードを書いておきましょう。

library(survival)
library(ggplot2)

set.seed(2901)

summary(fit.1 <- survreg(Surv(day, e1c0) ~ age, ds, dist="exponential"))
#fit.1$coefficients

# observed KMMC ----
summary(ds)

ds.kmmc <- data.frame(
  day=seq(7, 365, by=7)
)
ds.kmmc$age.mean <- NA

for (d in seq(nrow(ds.kmmc))) {
  age.ge <- ds$age[ds$day >= ds.kmmc$day[d]]
  ds.kmmc$age.mean[d] <- mean(age.ge)
}

g1 <- ggplot(ds.kmmc, aes(day, age.mean)) +
  geom_line()
#g1

# simulation ----
nsim <- 100

.sim.kmmc <- ds.kmmc
.sim.kmmc$isim <- NA
sim.kmmc <- data.frame()
for (isim in seq(nsim)) {
  ds.sim <- ds
  ds.sim$day <- round(rexp(nrow(ds), rate=exp(-(fit.1$coefficients[1]+fit.1$coefficients[2]*ds$age))), 0)

  for (d in seq(nrow(.sim.kmmc))) {
    age.ge <- ds.sim$age[ds.sim$day >= .sim.kmmc$day[d]]
    .sim.kmmc$age.mean[d] <- mean(age.ge)
  }
  .sim.kmmc$isim <- isim
  sim.kmmc <- rbind(sim.kmmc, .sim.kmmc)
}

res <- as.data.frame(matrix(unlist(by(sim.kmmc, sim.kmmc$day, FUN=function(d1){
  quantile(d1$age.mean, probs=c(0.05, 0.5, 0.95))
})), ncol=3, byrow=T))
names(res) <- c("P05","P50","P95")
res$day <- ds.kmmc$day

g2 <- ggplot(ds.kmmc, aes(day)) +
  geom_ribbon(aes(ymin=P05, ymax=P95), res, fill="green3", alpha=0.5) +
  geom_line(aes(y=P05), res, colour="black") +
  geom_line(aes(y=P95), res, colour="black") +
  geom_line(aes(y=age.mean), colour="blue", size=2) +
  ylab('Mean')
g2
ggsave('KMMC.png', dpi=600)