In ampel-leipzig/glmnettools: Additional Tools for glmnet

Introduction

The glmnettools package offers functions for repeated cross-validation and prediction of survival probabilities using the glmnet package [@friedman2010; @simon2011].

Repeated Cross-Validation

The example is taken from the cv.glmnet manual page (?"glmnet::cv.glmnet"):

library("glmnettools")

# Generating example data
set.seed(1010)
n <- 1000
p <- 100
nzc <- trunc(p/10)
x <- matrix(rnorm(n * p), n, p)
beta <- rnorm(nzc)
fx <- x[, seq(nzc)] %*% beta
eps <- rnorm(n) * 5
y <- drop(fx + eps)

set.seed(1011)
# nrepcv should usually be higher but to keep the runtime of the example low
# we choose 2 here, same for the nfolds
rcvob <- rcv.glmnet(x, y, nrepcv = 2, nfolds = 3)
plot(rcvob)
title("Gaussian Family", line = 2.5)
head(coef(rcvob))
predict(rcvob, newx = x[1:5, ], s = "lambda.min")

Prediction Survival Probabilities

The example data are taken from the survival package [@survival].

library("survival")

data("pbc")
pbc <- pbc[complete.cases(pbc),]
x <- data.matrix(pbc[!colnames(pbc) %in% c("id", "time", "status")])
y <- Surv(pbc$time, as.integer(pbc$status == 2))

rcvob <- rcv.glmnet(x, y, family = "cox", nrepcv = 2, nfolds = 3)

plot(rcvob)
plot(rcvob$glmnet.fit, label = TRUE)

nz <- predict(rcvob, type = "nonzero", s = "lambda.1se")[, 1]
nz
colnames(x)[nz]

## predicting survival probabilities
predict(
    rcvob,
    newx = x[1:5,], x = x, y = y, s = "lambda.1se",
    type = "survival", times = c(2, 5) * 365.25
)

## linear predictor
predict(rcvob, newx = x[1:5,], type = "link", s = "lambda.1se")

Convert `rcv.glmnet`/`coxnet` Objects into `coxph` or `cph` Objects

For further investigation/processing packages like survival and/or rms offer a lot of useful functions [@survival; @rms]. That's why glmnettools provides some helper functions to convert between these objects.

library("rms")
tm <- 2 * 365.2
cx <- as.cph(rcvob, x = x, y = y, s = "lambda.1se", time.inc = tm)

set.seed(123)
## B should be much higher in real life
cb <- calibrate(cx, B = 10, u = tm)
cb
plot(cb)