Inference of Multiscale Gaussian Graphical Model.

Cluster variables using L2 fusion penalty and simultaneously estimates a gaussian graphical model structure with the addition of L1 sparsity penalty.

Usage

mglasso(
  x,
  lambda1 = 0,
  fuse_thresh = 0.001,
  maxit = NULL,
  distance = c("euclidean", "relative"),
  lambda2_start = 1e-04,
  lambda2_factor = 1.5,
  precision = 0.01,
  weights_ = NULL,
  type = c("initial"),
  compact = TRUE,
  verbose = FALSE
)

Arguments

x: Numeric matrix ($n x p$). Multivariate normal sample with $n$ independent observations.
lambda1: Positive numeric scalar. Lasso penalty.
fuse_thresh: Positive numeric scalar. Threshold for clusters fusion.
maxit: Integer scalar. Maximum number of iterations.
distance: Character. Distance between regression vectors with permutation on symmetric coefficients.
lambda2_start: Numeric scalar. Starting value for fused-group Lasso penalty (clustering penalty).
lambda2_factor: Numeric scalar. Step used to update fused-group Lasso penalty in a multiplicative way..
precision: Tolerance for the stopping criterion (duality gap).
weights_: Matrix of weights.
type: If "initial" use classical version of MGLasso without weights.
compact: Logical scalar. If TRUE, only save results when previous clusters are different from current.
verbose: Logical scalar. Print trace. Default value is FALSE.

Value

A list-like object of class mglasso is returned.

out: List of lists. Each element of the list corresponds to a clustering level. An element named levelk contains the regression matrix beta and clusters vector clusters for a clustering in k clusters. When compact = TRUE out has as many elements as the number of unique partitions. When set to FALSE, the function returns as many items as the the range of values taken by lambda2.
l1: the sparsity penalty lambda1 used in the problem solving.

Details

Estimates a gaussian graphical model structure while hierarchically grouping variables by optimizing a pseudo-likelihood function combining Lasso and fuse-group Lasso penalties. The problem is solved via the COntinuation with NEsterov smoothing in a Shrinkage-Thresholding Algorithm (Hadj-Selem et al. 2018). Varying the fusion penalty $\lambda_2$ in a multiplicative fashion allow to uncover a seemingly hierarchical clustering structure. For $\lambda_2 = 0$, the approach is theoretically equivalent to the Meinshausen-Bühlmann (2006) neighborhood selection as resuming to the optimization of pseudo-likelihood function with $\ell_1$ penalty (Rocha et al., 2008). The algorithm stops when all the variables have merged into one cluster. The criterion used to merge clusters is the $\ell_2$-norm distance between regression vectors.

For each iteration of the algorithm, the following function is optimized : $$1/2 \sum_{i=1}^p || X^i - X^{\ i} \beta^i ||_2 ^2 + \lambda_1 \sum_{i = 1}^p || \beta^i ||_1 + \lambda_2 \sum_{i < j} || \beta^i - \tau_{ij}(\beta^j) ||_2.$$

where $\beta^i$ is the vector of coefficients obtained after regression $X^i$ on the others and $\tau_{ij}$ is a permutation exchanging $\beta_j^i$ with $\beta_i^j$.

Examples

# \donttest{
install_conesta()
#> mglasso requires the r-reticulate virtual environment. Attempting to create...
#> virtualenv: r-reticulate
#> Using virtual environment 'r-reticulate' ...
#> + '/home/runner/.virtualenvs/r-reticulate/bin/python' -m pip install --upgrade 'scipy == 1.7.1' 'scikit-learn'
#> Installing pylearn-parsimony
#> pylearn-parsimony is installed.
n = 50
K = 3
p = 9
rho = 0.85
blocs <- list()
for (j in 1:K) {
  bloc <- matrix(rho, nrow = p/K, ncol = p/K)
  for(i in 1:(p/K)) { bloc[i,i] <- 1 }
  blocs[[j]] <- bloc
}

mat.covariance <- Matrix::bdiag(blocs)
mat.covariance
#> 9 x 9 sparse Matrix of class "dgCMatrix"
#>                                                   
#>  [1,] 1.00 0.85 0.85 .    .    .    .    .    .   
#>  [2,] 0.85 1.00 0.85 .    .    .    .    .    .   
#>  [3,] 0.85 0.85 1.00 .    .    .    .    .    .   
#>  [4,] .    .    .    1.00 0.85 0.85 .    .    .   
#>  [5,] .    .    .    0.85 1.00 0.85 .    .    .   
#>  [6,] .    .    .    0.85 0.85 1.00 .    .    .   
#>  [7,] .    .    .    .    .    .    1.00 0.85 0.85
#>  [8,] .    .    .    .    .    .    0.85 1.00 0.85
#>  [9,] .    .    .    .    .    .    0.85 0.85 1.00

set.seed(11)
X <- mvtnorm::rmvnorm(n, mean = rep(0,p), sigma = as.matrix(mat.covariance))
X <- scale(X)

res <- mglasso(X, 0.1, lambda2_start = 0.1)
res$out[[1]]$clusters
#> [1] 1 2 3 4 5 6 7 8 9
res$out[[1]]$beta
#>              [,1]        [,2]         [,3]        [,4]         [,5]        [,6]
#>  [1,]  0.00000000  0.50942304  0.420856249  0.08059251  0.198555576 -0.15744283
#>  [2,]  0.57897509  0.00000000  0.375480194 -0.11129404 -0.024338821  0.03622193
#>  [3,]  0.50877541  0.39402762  0.000000000  0.09138873 -0.244767942  0.13197185
#>  [4,]  0.09644133 -0.11232346  0.096386587  0.00000000  0.380559101  0.53949430
#>  [5,]  0.34031463 -0.03189842 -0.367833554  0.55046923  0.000000000  0.28792248
#>  [6,] -0.23822144  0.05604748  0.156337318  0.67569926  0.247844203  0.00000000
#>  [7,]  0.07147603  0.00000000 -0.066314827  0.19093102  0.005503117 -0.21086078
#>  [8,] -0.17121934  0.13577447  0.006274291  0.04037628  0.017037923  0.01709042
#>  [9,]  0.00000000 -0.03662807  0.030397589 -0.27537816  0.024525025  0.21891289
#>               [,7]         [,8]        [,9]
#>  [1,]  0.088104950 -0.136441557  0.00000000
#>  [2,] -0.005820201  0.112868548 -0.04462052
#>  [3,] -0.078583570  0.004854681  0.03876315
#>  [4,]  0.274032064  0.031338363 -0.32834722
#>  [5,]  0.011236774  0.027884703  0.06051871
#>  [6,] -0.374671576  0.008336123  0.32264027
#>  [7,]  0.000000000  0.389473752  0.56690235
#>  [8,]  0.592191132  0.000000000  0.30687389
#>  [9,]  0.688368824  0.245147593  0.00000000
# }