### R code from vignette source 'Introduction.Rnw' ################################################### ### code chunk number 1: Introduction.Rnw:85-89 ################################################### options(width=80,useFancyQuotes="UTF-8") library(rstpm2) ################################################### ### code chunk number 2: Introduction.Rnw:234-242 ################################################### brcancer <- transform(brcancer, recyear=rectime / 365.24) fit <- stpm2(Surv(recyear,censrec==1)~hormon, data=brcancer, df=4) summary(fit) ## utility eform.coxph <- function(object) exp(cbind(coef(object),confint(object))) fit.cox <- coxph(Surv(recyear,censrec==1)~hormon, data=brcancer) rbind(cox=eform(fit.cox), eform(fit)[2,,drop=FALSE]) ################################################### ### code chunk number 3: Introduction.Rnw:250-255 ################################################### plot(fit, newdata=data.frame(hormon=0), xlab="Time since diagnosis (years)") lines(fit, newdata=data.frame(hormon=1), lty=2) lines(survfit(Surv(recyear,censrec==1)~hormon, data=brcancer), col="blue", lty=1:2) legend("topright", c("PH hormon=0","PH hormon=1","KM hormon=0","KM hormon=1"), lty=1:2, col=c("black","black","blue","blue")) ################################################### ### code chunk number 4: Introduction.Rnw:258-262 ################################################### plot(fit,newdata=data.frame(hormon=1), type="hazard", xlab="Time since diagnosis (years)", ylim=c(0,0.3)) lines(fit, newdata=data.frame(hormon=0), col=2, lty=2, type="hazard") legend("topright", c("hormon=1","hormon=0"),lty=1:2,col=1:2,bty="n") ################################################### ### code chunk number 5: Introduction.Rnw:269-273 ################################################### plot(fit,newdata=data.frame(hormon=0), type="hdiff", exposed=function(data) transform(data, hormon=1), main="hormon=1 compared with hormon=0", xlab="Time since diagnosis (years)") ################################################### ### code chunk number 6: Introduction.Rnw:275-279 ################################################### plot(fit,newdata=data.frame(hormon=0), type="sdiff", exposed=function(data) transform(data, hormon=1), main="hormon=1 compared with hormon=0", xlab="Time since diagnosis (years)") ################################################### ### code chunk number 7: Introduction.Rnw:290-299 ################################################### brcancer <- transform(brcancer, recyear=rectime / 365.24) fit <- stpm2(Surv(recyear,censrec==1)~hormon, data=brcancer, link.type="AH") summary(fit) confint(fit) plot(fit, newdata=data.frame(hormon=0), xlab="Time on study (years)") lines(fit, newdata=data.frame(hormon=1), lty=2) lines(survfit(Surv(recyear,censrec==1)~hormon, data=brcancer), col="blue", lty=1:2) legend("topright", c("AH hormon=0","AH hormon=1","KM hormon=0","KM hormon=1"), lty=1:2, col=c("black","black","blue","blue")) ################################################### ### code chunk number 8: Introduction.Rnw:303-311 ################################################### fit <- stpm2(Surv(recyear,censrec==1)~1, data=brcancer, link.type="AH", smooth.formula=~ns(sqrt(recyear),df=3)+hormon:ns(recyear,df=3)) summary(fit) plot(fit, newdata=data.frame(hormon=0), xlab="Time on study (years)") suppressWarnings(lines(fit, newdata=data.frame(hormon=1), lty=2)) lines(survfit(Surv(recyear,censrec==1)~hormon, data=brcancer), col="blue", lty=1:2) legend("topright", c("AH hormon=0","AH hormon=1","KM hormon=0","KM hormon=1"), lty=1:2, col=c("black","black","blue","blue")) ################################################### ### code chunk number 9: Introduction.Rnw:334-338 ################################################### options(width=80,useFancyQuotes="UTF-8") library(rstpm2) ################################################### ### code chunk number 10: Introduction.Rnw:340-353 ################################################### popmort2 <- transform(rstpm2::popmort,exitage=age,exityear=year,age=NULL,year=NULL) colon2 <- within(rstpm2::colon, { status <- ifelse(surv_mm>120.5,1,status) tm <- pmin(surv_mm,120.5)/12 exit <- dx+tm*365.25 sex <- as.numeric(sex) exitage <- pmin(floor(age+tm),99) exityear <- floor(yydx+tm) ##year8594 <- (year8594=="Diagnosed 85-94") }) colon2 <- merge(colon2,popmort2) ################################################### ### code chunk number 11: Introduction.Rnw:357-362 ################################################### fit0 <- stpm2(Surv(tm,status %in% 2:3)~I(year8594=="Diagnosed 85-94"), data=colon2, bhazard=colon2$rate, df=5) ################################################### ### code chunk number 12: Introduction.Rnw:364-371 ################################################### summary(fit <- stpm2(Surv(tm,status %in% 2:3)~I(year8594=="Diagnosed 85-94"), data=colon2, bhazard=colon2$rate, df=5,cure=TRUE)) predict(fit,head(colon2),se.fit=TRUE) ################################################### ### code chunk number 13: Introduction.Rnw:389-396 ################################################### newdata.eof <- data.frame(year8594 = unique(colon2$year8594), tm=10) predict(fit0, newdata.eof, type="fail", se.fit=TRUE) predict(fit, newdata.eof, type="fail", se.fit=TRUE) predict(fit, newdata.eof, type="haz", se.fit=TRUE) ################################################### ### code chunk number 14: Introduction.Rnw:400-415 ################################################### tms=seq(0,10,length=301)[-1] plot(fit0,newdata=data.frame(year8594 = "Diagnosed 85-94", tm=tms), ylim=0:1, xlab="Time since diagnosis (years)", ylab="Relative survival") plot(fit0,newdata=data.frame(year8594 = "Diagnosed 75-84",tm=tms), add=TRUE,line.col="red",rug=FALSE) ## warnings: Predicted hazards less than zero for cure plot(fit,newdata=data.frame(year8594 = "Diagnosed 85-94",tm=tms), add=TRUE,ci=FALSE,lty=2,rug=FALSE) plot(fit,newdata=data.frame(year8594="Diagnosed 75-84",tm=tms), add=TRUE,rug=FALSE,line.col="red",ci=FALSE,lty=2) legend("topright",c("85-94 without cure","75-84 without cure", "85-94 with cure","75-84 with cure"), col=c(1,2,1,2), lty=c(1,1,2,2), bty="n") ################################################### ### code chunk number 15: Introduction.Rnw:422-439 ################################################### plot(fit0,newdata=data.frame(year8594 = "Diagnosed 85-94", tm=tms), ylim=c(0,0.5), type="hazard", xlab="Time since diagnosis (years)",ylab="Excess hazard") plot(fit0,newdata=data.frame(year8594 = "Diagnosed 75-84", tm=tms), type="hazard", add=TRUE,line.col="red",rug=FALSE) plot(fit,newdata=data.frame(year8594 = "Diagnosed 85-94", tm=tms), type="hazard", add=TRUE,ci=FALSE,lty=2,rug=FALSE) plot(fit,newdata=data.frame(year8594="Diagnosed 75-84", tm=tms), type="hazard", add=TRUE,rug=FALSE,line.col="red",ci=FALSE,lty=2) legend("topright",c("85-94 without cure","75-84 without cure", "85-94 with cure","75-84 with cure"), col=c(1,2,1,2), lty=c(1,1,2,2), bty="n") ################################################### ### code chunk number 16: Introduction.Rnw:445-458 ################################################### newdata.eof <- data.frame(year8594 = unique(colon2$year8594), tm=10) test <- predictnl(fit, function(object,newdata=NULL) { lp1 <- predict(object, newdata.eof[1,], type="link") lp2 <- predict(object, newdata.eof[2,], type="link") lp1-lp2 }) with(test, c(fit=fit, se.fit=se.fit, statistic=fit/se.fit, p=2*pnorm(abs(fit/se.fit), lower.tail=FALSE)))