Comparing graph of stratified cox model to others in R

93 views Asked by At

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?

0

There are 0 answers