We reconsider the data from Pike (Biometrics, 1982) on times from insult with the carcinogen DMBA to mortality from vaginal cancer in rats. The groups were distinguished by a pretreatment regimen. In Lab 4 where we introduced these data, we noted limitations of classical methods for incidence and cumulative incidence to analyze these data. Note the advantages of the alternatives considered here.\
Days to vaginal cancer mortality in rats:
Plus signed observations are censored. Censorship for the four rats may have occurred because they died of causes unrelated to the application of the carcinogen and were free of tumor at death, or they may simply not have developed tumor at the time of the data analysis.
grp1 <- c(143, 164, 188, 188, 190, 192, 206, 209, 213, 216, 220, 227, 230, 234, 246, 265, 304 , 244 )
grp2 <- c(142, 156, 163, 198, 205, 232, 232, 233, 233, 233, 233, 239, 240, 261, 280, 280, 296, 296, 323 , 344)
sum(grp1<230)/length(grp1)
sum(grp2<230)/length(grp2)
prop.test(x=c(sum(grp1<230),sum(grp2<230)), n=c(length(grp1), length(grp2)))
status
is coded 1 if the rat dies and 0 if she is censored. Calculate the incidence rates for vaginal cancer mortality.rat <- read.table("rat.dat", quote="\"", comment.char="")
names(rat) <- c("group", "days", "status")
sum(rat$status[which(rat$group==1)])/sum(rat$days[which(rat$group==1)])
sum(rat$status[which(rat$group==2)])/sum(rat$days[which(rat$group==2)])
#Load data
rat.pois <- glm(status~ group + offset(log(days)), data=rat ,family="poisson")
summary(rat.pois)
exp(confint(rat.pois))
Comment on the different conclusions with regard to statistical significance from the two approaches. Which approach is more principled?
Now obtain Kaplan-Meier survival estimates for these two treatment groups. In question 1, we compared survival at 230 days between groups. Make this comparison based on estimated survival functions, and comment on the advantages of this approach.
suppressMessages(library(survival))
rats <- survfit(Surv(days, status)~group, data=rat)
summary(rats)
suppressMessages(library(survminer))
rats.haz <- survfit(Surv(days, status)~group, data=rat, type="fleming")
ggsurvplot(rats.haz, conf.int = TRUE,
risk.table = TRUE, risk.table.col = "strata",
fun = "cumhaz", legend.labs = c("Group 1", "Group 2"),
break.time.by=20)
survdiff(Surv(days, status)~group, data=rat)