pulsar (version 0.3.11)

pulsar-function: Graphical model functions for pulsar

Description

Correctly specify a function for graphical model estimation that is compatible with the pulsar package.

Details

It is easy to construct your own function for penalized model estimation that can be used with this package. The R function must have correctly specified inputs and outputs and is passed into the fun argument to pulsar or batch.pulsar. Any function that does not follow these rules will fail to give the desired output and may trigger an error.

These packages on CRAN have functions that work out of the box, so you won't need to construct a wrapper:

~function~	~package~
huge	huge
sugm	flare

Inputs:

The function may take arbitrary, named arguments but the first argument must be the data \(n*p\) data matrix with the \(n\) samples in rows and \(p\) features in the columns. At least one argument must be named "lambda", which is expected to be a decreasing numeric vector of penalties. The non-data arguments should be passed into pulsar or batch.pulsar as a named list (the names must match function arguments exactly) to the fargs argument.

Outputs:

The output from the function must be a list or another S3 object inherited from a list. At least one member must be named path. This path object itself must be a list of \(p*p\) adjacency matrices, one for each value of lambda. Each cell in the adjacency matrix contains a 1 or TRUE if there is an edge between two nodes or 0/FALSE otherwise. It is highly recommended (though not enforced by pulsar) that each adjacency matrix be a column-oriented, compressed, sparse matrix from the Matrix package. For example, dgCMatrix/dsCMatrix (general/symmetric numeric Matrix) or the 1-bit lgCMatrix/lsCMatrix classes. The function may return other named outputs, but these will be ignored.

References

Müller, C. L., Bonneau, R. A., & Kurtz, Z. D. (2016). Generalized Stability Approach for Regularized Graphical Models. arXiv: https://arxiv.org/abs/1605.07072.

See Also

pulsar, batch.pulsar, huge, Matrix

Examples

## Generate a hub example
 dat <- huge::huge.generator(100, 40, 'hub', verbose=FALSE)

## Simple correlation thresholding
corrthresh <- function(data, lambda) {
  S <- cor(data)
  path <- lapply(lambda, function(lam) {
    tmp <- abs(S) > lam
    diag(tmp) <- FALSE
    as(tmp, 'lMatrix')
  })
  list(path=path)
}

## Inspect output
lam <- getLamPath(getMaxCov(dat$sigmahat), 1e-4, 10)
out.cor  <- pulsar(dat$data, corrthresh, fargs=list(lambda=lam))
out.cor

if (FALSE) {
## Additional examples
## quic
library(QUIC)
quicr <- function(data, lambda, ...) {
    S    <- cov(data)
    est  <- QUIC(S, rho=1, path=lambda, msg=0, tol=1e-2, ...)
    est$path <-  lapply(seq(length(lambda)), function(i) {
                   ## convert precision array to adj list
                   tmp <- est$X[,,i]; diag(tmp) <- 0
                 as(tmp!=0, "lMatrix")
    })
    est
}
## clime
library(clime)
climer <- function(data, lambda, tol=1e-5, ...) {
     est <- clime(data, lambda, ...)
     est$path <- lapply(est$Omegalist, function(x) {
                     diag(x) <- 0
                     as(abs(x) > tol, "lMatrix")
                 })
     est
}

## inverse cov shrinkage Schafer and Strimmer, 2005
library(corpcor)
icovshrink <- function(data, lambda, tol=1e-3, ...) {
     path <- lapply(lambda, function(lam) {
                     tmp <- invcov.shrink(data, lam, verbose=FALSE)
                     diag(tmp) <- 0
                     as(abs(tmp) > tol, "lMatrix")
                 })
     list(path=path)
}

## Penalized linear model, only
library(glmnet)
lasso <- function(data, lambda, respind=1, family="gaussian", ...) {
         n <- length(lambda)
         tmp <- glmnet(data[,-respind], data[,respind],
                                   family=family, lambda=lambda, ...)
         path <-lapply(1:n, function(i) as(tmp$beta[,i,drop=FALSE], "lMatrix"))
         list(path=path)
}

## alternative stability selection (DIFFERENT from hdi package)
out <- pulsar(dat$data, lasso, fargs=list(lambda=lam))
mergmat <- do.call('cbind', tmp$stars$merge)
image(mergmat)
}