Merge pull request #1 from raredd/master

minor changes
raredd · Aug 13, 2019 · 9275a45 · 9275a45
2 parents 79519d9 + 20750bb
commit 9275a45
Show file tree

Hide file tree

Showing 8 changed files with 158 additions and 126 deletions.
diff --git a/R/cuminc2.R b/R/cuminc2.R
@@ -1,9 +1,8 @@
 ### formula method for cuminc
-# cuminc2, print.cuminc2, summary.cuminc2, split_cuminc, cuminc_pairs,
-# timepoints2
+# cuminc2, print.cuminc2, summary.cuminc2, cuminc_pairs, timepoints2
 # 
 # unexported:
-# get_events, gy_pval, gy_text, pw_pval, pw_text, name_or_index
+# get_events, gy_pval, gy_text, split_cuminc, pw_pval, pw_text, name_or_index
 ###
 
 
@@ -23,6 +22,10 @@
 #' the censoring indicator will be guessed from \code{formula} but may be
 #' overridden by \code{cencode}
 #' 
+#' @seealso
+#' \code{\link{summary.cuminc2}}; \code{\link{plot.cuminc2}};
+#' \code{\link[cmprsk]{cuminc}}; 
+#' 
 #' @examples
 #' tp <- within(transplant, {
 #'   event_ind <- as.integer(event) - 1L
@@ -56,7 +59,6 @@
 #'   with(tp, cuminc(futime, event_ind, sex, abo))
 #' )
 #' 
-#' 
 #' @export
 
 cuminc2 <- function(formula, data, rho = 0, cencode = NULL,
@@ -208,12 +210,15 @@ summary.cuminc2 <- function(object, times = NULL, digits = 5L, ...) {
   ## "atrisk"
   total_atrisk <- object$cuminc2[order(object$cuminc2$time), ]
   total_atrisk <- get_events(NULL, total_atrisk$time, times, TRUE)
-  # res <- c(total_events) - res
+  atrisk <- c(total_events) - res
+  # atrisk <- rbind(atrisk, Censored = total_atrisk - colSums(atrisk))
 
   l <- list(
     events = res, total_events = total_events,
     total_groups = c(table(object$cuminc2$group)),
-    total_atrisk = total_atrisk
+    atrisk = atrisk, atrisk_sum = colSums(atrisk),
+    total_atrisk = total_atrisk,
+    total_censored = total_atrisk - colSums(atrisk)
   )
 
   c(tp, l)
@@ -266,11 +271,11 @@ gy_text <- function(x, ..., details = TRUE) {
 
   # bquote(paste(chi^2, ' = ', .(txt)))
   if (details)
-    txt
-  else setNames(pvalr(x[, 'pv'], show.p = TRUE), rownames(x))
+    txt else setNames(pvalr(x[, 'pv'], show.p = TRUE), rownames(x))
 }
 
 split_cuminc <- function(x, wh = c('event', 'group'), ws_split = 'last') {
+  ## split (ie, hack) cuminc objects by event or group levels
   c2 <- inherits(x, 'cuminc2')
   sp_str <- switch(
     ws_split,
@@ -295,42 +300,44 @@ split_cuminc <- function(x, wh = c('event', 'group'), ws_split = 'last') {
 
   wh <- match.arg(wh)
 
-  switch(wh,
-         event = {
-           xx_ev <- lapply(ev, function(ii) {
-             tmp <- list(cuminc = xx[grep(ii, names(xx), fixed = TRUE)])
-
-             if (gy)
-               tmp$cuminc$Tests <- ci$Tests[grep(ii, rownames(ci$Tests)), , drop = FALSE]
-
-             if (c2) {
-               cc  <- as.character(x$cuminc2$cencode[1L])
-               dat <- x$cuminc2[x$cuminc2[, 'status'] %in% c(ii, cc), ]
-               tmp <- c(tmp, list(cuminc2 = droplevels(dat)))
-             }
-
-             if (c2)
-               structure(tmp, class = 'cuminc2') else tmp
-           })
-
-           setNames(xx_ev, ev)
-         },
-
-         group = {
-           xx_gr <- lapply(gr, function(ii) {
-             tmp <- list(cuminc = xx[grep(ii, names(xx), fixed = TRUE)])
-
-             if (c2) {
-               dat <- x$cuminc2[x$cuminc2[, 'group'] %in% ii, ]
-               tmp <- c(tmp, list(cuminc2 = droplevels(dat)))
-             }
-
-             if (c2)
-               structure(tmp, class = 'cuminc2') else tmp
-           })
-
-           setNames(xx_gr, gr)
-         }
+  switch(
+    wh,
+    event = {
+
+      xx_ev <- lapply(ev, function(ii) {
+        tmp <- list(cuminc = xx[grep(ii, names(xx), fixed = TRUE)])
+
+        if (gy)
+          tmp$cuminc$Tests <- ci$Tests[grep(ii, rownames(ci$Tests)), , drop = FALSE]
+
+        if (c2) {
+          cc  <- as.character(x$cuminc2$cencode[1L])
+          dat <- x$cuminc2[x$cuminc2[, 'status'] %in% c(ii, cc), ]
+          tmp <- c(tmp, list(cuminc2 = droplevels(dat)))
+        }
+
+        if (c2)
+          structure(tmp, class = 'cuminc2') else tmp
+      })
+
+      setNames(xx_ev, ev)
+    },
+
+    group = {
+      xx_gr <- lapply(gr, function(ii) {
+        tmp <- list(cuminc = xx[grep(ii, names(xx), fixed = TRUE)])
+
+        if (c2) {
+          dat <- x$cuminc2[x$cuminc2[, 'group'] %in% ii, ]
+          tmp <- c(tmp, list(cuminc2 = droplevels(dat)))
+        }
+
+        if (c2)
+          structure(tmp, class = 'cuminc2') else tmp
+      })
+
+      setNames(xx_gr, gr)
+    }
   )
 }
 
@@ -382,6 +389,7 @@ cuminc_pairs <- function(object, data = NULL, rho = 0, cencode = NULL,
     inherits(object, c('cuminc2', 'formula'))
   )
 
+  ## pairwise gray tests
   pwgray <- function(i, j, k) {
     force(k)
     data <- data[data[, 'group'] %in% c(unq[i], unq[j]), ]
@@ -500,10 +508,10 @@ pw_text <- function(formula, data, ..., details = TRUE, pFUN = NULL,
 
 name_or_index <- function(x, y = NULL) {
   ## return integer vector where x occurs in y
-  ## cmprsk2:::name_or_index(c('1', '3', 'e'))
-  ## cmprsk2:::name_or_index(c('a', 'c', 'e'), letters)
-  ## table is given priority over integer, eg, idx = 27 instead of 4
-  ## cmprsk2:::name_or_index(c('a', '4', 'e'), c(letters, '4'))
+  # cmprsk2:::name_or_index(c('1', '3', 'e'))
+  # cmprsk2:::name_or_index(c('a', 'c', 'e'), letters)
+  # table is given priority over integer, eg, idx = 27 instead of 4
+  # cmprsk2:::name_or_index(c('a', '4', 'e'), c(letters, '4'))
   suppressWarnings(
     ix <- as.integer(x)
   )

diff --git a/R/plot.R b/R/plot.R
@@ -19,14 +19,17 @@
 #' @param wh.events a character string giving the type of \code{events} table
 #' to show; one of \code{"events"} (cumulative number of events), \code{"est"}
 #' (estimates, see \code{\link[cmprsk]{timepoints}}), \code{"est.sd"} (estimate
-#' +/- standard deviation), or \code{"est.ci"} (estimate with confidence
-#' inverval)
+#' +/- standard deviation), \code{"est.ci"} (estimate with confidence
+#' inverval), or \code{"atrisk"} (event at-risk table with censored)
 #' @param events.lab heading for events table
 #' @param events.digits when estimates are shown in events table (see
 #' \code{wh.events}), number of digits past the decimal to show
 #' @param events.lines logical; draw lines next to groups in events table
 #' @param events.col logical or a vector with colors for events table text;
 #' if \code{TRUE}, \code{col.ci} will be used
+#' @param include_censored logical; if \code{TRUE}, censored patients are
+#' included in at-risk counts (may not be desired); the default (\code{FALSE})
+#' will only sum the at-risk table rows
 #' @param main title of plot
 #' @param xlab,ylab x- and y-axis labels
 #' @param groups.lab labels for each line in \code{events} table
@@ -60,7 +63,8 @@
 #' \code{tcl}, \code{cex.lab}, \code{xaxs}, etc) passed to \code{\link{par}}
 #' 
 #' @seealso
-#' \code{\link{cuminc2}}; \code{\link{summary.cuminc2}}
+#' \code{\link{ciplot_by}}; \code{\link{cuminc2}};
+#' \code{\link{summary.cuminc2}}; \code{\link[cmprsk]{cuminc}}
 #' 
 #' @examples
 #' tp <- within(transplant, {
@@ -92,10 +96,11 @@ ciplot <- function(x,
                    lty.ci = par('lty'), lwd.ci = par('lwd'),
 
                    events = TRUE, atrisk = TRUE, events.total = TRUE,
-                   wh.events = c('events', 'est', 'est.sd', 'est.ci'),
+                   wh.events = c('events', 'est', 'est.sd', 'est.ci', 'atrisk'),
                    events.lab = NULL,
                    events.digits = 3L,
                    events.lines = TRUE, events.col = FALSE,
+                   include_censored = FALSE,
 
                    main = NULL,
                    xlab = 'Time', ylab = 'Probability',
@@ -249,7 +254,8 @@ ciplot <- function(x,
           events = 'Cumulative events',
           est    = 'Estimate',
           est.sd = 'Estimate +/- Std. dev',
-          est.ci = 'Estimate [LCI, UCI]'
+          est.ci = 'Estimate [LCI, UCI]',
+          atrisk = 'Number at risk'
         )
       else events.lab
 
@@ -299,7 +305,8 @@ ciplot <- function(x,
       est.sd = timepoints2(x, times = events.at, sd = TRUE, html = FALSE,
                            ci = FALSE, digits = events.digits),
       est.ci = timepoints2(x, times = events.at, sd = FALSE, html = FALSE,
-                           ci = TRUE, digits = events.digits)
+                           ci = TRUE, digits = events.digits),
+      atrisk = summary(x, times = events.at)$atrisk
     )
 
     d1 <- data.frame(time = rep(as.numeric(colnames(ss)), each = nrow(ss)),
@@ -312,7 +319,9 @@ ciplot <- function(x,
     if (atrisk)
       d2 <- c(d2, list('At-risk' = data.frame(
         time = as.numeric(colnames(ss)),
-        n.risk = summary(x, times = events.at)$total_atrisk,
+        n.risk = if (include_censored)
+          summary(x, times = events.at)$total_atrisk
+        else summary(x, times = events.at)$atrisk_sum,
         strata = 'At-risk',
         stringsAsFactors = FALSE
       )))

diff --git a/R/rpart.R b/R/rpart.R
diff --git a/R/stat.R b/R/stat.R
@@ -221,31 +221,33 @@ terms.crr <- terms.crr2
 #' 
 #' ## model selection
 #' m2 <- with(bmt, crr(ftime, Status, cov1[, c(4:6)]))
-#' m3 <- with(bmt, crr(ftime, Status, cov1[, c(4:6,7)]))
-#' m4 <- with(bmt, crr(ftime, Status, cov1[, c(4:6,7,1)]))
-#' m5 <- with(bmt, crr(ftime, Status, cov1[, c(4:6,7,2)]))
-#' m6 <- with(bmt, crr(ftime, Status, cov1[, c(4:6,7,3)]))
+#' m3 <- with(bmt, crr(ftime, Status, cov1[, c(4:6, 7)]))
+#' m4 <- with(bmt, crr(ftime, Status, cov1[, c(4:6, 7, 1)]))
+#' m5 <- with(bmt, crr(ftime, Status, cov1[, c(4:6, 7, 2)]))
+#' m6 <- with(bmt, crr(ftime, Status, cov1[, c(4:6, 7, 3)]))
 #' 
 #' crrFits(m1, m2, m3, m4, m5, m6, p = 3)
 #' 
 #' 
-#' par(mfrow = c(2,2))
+#' par(mfrow = c(2, 2))
 #' with(m2, {
 #'   for (ii in seq.int(ncol(res)))
-#'     scatter.smooth(uftime, res[, ii],
-#'                    main = gsub('Phase', '', names(coef)[ii]),
-#'                    xlab = 'Failure time',
-#'                    ylab = 'Schoenfeld residuals')
+#'     scatter.smooth(uftime, res[, ii], main = names(coef)[ii],
+#'                    xlab = 'Failure time', ylab = 'Schoenfeld residuals')
 #' })
 #' 
 #' pred <- with(bmt, model.matrix(~ levels(Phase)))[, -1L]
 #' pred <- predict(m2, pred)
 #' 
-#' plot(pred, xlab = 'Failure time', ylab = 'CIF', col = 1:4, lty = 1,
-#'      ylim = c(0, 1))
-#' legend('top', lty = 1L, col = 1:4, horiz = TRUE, bty = 'n',
-#'        title = 'Phase', legend = levels(bmt$Phase),
-#'        x.intersp = 0.1, y.intersp = 0.5)
+#' plot(
+#'   pred, col = 1:4, lty = 1, ylim = c(0, 1),
+#'   xlab = 'Failure time', ylab = 'CIF'
+#' )
+#' legend(
+#'   'top', lty = 1L, col = 1:4, horiz = TRUE, bty = 'n',
+#'   title = 'Phase', legend = levels(bmt$Phase),
+#'   x.intersp = 0.1, y.intersp = 0.5
+#' )
 #' }
 #' 
 #' @export