In MGallow/CVglasso: Lasso Penalized Precision Matrix Estimation

knitr::opts_chunk$set(echo = TRUE)
library(magrittr)

Overview

CVglasso is an R package that estimates a lasso-penalized precision matrix via block-wise coordinate descent -- also known as the graphical lasso (glasso) algorithm. This package is a simple wrapper around the popular glasso package and extends and enhances its capabilities. These enhancements include built-in cross validation and visualizations.

A (possibly incomplete) list of functions contained in the package can be found below:

CVglasso() computes the estimated precision matrix
plot.CVglasso() produces a heat map or line graph for cross validation errors

See package website or manual.

Installation

# The easiest way to install is from CRAN
install.packages("CVglasso")

# You can also install the development version from GitHub:
# install.packages("devtools")
devtools::install_github("MGallow/CVglasso")

If there are any issues/bugs, please let me know: github. You can also contact me via my website. Pull requests are welcome!

Usage

library(CVglasso)
set.seed(123)

# generate data from a sparse oracle precision matrix
# we will try to estimate this matrix from the data

# first compute the oracle covariance matrix
S = matrix(0.7, nrow = 5, ncol = 5)
for (i in 1:5){
  for (j in 1:5){
    S[i, j] = S[i, j]^abs(i - j)
  }
}

# print oracle precision matrix (shrinkage might be useful)
(Omega = round(qr.solve(S), 3))

# generate 100 x 5 matrix with rows drawn from iid N_p(0, S)
set.seed(123)
Z = matrix(rnorm(100*5), nrow = 100, ncol = 5)
out = eigen(S, symmetric = TRUE)
S.sqrt = out$vectors %*% diag(out$values^0.5) %*% t(out$vectors)
X = Z %*% S.sqrt

# calculate sample covariance
sample = (nrow(X) - 1)/nrow(X)*cov(X)

# print sample precision matrix (perhaps a bad estimate)
round(qr.solve(cov(X)), 5)


# CVglasso (lam = 0.5)
CVglasso(S = sample, lam = 0.5)

# cross validation
(CV = CVglasso(X, trace = "none"))

# produce line graph for CV errors for CVGLASSO
plot(CV)

# produce CV heat map for CVGLASSO
plot(CV, type = "heatmap")