I'm looking at the bfeed dataset in R:
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")
and I made these three models for it:
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)
and I can plot the first two fine:
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)
But then when I try to plot the stratified cox model it doesn't show up on the graph. I'm sure I'm missing something simple, but I can't seem to catch it on my own. Here's what I'm trying:
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)
It doesn't show up as a line on the graph though?