#' The glmnet model postprocessor
#'
#' The function is based on the results returned by \code{\link[broom]{glmnet_tidiers}},
#' with an optional application of results filtering using the \code{\link{structural.changes.col.matrix}} function.
#'
#' @param glmnet.res the result of the \code{\link[glmnet]{glmnet}} function
#' @param nonzero.coeffs.only if \code{TRUE} only the non-zero model coefficients are returned in the \code{coeffs} field
#' @param structural.filter.type the type of the \code{steps} filter, selecting the steps from the \emph{lambda path} resulting in the change of the model structure:
#' \itemize{
#' \item{(\code{"none"}): all steps are returned (default)}
#' \item{(\code{"initial"}): return initial step of each subpath consisting of the same variables}
#' \item{(\code{"final"}): return final step of each subpath consisting of the same variables}
#' }
#' @param pseudoAIC \code{(logical)} should AIC be calculated for each model;
#' Remark: the value being calculated is \code{2*df -2*loglike + 2*loglike(NULL)}, where \code{df} is an approximation of df returned by glmnet
#'
#' @return the \code{list} with following fields:
#' \itemize{
#' \item{\code{steps}: a \code{data.frame} with columns:
#' \itemize{
#' \item{\code{step}}
#' \item{\code{lambda}}
#' \item{\code{dev.ratio}}
#' \item{\code{df}}
#' \item{\code{dev}}
#' \item{(optionally) \code{pseudoAIC}}
#' }
#' }
#' \item{\code{coeffs}: a \code{data.frame} with columns:
#' \itemize{
#' \item{\code{step}}
#' \item{\code{term}}
#' \item{\code{estimate}}
#' }
#' }
#' }
#' @export
glmnet.process.output <- function(glmnet.res, nonzero.coeffs.only = FALSE, structural.filter.type = c("none", "initial", "final"), pseudoAIC = FALSE){
glmnet.out.prc <- if (inherits(glmnet.res, "multnet") || !nonzero.coeffs.only) {
glmnet.process.output.reference(glmnet.res, nonzero.coeffs.only)
} else {
glmnet.process.output.nonzero(glmnet.res)
}
glmnet.out.steps <- glmnet.out.prc$steps
glmnet.out.steps$df <- glmnet.res$df
glmnet.out.steps$dev <- glmnet::deviance.glmnet(glmnet.res)
if (pseudoAIC) {
null.deviance <- glmnet.res$nulldev
# dev == 2*(loglike_sat - loglike)
# null.deviance = 2*(loglike_sat - loglike(NULL model)
# a "NULL model", is a model with icept only (or "empty model" - coxph)
# AIC = 2*no.params - 2*ln(Like)
# dev - null.deviance ==
# == 2*(loglike_sat - loglike) - 2*(loglike_sat - loglike(NULL) ==
# == -2*loglike + 2*loglike(NULL)
# pseudoAIC = 2*df + (dev - null.deviance) ==
# == 2*no.params -2*loglike + 2*loglike(NULL) == AIC + const
glmnet.out.steps <- glmnet.out.steps %>%
dplyr::mutate(
pseudoAIC = 2*df + (dev - null.deviance)
)
}
glmnet.out.coeffs <- glmnet.out.prc$coeffs
if (structural.filter.type[1] != "none") {
changes.idx <- structural.changes.col.matrix(
mtx = stats::coef(glmnet.res),
initial.substeps = ("initial" %in% structural.filter.type),
final.substeps = ("final" %in% structural.filter.type)
)
glmnet.out.steps <- glmnet.out.steps %>% dplyr::filter(step %in% changes.idx)
glmnet.out.coeffs <- glmnet.out.coeffs %>% dplyr::filter(step %in% changes.idx)
}
list(
steps = glmnet.out.steps,
coeffs = glmnet.out.coeffs
)
}
glmnet.process.output.reference <- function(glmnet.res, nonzero.coeffs.only){
glmnet.out <- broom::tidy(glmnet.res)
glmnet.out.steps <- glmnet.out %>%
dplyr::select(step, lambda, dev.ratio) %>%
dplyr::distinct() %>%
dplyr::arrange(step)
if (inherits(glmnet.res, "multnet")) {
glmnet.out.coeffs <- glmnet.out %>%
dplyr::select(step, class, term, estimate)
} else {
glmnet.out.coeffs <- glmnet.out %>%
dplyr::select(step, term, estimate)
}
if (nonzero.coeffs.only) {
glmnet.out.coeffs <- glmnet.out.coeffs %>%
dplyr::filter(estimate != 0)
}
list(
steps = glmnet.out.steps,
coeffs = glmnet.out.coeffs
)
}
glmnet.process.output.nonzero <- function(glmnet.res) {
beta <- stats::coef(glmnet.res)
# coeff names
beta.coeffs <- rownames(beta)
# summary for non-empty elems
beta.s <- Matrix::summary(beta)
# coefficients descriptor
model.coeffs <- data.frame(
step = beta.s$j,
term = beta.coeffs[beta.s$i],
estimate = beta.s$x,
stringsAsFactors = FALSE
)
model.steps <- data.frame(
step = 1:glmnet.res$dim[2],
lambda = glmnet.res$lambda,
dev.ratio = glmnet.res$dev.ratio
)
list(
steps = model.steps,
coeffs = model.coeffs
)
}
#' Extraction of critical steps from Lasso path
#'
#' Function extracts the steps of the LASSO path resulting in the structural changes of
#' the model coefficients matrix (adding/dropping coefficient)
#' as well as the preceeding ones
#'
#' @param mtx the coefficients
#' @param strict if \code{TRUE} only the points of the structural change will be returned
#' @return \code{(numeric)} vector containing the steps with structural changes
#' @export
structural.changes.col.matrix <- function(mtx, initial.substeps = TRUE, final.substeps = TRUE) {
if (ncol(mtx) == 1) {
return(c(1))
}
# getting the inclusion/exclusion matrix
coeff.incl <- (mtx != 0)
# first coeff set - s0 - is always novel, we compare starting with s1
coeff.incl.2 <- coeff.incl[, -1]
# dimensions must match, we do remove the last column
coeff.incl.3 <- coeff.incl[, -ncol(coeff.incl)]
coeff.incl.diff <- (coeff.incl.2 != coeff.incl.3)
coeff.incl.diff.reduced <- apply(coeff.incl.diff, 2, any)
# coeff.incl.diff.reduced[i] == FALSE <-> coeffs[i+1] == coeffs[i]
# coeff.incl.diff.reduced[i] == TRUE <-> coeffs[i+1] != coeffs[i]
# that is, coeffs[i+1] are DIFFERENT (***)
coeff.incl.diff.prv <- which(coeff.incl.diff.reduced)
# all steps with novelties
# adding the s0, adding (+1) to account for the (***)
subpaths.initial.elems <- c(1, 1 + coeff.incl.diff.prv)
subpaths.final.elems <- c(coeff.incl.diff.prv, ncol(coeff.incl))
res <- c()
if (initial.substeps) {
res <- subpaths.initial.elems
}
if (final.substeps) {
res <- c(res, subpaths.final.elems)
}
if (initial.substeps && final.substeps) {
# removing duplicates, which may occur
res <- sort(unique(res))
}
res
}
#' Extraction of critical points from Lasso path
#'
#' Function extracts the points of the structure change (adding/dropping coefficient)
#' as well as the preceeding ones, along with their deviance
#'
#' @param model.full.lasso the model fitted by \code{\link[glmnet]{glmnet}} function
#' @param strict if \code{TRUE} only the points of the structural change will be returned
#' @export
inference.lasso.path.structural.changes <- function(model.full.lasso, strict = FALSE){
# extracting pts of structural change of the model
# by comparing inclusion/exclusion matrices
coeff.mtx <- stats::coef(model.full.lasso)
changes.idx <- structural.changes.col.matrix(
mtx = coeff.mtx,
initial.substeps = TRUE,
final.substeps = strict
)
lasso.struct.coeffs <- coeff.mtx[, changes.idx]
Deviance <- stats::deviance(model.full.lasso)[changes.idx]
lasso.struct <- rbind(lasso.struct.coeffs, Deviance)
lasso.struct
}
#' @export
glmnet.predict.models <- function(glmnet.res, x, y, s, offset = NULL, exact = FALSE, ...) {
# since the update() of the glmnet not always works, let's expect it won't and do the job ;
# from the glmnet package:
if (exact && (!is.null(s))) {
which.s.computed <- match(s, glmnet.res$lambda, FALSE)
if (!all(which.s.computed > 0)) {
lambda <- unique(rev(sort(c(s, glmnet.res$lambda))))
glmnet.res <- glmnet::glmnet(x = x, y = y, lambda = lambda, offset = offset, ...)
# now all predict should work!
}
}
pred.values <- predict(
glmnet.res,
type = "response",
newx = x,
offset = offset,
exact = exact,
s = s
)
if (is.matrix(y)) {
mat.y.ok <- y
} else {
mat.y.ok <- matrix(y, ncol = 1)
}
bern.dev.1 <- t(log(pred.values)) %*% mat.y.ok
bern.dev.2 <- t(log(1 - pred.values)) %*% (1 - mat.y.ok)
model.llik <- (bern.dev.1 + bern.dev.2)
pred.coeffs <- predict(glmnet.res, type = "coefficients", s = s, exact = exact)
models.defs <- data.frame(s = s, loglik = model.llik[, 1])
models.defs$model.id <- 1:nrow(models.defs)
models.defs$df <- apply(pred.coeffs, 2, function(x) sum(x != 0))
models.coeffs <- broom::tidy(pred.coeffs) %>%
dplyr::rename(
term = row,
model.id = column,
estimate = value
)
models.defs2 <- models.defs %>%
dplyr::mutate(AIC = 2 * df - loglik)
list(
models = models.defs2,
coeffs = models.coeffs
)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.