feat: add plot.rcvglmnet(..., what = "path")

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: ameld
 Title: Data and Model of End-Stage Liver Disease used in the AMPEL Project
-Version: 0.0.24
+Version: 0.0.25
 Description:
     A dataset of patients evaluated for liver transplantation at the
     University Hospital Leipzig from November 2012 to June 2015.

NAMESPACE

@@ -47,6 +47,7 @@ importFrom(graphics,plot.new)
 importFrom(graphics,plot.window)
 importFrom(graphics,points)
 importFrom(graphics,segments)
+importFrom(graphics,strheight)
 importFrom(graphics,strwidth)
 importFrom(graphics,text)
 importFrom(graphics,title)

NEWS.md

@@ -2,6 +2,11 @@
 
 ## Changes in development
 
+## Changes in 0.0.25
+
+- Add `plot.rcvglmnet(..., what = "path")` to plot lambda path
+  (extends `glmnet::plot.glmnet`).
+
 ## Changes in 0.0.24
 
 - Pass `main` argument to `.plot.cal`, affected function: `plot.boot.glmnet`.

R/coxnet.R

@@ -18,6 +18,7 @@ basehaz <- function(fit, ...)UseMethod("basehaz")
 #' @param centered `logical(1)`, see [`survival::basehaz()`].
 #' @rdname basehaz
 #' @importFrom survival basehaz
+#' @aliases basehaz.coxph
 #' @method basehaz coxph
 #' @export
 basehaz.default <- basehaz.coxph <- function(fit, centered = TRUE, ...) {

R/rcv.glmnet.R

@@ -137,28 +137,78 @@ rcv.glmnet <- function(x, y, lambda = NULL, alpha = 1,
     out
 }
 
-#' Plot the cross-validation curve
+#' Plot the cross-validation curve/lambda path
 #'
 #' This functions plots the aggregated cross-validation curve produced by
-#' [`rcv.glmnet()`].
+#' [`rcv.glmnet()`] or the lambda path for the averaged glmnet.fit object.
 #'
-#' @param x `rcv.glmnet` object.
-#' @param \dots further arguments passed to `plot.cv.glmnet`.
+#' @param x `rcv.glmnet`, object.
+#' @param what `character`, uses [`glmnet::plot.cv.glmnet()`] for "cv" and
+#' [`glmnet::plot.glmnet()`] for "path".
+#' @param \dots further arguments passed to `plot.cv.glmnet` or `plot.glmnet`.
+#'
+#' @details
+#' For `what = "path"` the original `plot.glmnet` is extended by labelling the
+#' top `nlabel` (default: 9) labels on the right side of the plot.
 #'
 #' @author Sebastian Gibb
-#' @seealso [`glmnet::cv.glmnet()`]
-#' @importFrom graphics title
+#' @seealso [`glmnet::cv.glmnet()`], [`glmnet::plot.cv.glmnet()`],
+#' [`glmnet::plot.glmnet()`]
+#' @importFrom graphics title strheight
 #' @method plot rcv.glmnet
 #' @export
-plot.rcv.glmnet <- function(x, ...) {
-    NextMethod()
-    title(
-        sub = paste(
-            "Averaged across", x$nrepcv, "repeated cross-validations",
-            "each with", x$nfolds, "folds."
-        ),
-        adj = 0L
-    )
+plot.rcv.glmnet <- function(x, what = c("cv", "path"), ...) {
+
+    what <- match.arg(what)
+    if (what == "path") {
+        .plot.glmnet(x, ...)
+    } else {
+        NextMethod()
+        title(
+            sub = paste(
+                "Averaged across", x$nrepcv, "repeated cross-validations",
+                "each with", x$nfolds, "folds."
+            ),
+            adj = 0L
+        )
+    }
+}
+
+.plot.glmnet <- function(x, nlabel = 9, cex.lab = 1,
+                         col = viridisLite::cividis(nlabel),
+                         ...) {
+    beta <- x$glmnet.fit$beta
+    nr <- nrow(beta)
+    col <- rep_len(col, nrow(beta))
+    beta <- beta[.nonzero(beta), ncol(beta)]
+    o <- order(-abs(beta))
+    o <- o[seq_len(min(c(nlabel, length(o))))]
+    beta <- beta[o]
+
+    old.par <- par(no.readonly = TRUE)
+    on.exit(par(old.par))
+
+    mai <- par("mai")
+    w <- max(strwidth(names(beta), "inch") * cex.lab, na.rm = TRUE) + 1/8
+    if (mai[4L] < w)
+        mai[4L] <- mai[4L] + w # taken from dotchart
+    old.par <- par(mai = mai, no.readonly = TRUE)
+
+    plot(x$glmnet.fit, col = col, label = FALSE, ...)
+
+    abline(h = 0, lty = "dotted", col = "#808080")
+
+    beta <- .avoid_ylab_overlap(beta, strheight("X") * cex.lab)
+
+    for (i in seq_along(beta))
+        axis(
+            4L,
+            at = beta[i], labels = names(beta)[i],
+            las = 1,
+            cex.axis = cex.lab,
+            col.axis = col[o[i]],
+            col = col[o[i]]
+        )
 }
 
 #' Predictions for a `rcv.glmnet` object

R/utils-glmnet.R

@@ -1,3 +1,36 @@
+#' Avoid overlap in plotting ylabs
+#'
+#' Function to move ylabs up and down to avoid overlap of labels in lambda path
+#' for plot.rcv.glmnet(..., method = "path").
+#'
+#' @param y `numeric()`, y values.
+#' @param h `numeric(1)`, height of a single entry, `strheight` would be
+#' useful.
+#' @return `numeric()`, modified y.
+#' @noRd
+.avoid_ylab_overlap <- function(y, h) {
+    o <- order(y)
+    ys <- y[o]
+    n <- length(ys)
+    n1 <- n + 1L
+
+    step <- (h / 4)
+
+    d <- ys[-1L] - ys[-n]
+    iter <- 0L
+
+    while (any(d < h) && iter < 10L) {
+        w <- c(FALSE, d < h, FALSE)
+        ys[w[-1L]] <- ys[w[-1L]] - step
+        ys[w[-n1]] <- ys[w[-n1]] + step
+        d <- ys[-1L] - ys[-n]
+        iter <- iter + 1L
+    }
+
+    y[o] <- ys
+    y
+}
+
 #' Create balanced CV folds
 #'
 #' @param y `factor`, classes
@@ -43,6 +76,15 @@
     do.call(rbind, lapply(integer(nrep), function(i)f(y, nfolds = nfolds)))
 }
 
+#' Nonzero Coefficients
+#'
+#' @param x `matrix`.
+#' @return `integer` row indices with nonzero beta coefficients from rcv.glmnet.
+#' @noRd
+.nonzero <- function(x) {
+    which(as.vector(x %*% rep(1L, ncol(x))) != 0)
+}
+
 #' Convert s into its numeric equivalent
 #'
 #' This function converts s/lambda to its numeric equivalent.