JGL_AIC_widesearch: Fits a Joing Graphical Lasso model. The lambda 1 and lambda 2...

In GiulioCostantini/JGL2: JGL2: Joint Graphical Lasso, select the lambda parameters using AIC or crossvalidation

Description

Find the minimum value of a criteiron defined by aicfun (e.g., AIC) using 2-dimensional optimization. Since the problem is not convex, considering many starting values for l1 and l2 helps escaping local minima, therefore ncand values of lambda1 x ncand values of lambda2 are considered. This function should perform better than JGL_AIC_sequentialsearch, however it is also much more computationally intensive. The function is optimized for parallel computing on Windows machines. on Other OSs parallelization may not work, in this case, use ncores = 1.

Usage

JGL_AIC_widesearch(dat, splt, return.whole.theta = TRUE, l1min = 0, l1max = 1, l2min = 0, l2max = 1, ncand = 20, criterion = c("ebic","aic"), gamma = 0.5, optmethod = "CG", ncores = 1, ...)

Arguments

dat	A dataset that includes the variables on which the gaussian graphical models should be computed, plus an additional factor "splt" which defines different classes.
splt	Character string. The name of the variable in dat that defines different classes
return.whole.theta	Logical. Parameter passed directly to JGL, determines whether the whole concentration matrix should be returned in output. See JGL
l1min	Numeric. Minimum value of the lasso parameter lambda1
l1max	Numeric. Maximum value of the lasso parameter lambda1
l2min	Numeric. Minimum value of the lasso parameter lambda2
l2max	Numeric. Maximum value of the lasso parameter lambda2
ncand	Integer. number of values for lambda 1 and for lambda 2
criterion	The criterion for selecting the lambda values. Can be "aic" for the Akaike information criterion (AIC) or "ebic" for the Extended Bayes Information Criterion (EBIC).
gamma	The gamma value for the EBIC criterion. A value of 0 results in the BIC.
optmethod	character string. Determines the optimization method. See parameter method in function optim
ncores	Number of cores to use. The function is optimized for parallel computing in Windows, parallel computing may not work on other systems.
...	Other parameters for link{JGL}. This feature is still in development, it may not work as intended, especially if using parallel computing.

Value

jgl	The output of JGL
lambda1	the value of lambda 1 that minimizes the output of aicfun
l1theormax	the minimal value of lambda 1 that would result in at least one completely disconnected network (all missing edges). Values of lambda 1 > l1theormax are not considered
l2theormax	for each candidate value of lambda 1, the minimal value of lambda 2 that would make the networks in the different classes all equal to each other. For each lambda 1, values of lambda 2 > l2theormax are not considered
aic	The Akaike Information Criterion (or another criterion to minimize passed to aicfun) corresponding to the solution in jgl
critsFor3dPlot	Please ignore this output, only for internal verification purposes. A list of the lambda 1 and 2 values identified from each starting pair, as well as the corresponding AIC

Author(s)

Giulio Costantini

References

Danaher, P., Wang, P., and Witten, D. M. (2014). The joint graphical lasso for inverse covariance estimation across multiple classes. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 76(2), 373-397. http://doi.org/10.1111/rssb.12033

See Also

JGL, JGL_AIC_sequentialsearch, JGL_cv, JGL_AIC_surfaceplot

Examples

## Not run: 
N <- 1000 # sample size
sigma1 <- matrix(c(1, .5, 0, 0,
                   .5, 1, .2, 0,
                   0, .2, 1, 0,
                   0, 0, 0, 1), ncol = 4)

sigma2 <- matrix(c(1, .5, .4, .4,
                   .5, 1, .2, 0,
                   .4, .2, 1, 0,
                   .4, 0, 0, 1), ncol = 4)

dat <- list()
dat[[1]] <- MASS::mvrnorm(n = N, mu = rep(0, ncol(sigma1)), Sigma = sigma1)
dat[[2]] <- MASS::mvrnorm(n = N, mu = rep(0, ncol(sigma2)), Sigma = sigma2)
lapply(dat, function(x) corpcor::cor2pcor(cor(x)))
dat <- data.frame(rbind(dat[[1]], dat[[2]]))
dat$splt <- c(rep(1, N), rep(2, N))
system.time(
cg1 <- JGL_AIC_widesearch(dat = dat, splt = "splt", ncand = 10, ncores = 7, optmethod = "CG")
)

## End(Not run)

GiulioCostantini/JGL2 documentation built on May 6, 2019, 6:29 p.m.