In josue-rodriguez/glaxo: Relaxed Glasso

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%",
  message =  FALSE,
  warning = FALSE
)

glaxo

The glaxo package provides an implementation of the relaxed lasso for Gaussian graphical models based on the simple algorithm described in section 2.2 of @meinshausen2007relaxed. Because the classical lasso introduces bias into the parameters, the motivation behind the relaxed lasso is to debias the estimates.

Installation

You can install the released version of glaxo from CRAN with:

install.packages("glaxo")

And the development version from GitHub with:

# install.packages("devtools")
devtools::install_github("josue-rodriguez/glaxo")

Example

The main function in this package is called glaxo (relaxed + glasso). In it's most basic form, it takes a data.frame or matrix and estimates the relaxed lasso solution for the graphical model

library(glaxo)

Y <- psych::bfi[, 1:10]
Y <- na.omit(Y)

relaxed_glasso <- glaxo(Y, progress = FALSE)

relaxed_glasso

The resulting object can also be plotted using standard ggplot2 syntax

library(ggnetwork)

plot_edges(relaxed_glasso, layout = "circle") +
  scale_size(range = c(0,2)) +
  guides(size = FALSE) +
  scale_color_manual(values = c("#FAA43A", "#5DA5DA"),
                     name = "",
                     labels = c("Negative", "Positive")) +
  theme_facet() +
  theme(legend.position = "top")

Predictions can be made on new data using the relationship between Gaussian graphical models and linear regression [e.g., @stevens1998inverse]

preds <- predict(relaxed_glasso, newdata = Y[1001:2000, ])
str(preds)
round(head(preds$predictions), 3)
round(preds$beta_matrix, 3)

Relaxed lasso vs regular lasso

The intuition behind the relaxed lasso is that the lasso algorithm is run on the data once for an initial screening and then again to lower the false positive rate (i.e., 1 - specificity). This can result in better performance than the traditional lasso solution

library(GGMncv)

main <- GGMnonreg::gen_net()
n <- 5000
Y <- MASS::mvrnorm(n, mu = rep(0, 20), main$cors)

regular_glasso <- GGMncv::ggmncv(cor(Y),
                                 n = nrow(Y), 
                                 penalty = "lasso", 
                                 ic = "bic",
                                 progress = FALSE)

relaxed_glasso <- glaxo(Y, progress = FALSE, ic = "bic")


glaxo:::performance(Estimate = regular_glasso$adj, True = main$adj)
glaxo:::performance(Estimate = relaxed_glasso$adj, True = main$adj)