0

我正在查看 R 中的 bfeed 数据集：

library(survival)
library(KMsurv)
data(bfeed)
dput(head(bfeed))
structure(list(duration = c(16L, 1L, 4L, 3L, 36L, 36L), delta = c(1L, 
1L, 0L, 1L, 1L, 1L), race = c(1L, 1L, 1L, 1L, 1L, 1L), poverty = c(0L, 
0L, 0L, 0L, 0L, 0L), smoke = c(0L, 1L, 0L, 1L, 1L, 0L), alcohol = c(1L, 
0L, 0L, 1L, 0L, 0L), agemth = c(24L, 26L, 25L, 21L, 22L, 18L), 
    ybirth = c(82L, 85L, 85L, 85L, 82L, 82L), yschool = c(14L, 
    12L, 12L, 9L, 12L, 11L), pc3mth = c(0L, 0L, 0L, 0L, 0L, 0L
    )), row.names = c(NA, 6L), class = "data.frame")

我为它制作了这三个模型：

bf.surv <- Surv(duration,delta)
racef <- factor(race, labels = c("white","black","other"))
bf.moda <- aareg(bf.surv~smoke*alcohol+racef+poverty)
bf.mod1 <- coxph(bf.surv~smoke*alcohol+racef+poverty)
bf.mod2 <- coxph(bf.surv~strata(smoke)+strata(alcohol)+racef+poverty)

我可以很好地绘制前两个：

bf.new <- data.frame(rbind(c(smoke=0,alcohol=0),c(smoke=1,alcohol=0),c(smoke=0,alcohol=1),c(smoke=1,alcohol=1)),poverty=0,racef="white")
bf.fit <- survfit(bf.mod1,bf.new)
t <- c(0,bf.moda$times) # ties
temp <- rbind(0,bf.moda$coef)
temp.beta <- rowsum(temp,t,reorder=F)
Beta <- apply(temp.beta,2,cumsum)
T.s <- unique(t)

plot(T.s, exp(-Beta[,1]), type="s", xlab="Time", ylab="Estimated Survival Fuction", main="Comparing Models - No smoking, No Alcohol")
lines(bf.fit[1], type="s", col=2, conf.int=FALSE)

但是当我尝试绘制分层 cox 模型时，它并没有出现在图表上。我确定我错过了一些简单的东西，但我似乎无法自己抓住它。这是我正在尝试的：

bf.bh   <- basehaz(bf.mod2)
cond1   <- bf.bh$strata=="smoke=0,alcohol=0"
mod2.t1 <- bf.bh$time[cond1]
mod2.s1 <- bf.bh$hazard[cond1]

plot(T.s, exp(-Beta[,1]), type="s", xlab="Time", ylab="Estimated Survival Fuction", main="Comparing Models - No smoking, No Alcohol")
lines(bf.fit[1], type="s", col=2, conf.int=FALSE)
lines(mod2.t1,mod2.s1, type="s", col=3, conf.int=FALSE)

它不会在图表上显示为一条线吗？

4

0 回答 0