Introduction to rags2ridges

Package installation

The README details how to install the rags2ridges package. When installed, the package is loaded as seen below where we also define a function for adding pretty names to a matrix.

library(rags2ridges)

Small theoretical primer and package usage

The sample variance-covariance matrix, or simply covariance matrix, is well-known and ubiquitous. It is given by

$$ S = \frac{1}{n - 1}XX^T $$

where X is the n × p data matrix that is zero-centered with each p-dimensional observations in the rows. I.e. each row of X is an observation and each column is feature. Often high-dimensional data is organised this way (or transposed).

That X is zero-centered simply means that the column means has been subtracted the columns. The very similar estimate $S = \frac{1}{n}XX^T$ without Bessel’s correction is the maximum likelihood estimate in a multivariate normal model with mean 0 and covariance Σ. The likelihood function in this case is given by

ℓ(Ω; S) = ln |Ω| − tr(SΩ)

where Ω = Σ⁻¹ is the so-called precision matrix (also sometimes called the concentration matrix). It is precisely this Ω for which we seek an estimate we will denote P. Indeed, one can naturally try to use the inverse of S for this:

P = S⁻¹

Let’s try.

The createS() function can easily simulate covariance matrices. But we go a more verbose route for illustration:

p <- 6
n <- 20
X <- createS(n = n, p = p, dataset = TRUE)
head(X, n = 4) # Show 4 first of the n rows

##           A      B      C      D      E     F
## [1,]  0.812 -0.451  1.039 -0.832 -0.192 0.601
## [2,]  0.578  0.837 -0.541  1.199  0.935 0.727
## [3,]  0.139  0.361  1.433 -1.159  0.639 1.091
## [4,] -0.980  1.685 -0.581  0.448  0.971 1.810

Here the columns corresponds to features A, B, C, and so on.

When can then arrive a the MLE using covML() which centers X (subtracting the column means) and then computes the estimate:

S <- covML(X)
print(S)

##         A      B       C      D      E       F
## A  0.9974 -0.391  0.1427 -0.167 0.1557 -0.0166
## B -0.3907  1.153  0.1556  0.343 0.5601  0.2420
## C  0.1427  0.156  0.9868 -0.154 0.0426  0.3714
## D -0.1668  0.343 -0.1544  1.027 0.3172  0.3144
## E  0.1557  0.560  0.0426  0.317 1.0402  0.1828
## F -0.0166  0.242  0.3714  0.314 0.1828  0.8179

Using cov2cor() the well-known correlation matrix could be obtained.

By default, createS() simulates zero-mean i.i.d. normal variables (corresponding to Σ = Ω = I being the identity matrix), but it has plenty of possibilities for more intricate covariance structures. The S matrix could have been obtained directly had we omitted the dataset argument, leaving it to be the default FALSE. The rmvnormal() function is utilized by createS() to generate the normal sample.

We can obtain the precision estimate P using solve() to invert S:

P <- solve(S)
print(P)

##         A       B      C      D       E       F
## A  1.5225  0.8673 -0.319  0.124 -0.7255  0.0338
## B  0.8673  1.7590 -0.373 -0.181 -0.9993 -0.0404
## C -0.3191 -0.3729  1.487  0.485  0.1801 -0.7981
## D  0.1235 -0.1815  0.485  1.405 -0.2552 -0.6469
## E -0.7255 -0.9993  0.180 -0.255  1.6928 -0.0811
## F  0.0338 -0.0404 -0.798 -0.647 -0.0811  1.8646

That’s it! Everything goes well here only because n < p. However, when p is close to n, the estimate become unstable and varies wildly and when p exceeds n one can no longer invert S and this strategy fails:

p <- 25
S2 <- createS(n = n, p = p)  # Direct to S
P2 <- solve(S2)

## Error in solve.default(S2): system is computationally singular: reciprocal condition number = 2.57348e-19

Note that this is now a 25 × 25 precision matrix we are trying to estimate. Datasets where p > n are starting to be common, so what now?

To solve the problem, rags2ridges adds a so-called ridge penalty to the likelihood above — this method is also called L₂ shrinkage and works by “shrinking” the eigenvalues of S in a particular manner to combat that they “explode” when p ≥ n.

The core problem that rags2ridges solves is that

$$ \ell(\Omega; S) = \ln|\Omega| - \text{tr}(S\Omega) - \frac{\lambda}{2}|| \Omega - T||^2_2 $$ where λ > 0 is the ridge penalty parameter, T is a p × p known target matrix (which we will get back to) and || ⋅ ||₂ is the L₂-norm. The maximizing solution here is surprisingly on closed form, but it is rather complicated¹. Assume for now the target matrix is an all zero matrix and thus out of the equation.

The core function of rags2ridges is ridgeP which computes this estimate in a fast manner.

P2 <- ridgeP(S2, lambda = 1.17)
print(P2[1:7, 1:7]) # Showing only the 7 first cols and rows

##         A       B       C       D       E       F      G
## A  3.8762 -0.2152 -0.3045 -0.0924 -0.0926  0.2093 -0.149
## B -0.2152  3.9289 -0.3571  0.0133 -0.1166  0.2722  0.326
## C -0.3045 -0.3571  3.4616 -0.2600 -0.2379  0.0798  0.234
## D -0.0924  0.0133 -0.2600  3.5087 -0.1583 -0.2840  0.325
## E -0.0926 -0.1166 -0.2379 -0.1583  3.5242 -0.1522 -0.355
## F  0.2093  0.2722  0.0798 -0.2840 -0.1522  3.4746 -0.124
## G -0.1492  0.3265  0.2338  0.3247 -0.3547 -0.1239  3.238

And voilà, we have our estimate. We will now discuss the penalty parameters and target matrix and how to choose them.

The penalty parameter

The penalty parameter λ (lambda) shrinks the values of P such toward 0 (when T = 0) — i.e. very larges values of λ makes P “small” and more stable whereas smaller values of λ makes the P tend toward the (possibly non-existent) S⁻¹. So what lambda should you choose? One strategy for choosing λ is selecting it to be stable yet precise (a bias-variance trade-off). Automatic k-fold cross-validation can be done with optPenalty.kCVauto()is well suited for this:

Y <- createS(n, p, dataset = TRUE)
opt <- optPenalty.kCVauto(Y, lambdaMin = 0.001, lambdaMax = 100)
str(opt)

## List of 2
##  $ optLambda: num 0.159
##  $ optPrec  : 'ridgeP' num [1:25, 1:25] 3.3759 -0.3951 0.0732 1.0555 0.5301 ...
##   ..- attr(*, "lambda")= num 0.159
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:25] "A" "B" "C" "D" ...
##   .. ..$ : chr [1:25] "A" "B" "C" "D" ...

As seen, the function returns a list with the optimal penalty parameter and corresponding ridge precision estimate. By default, the the functions performs leave-one-out cross validation. See ?optPenalty.kCVauto` for more information.

The target matrix

The target matrix T is a matrix the same size as P which the estimate is “shrunken” toward — i.e. for large values of λ the estimate goes toward T. The choice of the target is another subject. While one might first think that the all-zeros T = [0] would be a default it is intuitively not a good target. This is because we’d like an estimate that is positive definite (the matrix-equivalent to at positive number) and the null-matrix is not positive definite.

If one has a very good prior estimate or some other information this might used to construct the target. E.g. the function kegg.target() utilizes the Kyoto Encyclopedia of Genes and Genomes (KEGG) database of gene and gene-networks together with pilot data to construct a target.

In the absence of such knowledge, the default could be a data-driven diagonal matrix. The function default.target() offers some different approaches to selecting this. A good choice here is often the diagonal matrix times the reciprocal mean of the eigenvalues of the sample covariance as entries. See ?default.target for more choices.

Gaussian graphical modeling and post processing

What is so interesting with the precision matrix anyway? I’m always interested in correlations and thus the correlation matrix.

As you may know, correlation does not imply causation. Nor does covariance imply causation. However, precision matrix provides stronger hints at causation. A relatively simple transformation of P maps it to partial correlations—much like how the sample covariance S easily maps to the correlation matrix. More precisely, the ijth partial correlation is given by

$$ \rho_{ij|\text{all others}} = \frac{- p_{ij}}{\sqrt{p_{ii}p_{jj}}} $$

where p_ij is the ijth entry of P.

Partial correlations measure the linear association between two random variables whilst removing the effect of other random variables; in this case, it is all the remaining variables. This somewhat absolves the issue in “regular” correlations where are often correlated but only indirectly; either by sort of ‘transitivity of correlations’ (which does not hold generally and is not² so³ simple⁴) or by common underlying variables.

OK, but what is graphical about the graphical ridge estimate?

In a multivariate normal model, p_ij = p_ji = 0 if and only if X_i and X_j are conditionally independent when condition on all other variables. I.e. X_i and X_j are conditionally independent given all X_k where k ≠ i and k ≠ j if and when the ijth and jith elements of P are zero. In real world applications, this means that P is often relatively sparse (lots of zeros). This also points to the close relationship between P and the partial correlations.

The non-zero entries of the a symmetric PD matrix can them be interpreted the edges of a graph where nodes correspond to the variables.

Graphical ridge estimation? Why not graphical Lasso?

The graphical lasso (gLasso) is the L1-equivalent to graphical ridge. A nice feature of the L1 penalty automatically induces sparsity and thus also select the edges in the underlying graph. The L2 penalty of rags2ridges relies on an extra step that selects the edges after P is estimated. While some may argue this as a drawback (typically due to a lack of perceived simplicity), it is often beneficial to separate the “variable selection” and estimation.

First, a separate post-hoc selection step allows for greater flexibility.

Secondly, when co-linearity is present the L1 penalty is “unstable” in the selection between the items. I.e. if 2 covariances are co-linear only one of them will typically be selected in a unpredictable way whereas the L2 will put equal weight on both and “average” their effect. Ultimately, this means that the L2 estimate is typically more stable than the L1.

At last point to mention here is also that the true underlying graph might not always be very sparse (or sparse at all).

How do I select the edges then?

The sparsify() functions lets you select the non-zero entries of P corresponding to edges. It supports a handful different approaches ranging from simple thresholding to false discovery rate based selection.

After edge select GGMnetworkStats() can be utilized to get summary statistics of the resulting graph topology.

References

1.. van Wieringen, W.N. and Peeters, C.F.W. (2016). Ridge Estimation of Inverse Covariance Matrices from High-Dimensional Data. Computational Statistics & Data Analysis, vol. 103: 284-303.

2. Bilgrau, A.E., Peeters, C.F.W., Eriksen, P.S., Boegsted, M., and van Wieringen, W.N. (2015). Targeted Fused Ridge Estimation of Inverse Covariance Matrices from Multiple High-Dimensional Data Classes. arXiv:1509.07982 [stat.ME].