0

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?

Cryena
  • 33
  • 4
  • `basehaz`? Per `r` tag (hover or click to see): specify all non-base packages with `library()` calls. – Parfait Nov 21 '20 at 01:10
  • I don't understand your comment? The bfeed data is from the KMSurv package, is that what you are asking? – Cryena Nov 21 '20 at 09:05
  • Basically if we copy and paste your code as is in an empty R environment will it be runnable and reproducible? Likely it will err on `basehaz not found`. Give it a try. We need all `library` lines. A quick lookup indicates `library(survival)`. Please assume nothing. No need to specify base R packages: `base`, `stats`, `graphics`, etc. – Parfait Nov 21 '20 at 13:23
  • `exp(coef(fit))` is a relative risk estimate, NOT the predicted survival whereas `basehaz(fit)` is, so they are values on incommensurate scales. @Parfait is right; but you should also have included `library(survival)` as well as `library(KMsurv)` – IRTFM Nov 21 '20 at 15:51
  • Alright I specified better, thank you. I'm still learning how to ask questions, sorry! – Cryena Nov 22 '20 at 08:17

0 Answers0