sparseiCov: Sparse/penalized estimators of covariance matrices

In zdk123/SpiecEasi: Sparse Inverse Covariance for Ecological Statistical Inference

Sparse/penalized estimators of covariance matrices

Description

This function estimates the sparse inverse covariance matrix/matrices given data (typically after clr transformation) and further arguments to huge package functions.

Usage

sparseiCov(data, method, npn = FALSE, verbose = FALSE, cov.output = TRUE, ...)

Arguments

data	data matrix with features/OTUs in the columns and samples in the rows. Should be transformed by clr for meaningful results, if the data is compositional
method	estimation method to use as a character string. Currently either 'glasso' or 'mb' (meinshausen-buhlmann)
npn	perform Nonparanormal (npn) transformation before estimation?
verbose	print progress to standard out
cov.output	return the covariance matrix as well.
...	further arguments to huge/estimation functions. See details.

Details

This is a wrapper function for sparse iCov estimations performed by glasso in the huge package.

Therefore, arguments ... should be named. Typically, these are for specifying a penalty parameter, lambda, or the number of penalties to use. By default 10 pentalties are used, ranging logarithmically between lambda.min.ratio*MAX and MAX. Max is the theoretical upper bound on lambda and us max|S|, the maximum absolute value in the data correlation matrix. lambda.min.ratio is 1e-3 by default. Lower values of lambda require more memory/cpu time to compute, and sometimes huge will throw an error.

The argument nlambda determines the number of penalties - somewhere between 10-100 is usually good, depending on how the values of empirical correlation are distributed.

Examples

# simulate data with 1 negative correlation
 set.seed(10010)
 Sigma <- diag(50)*2
 Sigma[1,2] <- Sigma[2,1] <- -.9
 data  <- exp(rmvnorm(50, runif(50, 0, 2), Sigma))

# normalize
 data.f   <- t(apply(data, 1, norm_to_total))
 data.clr <- t(clr(data.f, 1))

# estimate
 est.clr  <- sparseiCov(data.clr, method='glasso')
 est.f    <- sparseiCov(data.f, method='glasso')
 est.log  <- sparseiCov(log(data), method='glasso')

# visualize results
 par(mfrow=c(1,3))
 image(as.matrix(est.log$path[[3]][1:5,1:5]))
 image(as.matrix(est.clr$path[[3]][1:5,1:5]))
 image(as.matrix(est.f$path[[3]][1:5,1:5]))

zdk123/SpiecEasi documentation built on Oct. 20, 2023, 6:49 a.m.