## ----survtram-pkgs, echo = FALSE, results = "hide", message = FALSE, warning = FALSE---- # required packages pkgs <- c("mlt", "tram", "trtf", "SparseGrid", "ATR", "tramME", "multcomp", "coin", "TH.data", "survival", "colorspace", "xtable", "english") pkgs <- sapply(pkgs, require, character.only = TRUE) ## ----fail, results = "asis", echo = FALSE-------------------------------- if (any(!pkgs)) { cat(paste("Package(s)", paste(names(pkgs)[!pkgs], collapse = ", "), "not available, stop processing.", "\\end{document}\n")) knitr::knit_exit() } if (!interactive() && .Platform$OS.type != "unix") { cat(paste("Vignette only compiled under Unix alikes.", "\\end{document}\n")) knitr::knit_exit() } ## ----setup, echo = FALSE, results = "hide", message = FALSE, warning = FALSE---- set.seed(290875) ## plotting xlab <- "Time (in days)" lxlab <- paste0(xlab, " on log-scale") ylabS <- "Probability of survival" ylablHaz <- "Log-cumulative hazard" ylabcumHR <- expression(Lambda[1](t)/Lambda[0](t)) ylimS <- c(0, 1) ylimHR <- c(0, 1.6) q <- 0:2204 xlim <- c(0, max(q)) lwd <- 1.3 ## color acol <- sequential_hcl(6, "BluYl")[1:5] col <- acol[c(2, (length(acol)) - 1)] lcol <- lighten(col, amount = .4) ## lighten color for overlaid lines ## aux perm_test_biv.stram <- function(object, seed = 1) { stopifnot(inherits(object, "stram")) fixed <- c(trt = 0, scl = 0) lhs <- object$call[[2]][[3]] if (!(length(lhs) == 3 & lhs[[2]] == lhs[[3]])) stop("Bivariate score perm test not applicable") names(fixed) <- names(coef(object)) m0 <- update(object, fixed = fixed) ## uncond. model r <- resid(m0, what = "shifting") rs <- resid(m0, what = "scaling") set.seed(seed) formula <- as.formula(paste("r + rs ~", lhs[[2]])) pvalue(independence_test(formula, data = m0$data)) } ## formatting big.mark <- "'" frmt0 <- round frmt <- function(digits, x, math = FALSE) { if (!is.numeric(x)) return(x) ret <- formatC(round(x, digits), digits = digits, format = "f", big.mark = big.mark) if (math) ret <- paste("$", ret, "$") if (is.matrix(x)) { ret <- matrix(ret, nrow = nrow(x)) dimnames(ret) <- dimnames(x) } ret } frmt1 <- function(x, math = FALSE) frmt(1, x = x, math = math) frmt2 <- function(x, math = FALSE) frmt(2, x = x, math = math) frmt3 <- function(x, math = FALSE) frmt(3, x = x, math = math) ## logLik frmtll <- function(x, math = FALSE, mark = FALSE) { if (!inherits(x, "logLik") && !is.numeric(x) && all(!is.na(x))) x <- logLik(x) if (is.na(x)) return("") ret <- frmt2(abs(x), math = FALSE) if (x < 0) ret <- paste0(ifelse(math, "$-$", "-"), ret) if (mark) ret <- paste0("{\\color{darkgray}", ret, "}") ret } ## data load(system.file("rda", "Primary_endpoint_data.rda", package = "TH.data")) ## randomization arm levs <- levels(CAOsurv$randarm) trt <- with(CAOsurv, paste0("randarm", levs[2], collapse = "")) nd1 <- data.frame(randarm = factor(levs, levels = levs)) ## strata CAOsurv$strat <- with(CAOsurv, interaction(strat_t, strat_n)) slevs <- levels(CAOsurv$strat) nd2 <- data.frame(randarm = nd1$randarm[1], strat = factor(slevs, levels = slevs)) ## for pretty legends lslevs <- gsub("\\.", " : ", slevs) lslevs <- gsub("cT4", "cT4 ", lslevs) ## id CAOsurv$id <- factor(1:nrow(CAOsurv)) # ### lattice # trellis.par.set( # list( # plot.symbol = list(col = 1, pch = 20, cex = 0.7), # box.rectangle = list(col = 1), # box.umbrella = list(lty = 1, col = 1), # strip.background = list(col = "white") # ) # ) # # ltheme <- canonical.theme(color = FALSE) ## in-built B&W theme # ltheme$strip.background$col <- "transparent" ## change strip bg # lattice.options(default.theme = ltheme) ## knitr library("knitr") knitr::opts_chunk$set(results = "hide", echo = FALSE, purl = TRUE, tidy = FALSE, size = "small", error = FALSE, warning = FALSE, message = FALSE, fig.scap = NA, fig.align = "center", fig.width = 6, fig.height = 3.3, out.width = NULL, out.height = NULL, dev = c("pdf", "postscript")) opts_knit$set(global.par = TRUE) knitr::render_sweave() # use Sweave environments knitr::set_header(highlight = "") # do not \usepackage{Sweave} options(prompt = "R> ", continue = "+ ", useFancyQuotes = FALSE) # JSS style options(width = 74) ## format frmtI <- function(x, math = FALSE) { if (is.character(x)) return(x) ret <- trimws(formatC(x, format = "fg", width = 7, big.mark = big.mark)) if (x < 0) ret <- paste0(ifelse(math, "$-$", "-"), ret) if (math) ret <- paste("$", ret, "$") if (is.matrix(x)) { ret <- matrix(ret, nrow = nrow(x)) dimnames(ret) <- dimnames(x) } ret } frmtN <- function(x, bound = 10, all = TRUE) { ## => all because of consistency ret <- round(x) if (all | x <= bound) return(as.character(english::english(ret))) ret } frmtd <- function(date) format(date, format = "%b~%-d") frmtD <- function(date) format(date, format = "%B~%-d") toUpper <- function(x) gsub("\\b([[:lower:]])([[:lower:]]+)", "\\U\\1\\L\\2", x, perl = TRUE) ## print results frmtp <- function(pv) paste("$p$ =", frmt3(pv)) frmtCI <- function(x, phantom = FALSE, math = FALSE) { if (!isTRUE(nrow(x) > 1)) phantom <- FALSE if (phantom) phantom <- any(x[,2] < 0) ## switch all to phantom if one is negative FUN <- function(x) { ciL <- frmt3(x[1], math = FALSE) ciL <- ifelse(math & phantom & ciL > 0, paste0("$\\phantom{-}", ciL, "$"), paste0("$", ciL, "$")) ciU <- frmt3(x[2], math = FALSE) ciU <- ifelse(math & phantom & ciU > 0, paste0("$\\phantom{-}", ciU, "$"), paste0("$", ciU, "$")) paste(ciL, "to", ciU) } if (is.matrix(x)) ret <- apply(x, 1, FUN) else ret <- FUN(x) return(ret) } ## table with results from "tram" ret.tram <- function(object, seed = 25, math = TRUE) { mod0 <- object if (inherits(object, "Survreg")) object <- as.mlt(object) ## Estimate cf <- coef(object) ibs <- grep("Bs", names(cf)) if (length(ibs) > 1) cf <- cf[-ibs] if (inherits(object, "Compris")) cf <- cf[c(grep("^Event_EoI", names(cf), value = TRUE), "Event_DepC.Event_EoI.(Intercept)")] if (inherits(object, "fmlt")) cf <- c(cf, coef(object, addparm = TRUE)) ncf <- names(cf) ## SE if (inherits(object, "fmlt")) vc <- vcov(object, addparm = TRUE) else vc <- vcov(object) se <- sqrt(diag(vc))[ncf] ## CIs if (inherits(object, "tram") & !inherits(object, "stram")) { st <- score_test(object) ci <- st$conf.int nci <- "95\\%-Score CI" } else { ciWald <- function(cf, se) cf + c(-1, 1) * qnorm(1-0.05/2) * se ci <- t(sapply(1:length(cf), function(i) ciWald(cf[i], se[i]))) nci <- "95\\%-Wald CI" } # ci <- confint(object)[ncf, ] ncf <- sapply(ncf, function(x) switch(x, "(Intercept)" = "$\\eparm_1$", "log(iDFS)" = "$\\eparm_2$", "randarm5-FU + Oxaliplatin" = "$\\eshiftparm$", "Event_EoI.Event_EoI.randarm5-FU + Oxaliplatin" = "$\\eshiftparm_\\rY$", "scl_randarm5-FU + Oxaliplatin" = "$\\escaleparm$", "gamma1" = "$\\mparm$", "Event_DepC.Event_EoI.(Intercept)" = "$\\lparm$", "logrho" = "$\\log(\\lparm)$")) ret.cfs <- cbind("Coefficient" = ncf, "Estimate" = frmt3(cf, math = math), "Std. Error" = frmt3(se), "95\\%-CI" = frmtCI(ci, phantom = TRUE, math = math)) colnames(ret.cfs)[which(colnames(ret.cfs) == "95\\%-CI")] <- nci if (inherits(object, "tramME")) ret.cfs <- rbind(ret.cfs, c("$\\tau^2$", frmt3(VarCorr(object)[[1]]$var), "", "")) ll <- logLik(object) object <- mod0 if (inherits(object, "tram")) { p.lrt <- summary(object)$p.value if (inherits(object, "stram")) { p.pst <- perm_test_biv.stram(object) ret.mod <- cbind("Bivariate Permutation Score Test" = frmtp(p.pst), " " = NA) } else { if (inherits(mod0, "Survreg")) st <- score_test(mod0) p.st <- st$p.value p.pt <- (pt <- perm_test(object))$p.value ret.mod <- cbind("Score Test" = frmtp(p.st),"Permutation Score Test" = frmtp(p.pt)) } ret.mod <- cbind(cbind("Log-Likelihood" = frmtll(ll, math = math), "Likelihood Ratio Test" = frmtp(p.lrt)), ret.mod) } if (inherits(object, c("tramME", "mtram", "fmlt", "mmlt"))) ret.mod <- cbind("Log-Likelihood" = frmtll(ll, math = math), " " = NA, " " = NA, " " = NA) ret.mod <- rbind(colnames(ret.mod), ret.mod) ns <- colnames(ret.cfs) ret <- rbind(ret.cfs, ret.mod) colnames(ret) <- ns return(ret) } print.tram <- function(object) { ret <- ret.tram(object) cat("\\begin{center}") cat("\n") # cat("\\small") cat("\n") print(xtable(ret, align = "rrrrr"), add.to.row = list(pos = list(nrow(ret) - 2, nrow(ret) - 1, nrow(ret)), command = c("\\\\[-1.2ex] \\toprule ", "\\midrule ", "\\\\[-2ex] ")), booktabs = TRUE, floating = FALSE, sanitize.text.function = function(x){x}, include.rownames = FALSE, scale = .8) cat("\n") cat("\\end{center}") } ## Formatting is.neg <- function(x) x < 0 lHR <- function(model, frmt = function(x) {frmt3(x, math = FALSE)}) { stopifnot(inherits(model, c("tram", "mlt", "tramME"))) cf <- coef(model)[trt] if (inherits(model, "tram")) cf <- c(1, -1)[model$negative + 1] * cf frmt(cf) } HR <- function(model, frmt = function(x) {frmt3(x, math = FALSE)}) frmt(exp(lHR(model, frmt = function(x) {x}))) ## ----pars, include = FALSE---------------------------------------------- par_main <- expression(par(mgp = c(2.5, 1, 0), mar = c(4, 4, 1.5, 4), las = 1)) par_surv <- expression(par(mgp = c(2.5, 1, 0), mar = c(6, 6, .5, 4), las = 1)) ## ----packages, echo = FALSE--------------------------------------------- library("tram") ## ----risk-tab----------------------------------------------------------- risktab <- function(ti, st) { ## time-index and survival time nrisk <- NULL for (t in ti) nrisk <- c(nrisk, sum(st >= t)) return(nrisk) } plot.risktab <- function(tvar, ti = seq(min(q), max(q), by = 500), cex = .8, at = -450) { mtext(levs[1], 1, line = 4, at = at, cex = cex) mtext(risktab(ti, CAOsurv[CAOsurv$randarm == levs[1], tvar]), side = 1, line = 4, at = ti, cex = cex) mtext(levs[2], 1, line = 5, at = at, cex = cex) mtext(risktab(ti, CAOsurv[CAOsurv$randarm == levs[2], tvar]), side = 1, line = 5, at = ti, cex = cex) } ## ----surv-OS------------------------------------------------------------ surv_OS <- survfit(OS ~ randarm, data = CAOsurv) ## KM ## ----surv-iDFS---------------------------------------------------------- surv_iDFS <- survfit(iDFS ~ randarm, data = CAOsurv) ## Turnbull ## ----CAO-table, results = 'hide'---------------------------------------- tab <- xtabs( ~ strat + randarm, data = CAOsurv) tab <- rbind(tab, "Total" = colSums(tab)) ## ----surv-iDFS-plot----------------------------------------------------- eval(par_surv) plot(surv_iDFS, ylim = ylimS, xlim = xlim, col = lcol, lwd = lwd, xlab = xlab, ylab = ylabS) legend("bottomright", legend = levs, col = col, bty = "n", lty = 1, lwd = 1, cex = .8) plot.risktab(tvar = "iDFStime") ## ----surv-OS-plot------------------------------------------------------- plot(surv_OS, ylim = ylimS, xlim = xlim, col = lcol, lwd = lwd, xlab = xlab, ylab = ylabS) legend("bottomright", legend = levs, col = col, bty = "n", lty = 1, lwd = 1, cex = .8) plot.risktab(tvar = "OStime") ## ----WEI-model-fit, echo = FALSE, cache = TRUE-------------------------- mw <- Survreg(iDFS ~ randarm, data = CAOsurv, dist = "weibull") ## ----WEI-summary, cache = TRUE, results = "hide", fig.show = 'hide'----- summary(mw) coef(mw, as.survreg = TRUE) ## same interpretation as "survreg" score_test(mw) perm_test(mw) # plot(as.mlt(mw), type = "survivor", newdata = nd1, col = col) ## ----COX-model-fit, echo = FALSE, cache = TRUE-------------------------- mc <- Coxph(iDFS ~ randarm, data = CAOsurv, log_first = TRUE) ## ----COX-summary, cache = TRUE, results = "hide", fig.show = "hide"----- summary(mc) score_test(mc) perm_test(mc) # plot(as.mlt(mc), type = "survivor", newdata = nd1, col = col) ## ----COX-lHaz, echo = FALSE--------------------------------------------- ## confband object <- as.mlt(mc) newdata <- nd1[1,, drop = FALSE] K <- 20 cheat <- 200 y <- variable.names(object, "response") q <- mkgrid(object, n = K)[[y]] q[1] <- q[1] + 1 ## quick fix for log_first nd <- newdata[rep(1, length(q)), , drop = FALSE] nd[[y]] <- q X <- model.matrix(object$model$model, data = nd) cb <- confint(multcomp::glht(multcomp::parm(coef(object), vcov(object)), linfct = X))$confint q <- mkgrid(object, n = cheat)[[y]] q[1] <- q[1] + 1 ## quick fix for log_first nd <- newdata[rep(1, length(q)), , drop = FALSE] nd[[y]] <- q X <- model.matrix(object$model$model, data = nd) cb <- confint(multcomp::glht(multcomp::parm(coef(object), vcov(object)), linfct = X), calpha = attr(cb, "calpha"))$confint cb <- cbind(q, cb) ## ----COX-lHaz-plot------------------------------------------------------ eval(par_main) plot(cb[, "q"], cb[, "Estimate"], log = "x", type = "n", xlab = lxlab, ylab = ylablHaz, xlim = xlimlHaz <- range(cb[, "q"]), ylim = range(cb[, -1])) polygon(c(cb[, "q"], rev(cb[, "q"])), c(cb[, "lwr"], rev(cb[, "upr"])), border = NA, col = rgb(.1, .1, .1, .1)) lines(cb[, "q"], cb[, "Estimate"], lwd = lwd) ## ----STRAT-model-fit, cache = TRUE-------------------------------------- mcst <- Coxph(iDFS | strat ~ randarm, data = CAOsurv, log_first = TRUE) ## ----STRAT-summary, cache = TRUE, results = "hide", fig.show = "hide"---- summary(mcst) score_test(mcst) perm_test(mcst) ## ----STRAT-lHaz-plot---------------------------------------------------- plot(as.mlt(mcst), newdata = nd2, q = q, type = "trafo", log = "x", lty = lty <- 1:4, xlab = lxlab, ylab = ylablHaz, xlim = xlimlHaz, col = 1, lwd = lwd) legend("bottomright", legend = lslevs, title = "Stratum", lty = lty, lwd = lwd, col = 1, bty = "n") ## ----SCOX-model-fit, cache = TRUE--------------------------------------- mcs <- Coxph(iDFS ~ randarm | randarm, data = CAOsurv, log_first = TRUE) ## ----SCOX-summary, cache = TRUE, results = "hide", fig.show = "hide"---- summary(mcs) confint(mcs) perm_test_biv.stram(mcs) # plot(as.mlt(mcs), type = "survivor", newdata = nd1, col = col) ## ----SCOX-survplot, fig.show = "hide"----------------------------------- eval(par_surv) plot(surv_iDFS, ylim = ylimS, xlim = xlim, col = lcol, lwd = lwd, xlab = xlab, ylab = ylabS) legend("bottomright", legend = levs, col = col, bty = "n", lty = 1, lwd = 1, cex = .8) plot.risktab(tvar = "iDFStime") plot(as.mlt(mcs), type = "survivor", newdata = nd1, col = col, add = TRUE) ## ----SCOX-HR-plot, echo = FALSE----------------------------------------- qHR <- seq(50, max(q), by = 1) cumhaz <- predict(mcs, type = "cumhazard", newdata = nd1, q = qHR) cumhr <- unname(cumhaz[, 2] / cumhaz[, 1]) plot(qHR, cumhr, type = "l", ylab = ylabcumHR, xlab = xlab, ylim = ylimHR, xlim = xlimHR <- range(qHR), lwd = lwd) abline(h = exp(coef(mc)), lty = 2, lwd = 1) ## constant HR abline(h = 1, lty = 3) ## HR = 1 ## ----TCOX-model-fit, cache = TRUE--------------------------------------- mcv <- Coxph(iDFS | randarm ~ 1, data = CAOsurv, log_first = TRUE) ## ----TCOX-summary, cache = TRUE, results = "hide", fig.show = 'hide'---- logLik(mcv) ## ----TCOX-HR, cache = TRUE---------------------------------------------- mcv <- as.mlt(mcv) ## grid (was n = 500, qmvnorm failed) s <- mkgrid(mcv, 39) s$iDFS <- s$iDFS[s$iDFS >= min(xlimHR) & s$iDFS <= max(xlimHR)] nd3 <- expand.grid(s) ## confint K <- model.matrix(mcv$model, data = nd3) Kyes <- K[nd3$randarm == levels(nd3$randarm)[2],] Kyes[,grep("Intercept", colnames(Kyes))] <- 0 gh <- glht(parm(coef(mcv), vcov(mcv)), Kyes) ci <- exp(confint(gh)$confint) coxy <- s$iDFS ## confint for constant HR ci2 <- exp(confint(mc)) ## ----TCOX-HR-plot------------------------------------------------------- plot(coxy, ci[, "Estimate"], ylim = ylimHR, type = "n", xlim = xlimHR, xlab = xlab, ylab = ylabcumHR) polygon(c(coxy, rev(coxy)), c(ci[,"lwr"], rev(ci[, "upr"])), border = NA, col = rgb(.1, .1, .1, .1)) lines(coxy, ci[, "Estimate"], lty = 1, lwd = lwd) ## constant HR polygon(c(coxy[c(1, length(coxy))], rev(coxy[c(1, length(coxy))])), rep(ci2, c(2, 2)), border = NA, col = rgb(.1, .1, .1, .1)) abline(h = exp(coef(mc)), lty = 2, lwd = 1) ## HR = 1 abline(h = 1, lty = 3) ## ----DEPCENS-preproc, echo = FALSE-------------------------------------- ## DepC: loss of follow-up (everyone else is admin censored) Mail TH 23-06-12 patnr_lofu <-c(1012, 2003, 3002, 3003, 6018, 7001, 7003, 7005, 7008, 7012, 10003, 10012, 11018, 12003, 12014, 13028, 14002, 15001, 16001, 16004, 16005, 16007, 16009, 18016, 18025, 21011, 21013, 21014, 21022, 21023, 21026, 21027, 21029, 21043, 22003, 23008, 24008, 24021, 25001, 25004, 25005, 25006, 26005, 26018, 27005, 27030, 27034, 29002, 30006, 30011, 31003, 31004, 31005, 34001, 35011, 35014, 36017, 41004, 42001, 42003, 42005, 42007, 42010, 44004, 44005, 45002, 45003, 45009, 45011, 46003, 49001, 49003, 49011, 49012, 49015, 50001, 50003, 50004, 50007, 50011, 52004, 54004, 56006, 56008, 59002, 59005, 68001, 70010, 71002, 73009, 74004, 75002, 75004, 75005, 80003, 81001, 84005, 84007, 86002) ilofu <- with(CAOsurv, which(patnr %in% patnr_lofu)) CAOsurv$DepCevent <- CAOsurv$OSevent CAOsurv$DepCevent <- factor(as.numeric(CAOsurv$DepCevent), levels = 0:2, labels = c("AdminC", "EoI", "DepC")) CAOsurv$DepCevent[ilofu] <- "DepC" ## ----DEPCENS-table, results = 'hide'------------------------------------ CAOsurv$nDepCevent <- factor(as.character(CAOsurv$DepCevent), levels = c("AdminC", "EoI", "DepC"), labels = c("Administrative censoring", "Event of interest", "Loss of follow-up")) tab <- xtabs(~ nDepCevent + randarm, data = CAOsurv) tab ## ----DEPCENS-model-fit, cache = TRUE------------------------------------ md <- Coxph(Surv(OStime, event = DepCevent) ~ randarm, data = CAOsurv, log_first = TRUE) ## ----DEPCENS-summary, cache = TRUE, results = "hide", fig.show = 'hide'---- summary(md) confint(md) ## ----COXME-install, echo = TRUE, eval = FALSE--------------------------- # install.packages("tramME") # library("tramME") ## ----COXME-load--------------------------------------------------------- library("tramME") ## ----COXME-model-fit, cache = TRUE-------------------------------------- mcME <- CoxphME(iDFS ~ randarm + (1 | Block), data = CAOsurv, log_first = TRUE) ## ----COXME-summary, cache = TRUE, results = "hide", fig.show = 'hide'---- summary(mcME) confint(mcME) ## ----COXME-margsurv, eval = FALSE--------------------------------------- # ## computationally intensive # mod <- mcME # # ## A function to evaluate the joint cdf of the response and the random effects: # ## Takes a vector of random effect and covariates values, evaluates the conditional # ## distribution at these values and multiplies it with the pdf of the random effects # jointCDF <- function(re, nd, mod) { # nd <- nd[rep(1, length(re)), ] # nd$Block <- seq(nrow(nd)) ## to take vector-valued REs # pr <- predict(mod, newdata = nd, ranef = re, type = "distribution") * # dnorm(re, 0, sd = sqrt(varcov(mod)[[1]][1, 1])) # c(pr) # } # ## Marginalize the joint cdf by integrating out the random effects # ## using adaptive quadrature # marginalCDF <- function(nd, mod) { # nd$cdf <- integrate(jointCDF, lower = -Inf, upper = Inf, nd = nd, mod = mod)$value # nd # } # ## Set up the grid on which we evaluate the marginal distribution # nd <- expand.grid(iDFS = 1:max(CAOsurv$DFStime), randarm = unique(CAOsurv$randarm)) # ## Calls marginalCDF on each row of nd # ## (done in parallel to speed up computations) # mp <- parallel::mclapply(split(nd, seq(nrow(nd))), # marginalCDF, mod = mod, mc.cores = 1) # mp <- do.call("rbind", mp) # mp$surv <- with(mp, 1 - cdf) # # <> # with(mp[mp$randarm == levs[1], ], lines(iDFS, surv, col = col[1], lwd = lwd)) # with(mp[mp$randarm == levs[2], ], lines(iDFS, surv, col = col[2], lwd = lwd)) # legend("bottomright", legend = levs, col = col, bty = "n", lty = 1, lwd = 1, cex = .8) ## ----MCOX-preproc, echo = FALSE----------------------------------------- ### convert "exact" event dates to interval-censoring (+/- two days) tmp <- CAOsurv$iDFS exact <- tmp[, 3] == 1 tmp[exact, 2] <- tmp[exact, 1] + 2 tmp[exact, 1] <- pmax(tmp[exact, 1] - 2, 0) tmp[exact, 3] <- 3 CAOsurv$iDFS2 <- tmp ## ----MCOX-model-fit, cache = TRUE--------------------------------------- mmc <- mtram(Coxph(iDFS2 ~ randarm, data = CAOsurv, log_first = TRUE), formula = ~ (1 | Block), data = CAOsurv) ## ----MCOX-FUN----------------------------------------------------------- ## marginal HR from "mtram" ## reset seed on.exit mHR.mtram <- function(object, with_confint = FALSE, seed = 1) { stopifnot(inherits(object, "mtram")) cf <- coef(object) cf <- cf[-grep("Bs", names(cf))] stopifnot(length(cf) == 2) mlHR <- cf[1] / sqrt(1 + cf["gamma1"]^2) ret <- mHR <- exp(mlHR) if (with_confint) { set.seed(seed) S <- vcov(object) rbeta <- rmvnorm(10000, mean = coef(object), sigma = S) s <- rbeta[,ncol(rbeta)] rbeta <- rbeta[,-ncol(rbeta)] / sqrt(s^2 + 1) ci <- quantile(exp(rbeta[, ncol(rbeta)]), prob = c(.025, .975)) ret <- c(mHR, ci) ret <- as.array(t(ret)) } return(ret) } ## ----MCOX-summary, cache = TRUE, results = "hide", fig.show = 'hide'---- coef(mmc) sqrt(diag(vcov(mmc))) (ci_MCOX <- mHR.mtram(mmc, with_confint = TRUE)) ## ----HTECOX-model-fit, cache = TRUE------------------------------------- ma <- CoxphME(iDFS ~ randarm + s(age, by = as.ordered(randarm), fx = TRUE, k = 6), data = CAOsurv, log_first = TRUE) nd <- model.frame(ma)[rep(2, 100), ] nd$age <- seq(min(CAOsurv$age), max(CAOsurv$age), length.out = 100) xx <- model.matrix(ma, data = nd, type = "X", keep_sign = FALSE)$X ip <- grep("randarm", names(bb <- coef(ma, with_baseline = TRUE))) vc <- vcov(ma, parm = names(bb)[ip]) bb <- bb[ip] ## NOTE: unadjusted cb <- exp(confint(multcomp::glht(multcomp::parm(bb, vc), linfct = xx), calpha = univariate_calpha())$confint) ## ----HTECOX-summary, cache = TRUE, results = "hide", fig.show = 'hide'---- summary(ma) ## ----HTECOX-HR-plot----------------------------------------------------- ## Plot HR plot(nd$age, cb[, "Estimate"], type = "n", ylab = "Hazard ratio", xlab = "Age (in years)", ylim = ylimHR) polygon(c(nd$age, rev(nd$age)), c(cb[, "lwr"], rev(cb[, "upr"])), border = NA, col = rgb(.1, .1, .1, .1)) lines(nd$age, cb[, "Estimate"], lwd = lwd) abline(h = 1, lty = 3) rug(CAOsurv$age, lwd = 2, col = rgb(.1, .1, .1, .1)) ## ----TRT-load-install, echo = TRUE, eval = FALSE------------------------ # install.packages("trtf") # library("trtf") ## ----TRT-load----------------------------------------------------------- library("trtf") set.seed(4) ## ----TRTF-model-fit, cache = TRUE, warning=TRUE------------------------- tr <- trafotree(Coxph(iDFS ~ randarm, data = CAOsurv, log_first = TRUE), formula = iDFS ~ randarm | age, data = CAOsurv, control = ctree_control(teststat = "maximum", minbucket = 40)) ## ----TRT-results, results = "hide", fig.show = 'hide'------------------- logLik(tr) ## ----TRTF-surv-plot, fig.width = 10, fig.height = 6--------------------- library("ATR") plot(rotate(tr), tp_args = list(newdata = nd1, type = "survivor", col = col, lwd = lwd), terminal_panel = trtf:::node_mlt) ## ----FRAILTY-model-fit, cache = TRUE------------------------------------ mf <- Coxph(iDFS ~ randarm, data = CAOsurv, log_first = TRUE, frailty = "Gamma") ## ----FRAILTY-summary, cache = TRUE, results = "hide", fig.show = 'hide'---- logLik(mf) coef(mf)[trt] coef(mf, addparm = TRUE) confint(mf, parm = c(trt, "logrho")) ## ----LOGIT-model-fit---------------------------------------------------- ml <- Colr(iDFS ~ randarm, data = CAOsurv, log_first = TRUE) ## ----SS, echo = TRUE---------------------------------------------------- m <- as.mlt(Survreg(OS ~ randarm + age, data = CAOsurv, dist = "weibull", support = c(.1, 80 * 365))) ## ----SS-nd-------------------------------------------------------------- N <- 500 nd <- with(CAOsurv, data.frame(randarm = gl(2, N, labels = levels(randarm)), age = rnorm(N, mean = mean(age), sd = sd(age)))) ## ----SS-setup, echo = FALSE--------------------------------------------- `coef<-` <- mlt::`coef<-` ## masked by tramME? ## ----SS-lOR, echo = TRUE------------------------------------------------ cf <- coef(m) cf["randarm5-FU + Oxaliplatin"] <- .25 coef(m) <- cf nd$T <- as.Surv(simulate(m, newdata = nd, K = 1000)) ## ----SS-simC, echo = TRUE----------------------------------------------- cf["(Intercept)"] <- cf["(Intercept)"] + qlogis(.8) coef(m) <- cf nd$C <- as.Surv(simulate(m, newdata = nd, K = 1000)) ## ----SS-OS, echo = TRUE, results = "asis"------------------------------- nd$nOS <- with(nd, Surv(time = pmin(T[, "time"], C[, "time"]), event = T[,"time"] < C[,"time"])) ## ----SS-refit, eval = FALSE, results = "hide", fig.show = "hide"-------- # table(nd$nOS[,"status"]) # levels(nd$randarm)[2] <- "innovative" # plot(survfit(nOS ~ randarm, data = nd), col = col, xlab = "Time (in days)", # ylab = "Probability of survival") # legend("topright", legend = levels(nd$randarm), col = col, lty = 1, bty = "n") # Survreg(nOS ~ randarm + age, data = nd, dist = "weibull")