pulsar-function | R Documentation |
Correctly specify a function for graphical model estimation that is compatible with the pulsar package.
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.
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.
pulsar
, batch.pulsar
, huge, Matrix
## 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
## Not run:
## 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)
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.