Modifications in Cox model visualization

R/ggcoxfunctional.R

@@ -3,18 +3,19 @@
 #' @importFrom stats approx
 #' @importFrom stats resid
 #' @importFrom survival coxph
-#' @importFrom magrittr %>%
 NULL
 #' Functional Form of Continuous Variable in Cox Proportional Hazards Model
 #'@description Displays graphs of continuous explanatory variable against martingale residuals of null
 #'cox proportional hazards model, for each term in of the right side of \code{formula}. This might help to properly
 #'choose the functional form of continuous variable in cox model (\link{coxph}). Fitted lines with \link{lowess} function
 #'should be linear to satisfy cox proportional hazards model assumptions.
+#'
 #'@param fit an object of class \link{coxph.object} - created with \link{coxph} function.
-#'@param formula a formula object, with the response on the left of a ~ operator, and the terms on the right. The response must be a survival object as returned by the \link{Surv} function.
+#'@param formula (deprecated) a formula object, with the response on the left of a ~ operator, and the terms on the right. The response must be a survival object as returned by the \link{Surv} function.
 #'@param data a \code{data.frame} in which to interpret the variables named in the formula,
 #'@param iter parameter of \link{lowess}.
 #'@param f parameter of \link{lowess}.
+#'@param vars.keep character vector of variables for which we want to draw the plot.
 #'@param xlim,ylim x and y axis limits e.g. xlim = c(0, 1000), ylim = c(0, 1).
 #'@param ylab y axis label.
 #'@param title the title of the final \link{grob} (\code{top} in \link{arrangeGrob})
@@ -31,16 +32,17 @@ NULL
 #'
 #' library(survival)
 #' data(mgus)
-#' res.cox <- coxph(Surv(futime, death) ~ mspike + log(mspike) + I(mspike^2) +
+#' res.cox <- coxph(Surv(futime, death) ~ mspike + sex + log(mspike) + I(mspike^2) +
 #'     age + I(log(age)^2) + I(sqrt(age)), data = mgus)
-#' ggcoxfunctional(res.cox,  data = mgus, point.col = "blue", point.alpha = 0.5)
+#' ggcoxfunctional(res.cox, data = mgus, point.col = "blue", point.alpha = 0.5)
+#' ggcoxfunctional(res.cox, data = mgus, vars.keep=c("mspike", "log(mspike)"))
 #' ggcoxfunctional(res.cox, data = mgus, point.col = "blue", point.alpha = 0.5,
 #'                 title = "Pass the title", caption = "Pass the caption")
 #'
 #'
 #'@describeIn ggcoxfunctional Functional Form of Continuous Variable in Cox Proportional Hazards Model.
 #'@export
-ggcoxfunctional <- function (formula, data = NULL, fit, iter = 0, f = 0.6,
+ggcoxfunctional <- function (formula, data = NULL, fit, iter = 0, f = 0.6, vars.keep,
                              point.col = "red", point.size = 1, point.shape = 19, point.alpha = 1,
                              xlim = NULL, ylim = NULL,
                              ylab = "Martingale Residuals \nof Null Cox Model",
@@ -58,9 +60,20 @@ ggcoxfunctional <- function (formula, data = NULL, fit, iter = 0, f = 0.6,
   }
   formula <- fit$formula
   data <- .get_data(fit, data)
+  remov = sapply(attr(stats::terms(formula), "term.labels"),
+                 function(x) {!is.numeric(with(data, eval(parse(text=x))))})
+  if(any(remov))
+    formula <- stats::drop.terms(stats::terms(formula), which(remov), keep.response=TRUE)
+
+  explanatory.variables.names <- attr(stats::terms(formula), "term.labels")
+
+  if(!missing(vars.keep)){
+    if(!all(vars.keep %in% explanatory.variables.names)) stop("Unknown or non-numeric `vars.keep`")
+    explanatory.variables.names <- vars.keep
+  }
 
-  attr(stats::terms(formula), "term.labels") -> explanatory.variables.names
-  stats::model.matrix(formula, data = data) -> explanatory.variables.values
+  explanatory.variables.values <- stats::model.matrix(formula, data = data)
+  data = data[rownames(explanatory.variables.values), ]
   SurvFormula <- deparse(formula[[2]])
   martingale_resid <- lowess_x <- lowess_y <- NULL
   lapply(explanatory.variables.names, function(i){

R/ggcoxzph.R

@@ -1,7 +1,7 @@
 #'Graphical Test of Proportional Hazards with ggplot2
 #'@description Displays a graph of the scaled Schoenfeld residuals, along with a
-#'  smooth curve using \pkg{ggplot2}. Wrapper around \link{plot.cox.zph}.
-#'@param fit an object of class \link{cox.zph}.
+#'  smooth curve using \pkg{ggplot2}. Wrapper around \link{plot.cox.zph}. If fit is a \code{coxph}, the function will also plot the coefficient estimated by the model as a horizontal line.
+#'@param fit an object of class \link{coxph} or \link{cox.zph}.
 #'@param resid	a logical value, if TRUE the residuals are included on the plot,
 #'  as well as the smooth fit.
 #'@param se a logical value, if TRUE, confidence bands at two standard errors
@@ -11,8 +11,11 @@
 #'@param nsmo	number of points used to plot the fitted spline.
 #'@param var the set of variables for which plots are desired. By default, plots
 #'  are produced in turn for each variable of a model.
+#'@param var_pval the max Grambsch-Therneau test pvalue for which plots are desired. Use only one of \code{var} or \code{var_pval}.
 #'@param point.col,point.size,point.shape,point.alpha color, size, shape and visibility to be used for points.
 #'@param caption the caption of the final \link{grob} (\code{bottom} in \link{arrangeGrob})
+#'@param zph.transform the argument to pass to \code{survival::cox.zph} if \code{fit} is a \code{coxph} object
+#'@param hline_size the argument to pass to \code{survival::cox.zph} if \code{fit} is a \code{coxph} object
 #'@param ggtheme function, ggplot2 theme name.
 #'  Allowed values include ggplot2 official themes: see \code{\link[ggplot2]{theme}}.
 #'@param ... further arguments passed to either the print() function or to the \code{\link[ggpubr]{ggpar}} function for customizing the plot (see Details section).
@@ -37,22 +40,43 @@
 #' cox.zph.fit <- cox.zph(fit)
 #' # plot all variables
 #' ggcoxzph(cox.zph.fit)
+#' ggcoxzph(fit)
+#' # plot all variables for which the Grambsch-Therneau test pvalue is less than 0.55
+#' ggcoxzph(fit, var_pval=0.55)
 #' # plot all variables in specified order
 #' ggcoxzph(cox.zph.fit, var = c("ecog.ps", "rx", "age"), font.main = 12)
 #' # plot specified variables in specified order
 #' ggcoxzph(cox.zph.fit, var = c("ecog.ps", "rx"), font.main = 12, caption = "Caption goes here")
 #'
 #'@describeIn ggcoxzph Graphical Test of Proportional Hazards using ggplot2.
 #'@export
-ggcoxzph <- function (fit, resid = TRUE, se = TRUE, df = 4, nsmo = 40, var,
+ggcoxzph <- function (fit, resid = TRUE, se = TRUE, df = 4, nsmo = 40, var, var_pval, 
                       point.col = "red", point.size = 1, point.shape = 19, point.alpha = 1,
-                      caption = NULL,
-                      ggtheme = theme_survminer(), ...){
-
-  x <- fit
-  if(!methods::is(x, "cox.zph"))
-    stop("Can't handle an object of class ", class(x))
-
+                      caption = NULL, zph.transform="km", hline_size = 1.25,
+                      ggtheme = theme_survminer(), ...){ 
+
+
+  if(!methods::is(fit, "cox.zph") && !methods::is(fit, "coxph"))
+    stop("Can't handle an object of class ", class(fit))
+
+  if(methods::is(fit, "coxph")){
+    COX_FIT <- TRUE
+    x <- survival::cox.zph(fit, transform = zph.transform)
+  } else{
+    COX_FIT <- FALSE
+    x <- fit
+  }
+
+  if(!missing(var_pval)){
+    if(!is.numeric(var_pval))
+      stop("var_pval should be a numeric vector")
+    if(!missing(var))
+      stop("Can't handle both var and var_pval")
+    tmp <- x$table[,"p"]
+    var <- names(tmp[tmp<var_pval])
+    var <- var[var!="GLOBAL"]
+  }
+
   xx <- x$x
   yy <- x$y
   d <- nrow(yy)
@@ -146,15 +170,33 @@ ggcoxzph <- function (fit, resid = TRUE, se = TRUE, df = 4, nsmo = 40, var,
   names(plots) <- var
   class(plots) <- c("ggcoxzph", "ggsurv", "list")
 
+  if(COX_FIT){
+    if(missing(var)||is.null(var)||is.na(var)){
+      betas <- fit$coefficients
+    } else {
+      betas <- fit$coefficients[var]
+    }
+
+    if (length(betas) != length(plots))
+      stop("Error: length of plot list (", length(plots),
+           ") is different from length of coefficients (", length(betas), ")")
+    for (i in 1:length(betas)) {
+      p <- plots[[i]]
+      p <- p + geom_hline(yintercept = betas[i], color = "steel blue",
+                          size = hline_size)
+      plots[[i]] <- p
+    }
+  }
+
   if("GLOBAL" %in% rownames(x$table)) # case of multivariate Cox
     global_p <- x$table["GLOBAL", 3]
   else global_p <- NULL # Univariate Cox
   attr(plots, "global_pval") <- global_p
   attr(plots, "caption") <- caption
   plots
-
 }
 
+
 #' @param x an object of class ggcoxzph
 #' @param newpage open a new page. See \code{\link{grid.arrange}}.
 #' @method print ggcoxzph

fork_survminer.Rproj

@@ -0,0 +1,18 @@
+Version: 1.0
+
+RestoreWorkspace: Default
+SaveWorkspace: Default
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 2
+Encoding: UTF-8
+
+RnwWeave: Sweave
+LaTeX: pdfLaTeX
+
+BuildType: Package
+PackageUseDevtools: Yes
+PackageInstallArgs: --no-multiarch --with-keep.source
+PackageRoxygenize: rd,collate,namespace