## ----------------------------------------------------------------------------- evaluate_vignette <- requireNamespace("eha", quietly = TRUE) && requireNamespace("visreg", quietly = TRUE) && requireNamespace("splines", quietly = TRUE) knitr::opts_chunk$set(eval = evaluate_vignette) ## ----------------------------------------------------------------------------- # set.seed(12345) # library(casebase) # library(survival) # data(veteran) # table(veteran$status) # # evtimes <- veteran$time[veteran$status == 1] # hist(evtimes, nclass = 30, main = '', xlab = 'Survival time (days)', # col = 'gray90', probability = TRUE) # tgrid <- seq(0, 1000, by = 10) # lines(tgrid, dexp(tgrid, rate = 1.0/mean(evtimes)), # lwd = 2, lty = 2, col = 'red') ## ----------------------------------------------------------------------------- # veteran$prior <- factor(veteran$prior, levels = c(0, 10), labels = c("no","yes")) # veteran$celltype <- factor(veteran$celltype, # levels = c('large', 'squamous', 'smallcell', 'adeno')) # veteran$trt <- factor(veteran$trt, levels = c(1, 2), labels = c("standard", "test")) ## ----------------------------------------------------------------------------- # library(eha) # y <- with(veteran, Surv(time, status)) # # model1 <- weibreg(y ~ karno + diagtime + age + prior + celltype + trt, # data = veteran, shape = 1) # summary(model1) ## ----------------------------------------------------------------------------- # model2 <- weibreg(y ~ karno + diagtime + age + prior + celltype + trt, # data = veteran, shape = 0) # summary(model2) ## ----------------------------------------------------------------------------- # model3 <- coxph(y ~ karno + diagtime + age + prior + celltype + trt, # data = veteran) # summary(model3) ## ----------------------------------------------------------------------------- # # create popTime object # pt_veteran <- popTime(data = veteran) # # class(pt_veteran) # # # plot method for objects of class 'popTime' # plot(pt_veteran) ## ----------------------------------------------------------------------------- # model4 <- fitSmoothHazard(status ~ time + karno + diagtime + age + prior + # celltype + trt, data = veteran, ratio = 100) # summary(model4) ## ----results='show', R.options=list(max.print=1)------------------------------ # library(visreg) # plot(model4, # hazard.params = list(alpha = 0.05)) ## ----------------------------------------------------------------------------- # absRisk4 <- absoluteRisk(object = model4, time = 90) # mean(absRisk4) # # ftime <- veteran$time # mean(ftime <= 90) ## ----------------------------------------------------------------------------- # model5 <- fitSmoothHazard(status ~ log(time) + karno + diagtime + age + prior + # celltype + trt, data = veteran, ratio = 100) # summary(model5) ## ----------------------------------------------------------------------------- # # Fit a spline for time # library(splines) # model6 <- fitSmoothHazard(status ~ bs(time) + karno + diagtime + age + prior + # celltype + trt, data = veteran, ratio = 100) # summary(model6) # # str(absoluteRisk(object = model6, time = 90)) ## ----------------------------------------------------------------------------- # linearRisk <- absoluteRisk(object = model4, time = 90, newdata = veteran) # splineRisk <- absoluteRisk(object = model6, time = 90, newdata = veteran) # # plot.default(linearRisk, splineRisk, # xlab = "Linear", ylab = "Splines", pch = 19) # abline(a = 0, b = 1, lty = 2, lwd = 2, col = 'red') ## ----echo = FALSE------------------------------------------------------------- # table_coef <- cbind(coefficients(model3), # coefficients(model4)[-(1:2)], # coefficients(model5)[-(1:2)], # coefficients(model6)[-(1:4)]) # knitr::kable(table_coef, format = "html", digits = 4, # col.names = c("Cox model", # "CB linear", # "CB log-linear", # "CB splines")) ## ----------------------------------------------------------------------------- # # define a specific covariate profile # new_data <- data.frame(trt = "test", # celltype = "adeno", # karno = median(veteran$karno), # diagtime = median(veteran$diagtime), # age = median(veteran$age), # prior = "no") # # # calculate cumulative incidence using casebase model # smooth_risk <- absoluteRisk(object = model4, # time = seq(0,300, 1), # newdata = new_data) # # cols <- c("#8E063B","#023FA5") # # # cumulative incidence function for the Cox model # plot(survfit(model3, newdata = new_data), # xlab = "Days", ylab = "Cumulative Incidence (%)", fun = "event", # xlim = c(0,300), conf.int = F, col = cols[1], # main = sprintf("Estimated Cumulative Incidence (risk) of Lung Cancer\ntrt = test, celltype = adeno, karno = %g,\ndiagtime = %g, age = %g, prior = no", median(veteran$karno), median(veteran$diagtime), # median(veteran$age))) # # # add casebase curve with legend # plot(smooth_risk, add = TRUE, col = cols[2], gg = FALSE) # legend("bottomright", # legend = c("semi-parametric (Cox)", "parametric (casebase)"), # col = cols, # lty = c(1, 1), # bg = "gray90") ## ----------------------------------------------------------------------------- # smooth_risk <- absoluteRisk(object = model4, # time = seq(0,300, 1), # newdata = new_data, # type = "survival") # # plot(survfit(model3, newdata = new_data), # xlab = "Days", ylab = "Survival Probability (%)", # xlim = c(0,300), conf.int = F, col = cols[1], # main = sprintf("Estimated Survival Probability of Lung Cancer\ntrt = test, celltype = adeno, karno = %g,\ndiagtime = %g, age = %g, prior = no", median(veteran$karno), median(veteran$diagtime), # median(veteran$age))) # # # add casebase curve with legend # plot(smooth_risk, add = TRUE, col = cols[2], gg = FALSE) # legend("topright", # legend = c("semi-parametric (Cox)", "parametric (casebase)"), # col = cols, # lty = c(1, 1), # bg = "gray90") ## ----echo=FALSE, eval=TRUE---------------------------------------------------- print(sessionInfo(), locale = F)