Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
58,268 characters · 17 sections · 35 citation commands
Tree-based Node Aggregation in Sparse Graphical Models
\paragraph{Abstract.} High-dimensional graphical models are often estimated using regularization that is aimed at reducing the number of edges in a network. In this work, we show how even simpler networks can be produced by aggregating the nodes of the graphical model. We develop a new convex regularized method, called the {\em tree-aggregated graphical lasso} or tag-lasso, that estimates graphical models that are both edge-sparse and node-aggregated. The aggregation is performed in a data-driven fashion by leveraging side information in the form of a tree that encodes node similarity and facilitates the interpretation of the resulting aggregated nodes. We provide an efficient implementation of the tag-lasso by using the locally adaptive alternating direction method of multipliers and illustrate our proposal's practical advantages in simulation and in applications in finance and biology.
\paragraph{Keywords.} aggregation, graphical model, high-dimensionality, regularization, sparsity
Graphical models are greatly useful for understanding the relationships among large numbers of variables. Yet, estimating graphical models with many more parameters than observations is challenging, which has led to an active area of research on high-dimensional inverse covariance estimation. Numerous methods attempt to curb the curse of dimensionality through regularized estimation procedures (e.g., meinshausen2006high, yuan2007model, banerjee2008model, friedman2008sparse, rothman2008sparse, peng2009partial, yuan2010high, cai2011constrained, cai2016estimating). Such methods aim for sparsity in the inverse covariance matrix, which corresponds to graphical models with only a small number of edges. A common method for estimating sparse graphical models is the graphical lasso (glasso) yuan2007model, banerjee2008model, rothman2008sparse, friedman2008sparse, which adds an $\ell_1$-penalty to the negative log-likelihood of a sample of multivariate normal random variables. While this and many other methods focus on the edges for dimension reduction, far fewer contributions (e.g., tan2015cluster, eisenach2020high, pircalabelu2020community) focus on the {\em nodes} as a guiding principle for dimension reduction.
Nonetheless, node dimension reduction is becoming increasingly relevant in many areas where data are being measured at finer levels of granularity. For instance, in biology, modern high-throughput sequencing technologies provide low-cost microbiome data at high resolution; in neuroscience, brain activity in hundreds of regions of interest can be measured; in finance, data at the individual company level at short time scales are routinely analyzed; and in marketing, joint purchasing data on every stock-keeping-unit (product) is recorded. The fine-grained nature of this data brings new challenges. The sheer number of fine-grained, often noisy, variables makes it difficult to detect dependencies. Moreover, there can be a mismatch between the resolution of the measurement and the resolution at which natural meaningful interpretations can be made. The purpose of an analysis may be to draw conclusions about entities at a coarser level of resolution than happened to be measured. Because of this mismatch, practitioners are sometimes forced to devise ad hoc post-processing steps involving, for example, coloring the nodes based on some classification of them into groups in an attempt to make the structure of an estimated graphical model more interpretable and the domain-specific takeaways more apparent (e.g., millington2019quantifying).
Our solution to this problem is to incorporate the side information about the relationship between nodes directly into the estimation procedure. In our framework, this side information is encoded as a tree whose leaves correspond to the measured variables. Such tree structures are readily available in many domains (e.g., taxonomies in biology and hierarchical classifications of jobs, companies, and products in business) and is well-suited to expressing multi-resolution structure that is present in many problems. We propose a new convex regularization procedure, called tag-lasso, which stands for tree-aggregated-graphical-lasso. This procedure combines node (or variable) aggregation with edge-sparsity. The tree-based aggregation serves to both amplify the signal of similar, low-level variables and render a graphical model involving nodes at an appropriate level of scale to be relevant and interpretable. The edge-sparsity encourages the graphical model involving the aggregated nodes has a sparse network structure.
Our procedure is based on a tree-based parameterization strategy that translates the node aggregation problem into a sparse modeling problem, following an approach previously introduced in the regression setting yan2018rare. In Figure (ref) (to be discussed more thoroughly in Section (ref)), we see that tag-lasso is able to recover the aggregated, sparse graph structure. By doing so, it yields a more accurate estimate of the true graph, and its output is easier to interpret than the full, noisy graph obtained by the glasso.
The rest of the paper is organized as follows. Section (ref) introduces the tree-based parameterization structure for nodewise aggregation in graphical models. Section (ref) introduces the tag-lasso estimator, formulated as a solution to a convex optimization problem, for which we derive an efficient algorithm. Section (ref) presents the results of a simulation study. Section (ref) illustrates the practical advantages of the tag-lasso on financial and microbiome data sets. Section (ref) concludes.
Let $\bf S$ be the empirical covariance matrix based on $n$ multivariate normal observations of dimension $p$, with mean vector $\boldsymbol\mu$ and covariance matrix $\boldsymbol \Sigma$. The target of estimation is the precision matrix $\boldsymbol \Omega = \boldsymbol \Sigma^{-1}$, whose sparsity pattern provides the graph structure of the Gaussian graphical model, since $\Omega_{jk}=0$ is equivalent to variables $j$ and $k$ being conditionally independent given all other variables. To estimate the precision matrix, it is common to use a convex penalization method of the form
where $\text{tr}(\cdot)$ denotes the trace, $\mathcal{P}(\cdot)$ is a convex penalty function, and $\lambda>$ is a tuning parameter controlling the degree of penalization. Choosing the $\ell_1$-norm
where $\boldsymbol \Omega^{-\text{diag}}$ contains the unique off-diagonal elements, yields the {\em graphical lasso} (glasso) friedman2008sparse,yuan2007model, banerjee2008model, rothman2008sparse. It encourages $\widehat{\boldsymbol\Omega}$ to be sparse, corresponding to a graphical model with few edges.
However, when $\boldsymbol\Omega$ is not sparse, demanding sparsity in $\widehat{\boldsymbol\Omega}$ may not be helpful, as we will show in Section (ref). Such settings can arise when data are measured and analyzed at ever higher resolutions (a growing trend in many areas, see e.g. callahan2017exact). A tree is a natural way to represent the different scales of data resolution, and we introduce a new choice for $\mathcal{P}$ that uses this tree to guide node aggregation, thereby allowing for a data adaptive choice of data scale for capturing dependencies. Such tree-based structures are available in many domains. For instance, companies can be aggregated according to hierarchical industry classification codes; products can be aggregated from brands towards product categories; brain voxels can be aggregated according to brain regions; microbiome data can be aggregated according to taxonomy. The resulting penalty function then encourages a more general and yet still highly interpretable structure for $\widehat{\boldsymbol\Omega}$. In the following subsection, we use a toy example to illustrate the power of such an approach.
Consider a toy example with $p$ variables
where $\varepsilon_1,\ldots,\varepsilon_p$ are independent standard normal random variables. By construction, it is clear that there is a very simple relationship between the variables: The first two variables both depend on the sum of the other $p-2$ variables. However, a standard graphical model on the $p$ variables does not naturally express this simplicity. The first row of Table (ref) shows the covariance and precision matrices for the full set of variables $X_1, \ldots, X_p$. The graphical model on the full set of variables is extremely dense $O(p^2)$ edges. Imagine if instead we could form a graphical model with only three variables: $X_1, X_2, \widetilde{X}$, where the last variable $\widetilde{X}=\sum_{j=3}^p X_j$ aggregates all but the first two variables. The bottom row of Table (ref) results in a graphical model that matches the simplicity of the situation.
The lack of sparsity in the $p$-node graphical model means that the graphical lasso will not do well. Nonetheless, a method that could perform node aggregation would be able to yield a highly-interpretable aggregated sparse graphical model since $X_1$ and $X_2$ are conditionally independent given the aggregated variable $\widetilde{X}$.
It is useful to map from the small aggregated graphical model to the original $p$-node graphical model. One does so by writing the precision matrix in “$G$-block" format (bunea2020model, although they introduce this terminology in the context of the covariance matrix, not its inverse) for a given partition $G=\{G_1, ..., G_K\}$ of the nodes $\{1,\ldots, p\}$ and corresponding $p\times K$ membership matrix ${\bf M}$, with entries $M_{jk}=1$ if $j \in G_k$, and $M_{jk}=0$ otherwise. In particular, there exists a $K\times K$ symmetric matrix ${\bf C}$ and a $p\times p$ diagonal matrix ${\bf D}$ such that the precision matrix can be written as $\boldsymbol\Omega= {\bf M} {\bf C} {\bf M}^\top + {\bf D}$. The block-structure of $\boldsymbol \Omega$ is captured by the first part of the decomposition, the aggregated $K\times K$ precision matrix on the set of aggregated nodes can then be written as $ \boldsymbol\Omega_{\text{agg}} = {\bf C} + {\bf D}_{\text{agg}}, $ where ${\bf D}_{\text{agg}} = ({\bf M}^\top {\bf D}^{-1}{\bf M})^{-1}$ is diagonal. In the above example, $K=3$, $G_1 = \{1\}, \ G_2 = \{2\}, \ G_3 = \{3, \ldots, p\}$ and ${\bf M} {\bf C} {\bf M}^\top$ has only three distinct rows/columns since the aggregated variables $j=3, \ldots, p$ share all their entries. In the presence of node aggregation and edge sparsity, the graphical model corresponding to the aggregated precision matrix is far more parsimonious than the graphical model on the full precision matrix (see Table (ref)).
As motivated by this example, our main goal is to estimate the precision matrix in such a way that we can navigate from a $p$-dimensional problem to a $K$-dimensional problem whose corresponding graphical model provides a simple description of the conditional dependency structure among $K$ aggregates of the original variables. In the following proposition, we show that this can be accomplished by looking for a precision matrix that has a $G$-block structure. The proof of the proposition is included in Appendix (ref)
While Proposition (ref) gives us the desired interpretation in the graphical model with $K$ aggregated nodes, \color{black} in practice, the partition $G$, its size $K$, and corresponding membership matrix ${\bf M}$ are, however, unknown. Rather than considering arbitrary partitions of the variables, we constrain ourselves specifically to partitions guided by a known tree. In so doing, we allow ourselves to exploit side information and help ensure that the aggregated nodes will be easily interpretable. To this end, we introduce a tree-based parameterization strategy that allows us to embed the node dimension reduction into a convex optimization framework.
Our aggregation procedure assumes that we have, as side information, a tree that represents the closeness (or similarity) of variables. We introduce here a matrix-valued extension of the tree-based parameterization developed in yan2018rare for the regression setting. We consider a tree $\mathcal{T}$ with $p$ leaves $\boldsymbol \Omega_1, \ldots, \boldsymbol \Omega_p$ where $\boldsymbol \Omega_j$ denotes column $1 \leq j \leq p$ of $\boldsymbol \Omega$. We restrict ourselves to partitions that can be expressed as a collection of branches of $\mathcal{T}$. Newly aggregated nodes are then formed by summing variables within branches. To this end, we assign a $p$-dimensional parameter vector ${\boldsymbol\gamma}_u$ to each node $u$ in the tree $\mathcal{T}$ (see Figure (ref) for an example). Writing the set of nodes in the path from the root to the $j^{\text{th}}$ leaf (variable) as $\text{ancestor}(j) \cup \{j\}$, we express each column/row in the precision matrix as
where we sum over all the $\boldsymbol \gamma_u$'s along this path, and ${\bf e}_j$ denotes the $p$-dimensional vector with all zeros except for its $j^{\text{th}}$ element that is equal to one. In the remainder, we will make extensive use of the more compact notation $ \boldsymbol \Omega = {\bf A}\boldsymbol{\Gamma} + {\bf D},$ where ${\bf A} \in \{0,1\}^{p \times |\mathcal{T}|}$ is a binary matrix with $A_{jk} = 1{\{u_k \in \text{ancestor}(j) \cup \{j\}\}} = 1{\{j \in \text{descendant}(u_k) \cup \{u_k\}\}}$, $\boldsymbol\Gamma$ is a $|\mathcal{T}|\times p$ parameter matrix collecting the $\boldsymbol\gamma_u$'s in its rows and ${\bf D}$ is a diagonal parameter matrix with elements $d_1, \ldots, d_p$.
By zeroing out $\boldsymbol \gamma_u$'s, certain nodes will be aggregated, as can be seen from the illustrative example in Figure (ref). More precisely, let ${\mathcal{Z}} = \{u : {{\boldsymbol \gamma}}_u \neq {\bf 0} \}$ denote the set of non-zero rows in ${\boldsymbol \Gamma}$ and let ${\bf A}_{{\mathcal{Z}}}$ be the sub-matrix of ${\bf A}$ where only the columns corresponding to the non-zeros rows in ${{\boldsymbol \Gamma}}$ are kept. The number of blocks $K$ in the aggregated network is then given by the number of unique rows in ${\bf A}_{{\mathcal{Z}}}$. The membership matrix $\bf M$ (Section (ref)), and hence the set of aggregated nodes, can then be derived from the variables (rows) in the matrix ${\bf A}_{{\mathcal{Z}}}$ that share all their row-entries.
We are now ready to introduce the tag-lasso, which is based on this parameterization.
To achieve dimension reduction via node aggregation and edge sparsity simultaneously, we extend optimization problem (ref) by incorporating the parameterization introduced above. Our estimator, called the tag-lasso, is defined as
with $\|\boldsymbol{\Gamma}_{-r}\|_{2,1} = \sum_{u \in \mathcal{T}_{-r}}\|\boldsymbol \gamma_u\|_2$ and $\mathcal{T}_{-r}$ being the set of all nodes in $\mathcal T$ other than the root. This norm induces row-wise sparsity on all non-root rows of ${\boldsymbol\Gamma}$. This row-wise sparsity, in turn, induces node aggregation as explained in Section (ref). The root is excluded from this penalty term so that in the extreme of large $\lambda_1$ one gets complete aggregation but not necessarily sparsity (in this extreme, all off-diagonal elements of $\widehat{\boldsymbol\Omega}$ are equal to the scalar $\gamma$ that appears in the equality constraint involving $\boldsymbol\gamma_r$). While $\lambda_1$ controls the degree of node aggregation, $\lambda_2$ controls the degree of edge sparsity. When $\lambda_1=0$, the optimization problem in (ref) reduces to the glasso.
Finally, note that optimization problem (ref) fits into the general formulation of penalized graphical models given in (ref) since it can be equivalently expressed as $$ \widehat{\boldsymbol \Omega} = \underset{\boldsymbol \Omega}{\operatorname{argmin}} \{- \text{logdet}(\boldsymbol\Omega) + \text{tr}({\bf S}\boldsymbol{\Omega}) + \lambda_1 \mathcal{P}_\text{aggregate}(\boldsymbol\Omega) + \lambda_2 \mathcal{P}_\text{sparse}(\boldsymbol\Omega) \ \text{s.t.} \ \boldsymbol{\Omega}=\boldsymbol{\Omega}^\top, \boldsymbol \Omega \succ {\bf 0}\}, $$ where $$ \mathcal{P}_\text{aggregate}(\boldsymbol\Omega) = \underset{\boldsymbol\Gamma, {\bf D}}{\operatorname{min}} \ \{\|\boldsymbol{\Gamma}_{-r}\|_{2,1} \ \text{s.t.} \ \boldsymbol\gamma_r = \gamma {\bf 1}_p, \ \boldsymbol \Omega = {\bf A}\boldsymbol{\Gamma} + {\bf D}, \ {\bf D} \ \text{diag}, \ D_{jj} \geq 0 \ \text{for} \ j=1,\ldots,p\} $$ and $\mathcal{P}_\text{sparse}(\boldsymbol\Omega)$ is the $\ell_1$-norm defined in (ref).
We develop an {\em alternating direction method of multipliers} (ADMM) algorithm Boyd11, specifically tailored to solving (ref). Our ADMM algorithm is based on solving this equivalent formulation of (ref):
Additional copies of $\boldsymbol\Omega$ and $\boldsymbol\Gamma$ are introduced to efficiently decouple the optimization problem.
Furthermore, we use an extension called {\em locally adaptive-ADMM} (LA-ADMM, xu2017admm) with adaptive penalization to improve performance. The full details of the algorithm are provided in Appendix (ref).
To select the tuning parameters $\lambda_1$ and $\lambda_2$, we form a $10\times 10$ grid of ($\lambda_1, \lambda_2$) values and find the pair that minimizes a 5-fold cross-validated likelihood-based score,
where $\widehat{\boldsymbol\Omega}_{-\mathcal F_k}$ is an estimate of the precision matrix trained while withholding the samples in the $k^{\text{th}}$ fold and ${\bf S}_{\mathcal F_k}$ is the sample covariance matrix computed on the $k^{\text{th}}$ fold. In particular, we take $\widehat{\boldsymbol\Omega}_{-\mathcal F_k}$ to be a re-fitted version of our estimator (e.g., belloni2013least). After fitting the tag-lasso, we obtain $\widehat{\mathcal{Z}} = \{u : \widehat{{\boldsymbol \gamma}}_u \neq {\bf 0} \},$ the set of non-zero rows in $\widehat{{\boldsymbol \Gamma}}$, which suggests a particular node aggregation; and $ \widehat{\mathcal{P}} = \{ (i,j) : \widehat{{\boldsymbol \Omega}}_{ij} \neq {\bf 0} \},$ the set of non-zero elements in $\widehat{{\boldsymbol \Omega}}$, which suggests a particular edge sparsity structure. We then re-estimate $\boldsymbol{\Omega}$ by maximizing the likelihood subject to these aggregation and sparsity constraints:
We solve this with an LA-ADMM algorithm similar to what is described in Section (ref) and Appendix (ref).
Combined forms of dimension reduction in graphical models can be found in, amongst others, chandrasekaran2010latent, tan2015cluster, eisenach2020high, brownlees2020community, pircalabelu2020community.
chandrasekaran2010latent consider a blend of principal component analysis with graphical modeling by combining sparsity with a low-rank structure. tan2015cluster and eisenach2020high both propose two-step procedures that first cluster variables in an initial dimension reduction step and subsequently estimate a cluster-based graphical model. brownlees2020community introduce partial correlation network models with community structures but rely on the sample covariance matrix of the observations to perform spectral clustering. Our procedure differs from these works by introducing a single convex optimization problem that simultaneously induces aggregation and edge sparsity for the precision matrix.
Our work is most closely related to pircalabelu2020community who estimate a penalized graphical model and simultaneously classify nodes into communities. However, pircalabelu2020community do not use tree-based node-aggregation. Our approach, in contrast, considers the tree $\mathcal{T}$ as an important part of the problem to help determine the extent of node aggregation, and as a consequence the number of aggregated nodes (i.e. clusters, communities or blocks) $K$, in a data-driven way through guidance of the tree-based structure on the nodes.
We investigate the advantages of jointly exploiting node aggregation and edge sparsity in graphical models. To this end, we compare the performance of the tag-lasso to two benchmarks:
All simulations were performed using the simulator package bien2016simulator in R Rcoreteam. We evaluate the estimators in terms of three performance metrics: estimation accuracy, aggregation performance, and sparsity recovery. We evaluate estimation accuracy by averaging over many simulation runs the Kullback-Leibler (KL) distance $$\text{KL} = - \text{logdet}(\boldsymbol{\Sigma}\widehat{\boldsymbol\Omega}) + \text{tr}(\boldsymbol{\Sigma}\widehat{\boldsymbol\Omega}) - p ,$$ where $\boldsymbol{\Sigma}=\boldsymbol\Omega^{-1}$ is the true covariance matrix. Note that the KL distance is zero if the estimated precision matrix equals the true precision matrix.
To evaluate aggregation performance, we use two measures: the Rand index rand1971objective and the adjusted Rand index hubert1985comparing. Both indices measure the degree of similarity between the true partition on the set of nodes ${1,\ldots,p}$ and the estimated partition. The Rand index ranges from zero to one, where one means that both partitions are identical. The adjusted Rand index performs a re-scaling to account for the fact that random chance will cause some variables to occupy the same group.
Finally, to evaluate sparsity recovery, we use the false positive and false negative rates
The FPR reports the fraction of truly zero components of the precision matrix that are estimated as nonzero. The FNR gives the fraction of truly nonzero components of the precision matrix that are estimated as zero.
Data are drawn from a multivariate normal distribution with mean zero and covariance matrix $\boldsymbol\Sigma=\boldsymbol\Omega^{-1}$. We take $p=15$ variables and investigate the effect of increasing the number of variables in Section (ref). We consider four different simulation designs, shown in Figure (ref), each having a different combination of aggregation and sparsity structures for the precision matrix $\boldsymbol\Omega$.
Aggregation is present in the first three structures. The precision matrix has a $G$-block structure with $K=3$ blocks. In Section (ref), we investigate the effect of varying the number of blocks. In the chain graph, adjacent aggregated groups are connected through an edge. This structure corresponds to the motivating example of Section (ref). In the random graph, one non-zero edge in the aggregated network is chosen at random. In the unbalanced graph, the clusters are of unequal size. In the unstructured graph, no aggregation is present.
Across all designs, we take the diagonal elements of $\boldsymbol \Omega$ to be $1$, the elements within a block of aggregated variables to be $0.5$, and the non-zero elements across blocks to be $0.25$. We generate $100$ different data sets for every simulation design and use a sample size of $n=120$. The number of parameters ($p+p(p-1)/2 = 120$) equals the sample size.
The tag-lasso estimator relies on the existence of a tree to perform node dimension reduction. We consider two different tree structures throughout the simulation study. First, we use an “ideal" tree which contains the true aggregation structure as the sole aggregation level between the leaves and the root of the tree. As an example, the true aggregation structure for the chain graph structure is shown in the left panel of Figure (ref). We form ${\bf A}$ corresponding to this oracle tree to obtain the “tag-lasso ideal" estimator.
We also consider a more realistic tree, shown in the right panel of Figure (ref), following a construction similar to that of yan2018rare. The tree is formed by performing hierarchical clustering of $p$ latent points chosen to ensure that the tree contains the true aggregation structure and that these true clusters occur across a variety of depths. In particular, we generate $K$ cluster means $\mu_1, \ldots, \mu_K$ with $\mu_i=1/i$. We set the number of latent points associated with each of the $K$ means equal to the cluster sizes from Figure (ref). These latent points are then drawn independently from $N(\mu_i,[0.05\cdot\min_j(\mu_i - \mu_j)]^2)$. Finally, we form ${\bf A}$ corresponding to this tree to obtain the “tag-lasso realistic" estimator.
We subsequently discuss the results on estimation accuracy, aggregation performance, and sparsity recovery.
\paragraph{Estimation Accuracy.} Boxplots of the KL distances for the three estimators (tag-lasso ideal, tag-lasso realistic and glasso) relative to the oracle are given in Figure (ref). The first three panels correspond to simulation designs with aggregation structures. In these settings, the tag-lasso estimators considerably outperform the glasso, on average by a factor five. The tag-lasso ideal method performs nearly as well as the oracle. Comparing the tag-lasso realistic method to the tag-lasso ideal method suggests a minimal price paid for using a more realistic tree.
The “unstructured” panel of Figure (ref) shows a case in which there is sparsity but no aggregation in the true data generating model. As expected, the glasso performs best in this case; however, we observe minimal cost to applying the tag-lasso approaches (which encompass the glasso as a special case when $\lambda_1 = 0$).
\paragraph{Aggregation Performance.} Table (ref) summarizes the aggregation performance of the three estimators in terms of the Rand index (RI) and adjusted Rand index (ARI). No results on the ARI in the unstructured simulation design are reported since it cannot be computed for a partition consisting of singletons. The tag-lasso estimators perform very well. If one can rely on an oracle tree, the tag-lasso perfectly recovers the aggregation structure, as reflected in the perfect (A)RI values of the tag-lasso ideal method. Even when the tag-lasso uses a more complex tree structure, it recovers the correct aggregation structure in the vast majority of cases. The glasso returns a partition of singletons as it is unable to perform dimension reduction through aggregation, as can be seen from its zero values on the ARI.
\paragraph{Sparsity Recovery.} Table (ref) summarizes the results on sparsity recovery (FPR and FNR). The tag-lasso estimators enjoy favorable FPR and FNR, mostly excluding the irrelevant conditional dependencies (as reflected by their low FPR) and including the relevant conditional dependencies (as reflected by their low FNR). In the simulation designs with aggregation, the glasso pays a big price for not being able to reduce dimensionality through aggregation, leading it to include too many irrelevant conditional dependencies, as reflected through its large FPRs. In the unstructured design, the rates of all estimators are, overall, low.
We investigate the sensitivity of our results to an increasing number of variables $p$. We focus on the chain simulation design from Section (ref) and subsequently double $p$ from 15 to 30, 60 and 120 while keeping the number of blocks $K$ fixed at three. The sample size $n$ is set proportional to the complexity of the model, as measured by $Kp + p$. Hence, the sample sizes corresponding to the increasing values of $p$ are respectively, $n=120,240,480,960$, thereby keeping the ratio of the sample size to the complexity fixed at two. In each setting, the number of parameters to be estimated is large, equal to 120, 465, 1830, 7260, respectively; thus increasing relative to the sample size.
The left panel of Figure (ref) shows the mean KL distance (on a log-scale) of the four estimators as a function of $p$. As the number of nodes increases, the estimation accuracy of the tag-lasso estimators and the oracle increases slightly. For fixed $K$ and increasing $p$, the aggregated nodes---which can be thought of as the average of $p/K$ random variables---may be stabler, thereby explaining why the problem at hand does not get harder when increasing $p$ for the methods with node aggregation. By contrast, the glasso---which is unable to exploit the aggregation structure---performs worse as $p$ increases. For $p=120$, for instance, the tag-lasso estimators outperform the glasso by a factor 50.
Results on aggregation performance and sparsity recovery are presented in Figure (ref) of Appendix (ref). The tag-lasso ideal method perfectly recovers the aggregation structure for all values of $p$. The realistic tag-lasso's aggregation performance is close to perfect and remains relatively stable as $p$ increases. The glasso is unable to detect the aggregation structure, as expected and reflected through its zero ARIs. The tag-lasso estimators also maintain a better balance between the FPR and FNR than the glasso. While their FPRs increase as $p$ increases, their FNRs remain close to perfect, hence all relevant conditional dependencies are recovered. The glasso, in contrast, fails to recover the majority of relevant conditional dependencies when $p=60,120$, thereby explaining its considerable drop in estimation accuracy.
Finally, we investigate the effect of increasing the number of blocks $K$. We take the chain simulation design from Section (ref) and increase the number of blocks from $K=3$ to $K=5,6,10$, while keeping the number of variables fixed at $p=30$. The right panel of Figure (ref) shows the mean KL distance (on a log-scale) of the four estimators as a function of $K$. As one would expect, the difference between the aggregation methods and the glasso decreases as $K$ increases. However, for all $K$ considered, the glasso does far less well than the aggregation based methods.
Similar conclusions hold in terms of aggregation and sparsity recovery performance. Detailed results are presented in Figure (ref) of Appendix (ref). The tag-lasso ideal method performs as well as the oracle in terms of capturing the aggregation structure; the tag-lasso realistic method performs close to perfect and its aggregation performance improves with increasing $K$. In terms of sparsity recovery, the tag-lasso estimators hardly miss relevant conditional dependencies and only include a small number of irrelevant conditional dependencies. The glasso's sparsity recovery performance is overall worse but does improve with increasing $K$.
We demonstrate our method on a financial data set containing daily realized variances of $p=31$ stock market indices from across the world in 2019 ($n=254$). Daily realized variances based on five minute returns are taken from the Oxford-Man Institute of Quantitative Finance (publicly available at http://realized.oxford-man.ox.ac.uk/data/download). Following standard practice, all realized variances are log-transformed. An overview of the stock market indices is provided in Appendix (ref). We encode similarity between the 31 stock market indices according to geographical region, and use the tree shown in Figure (ref) to apply the tag-lasso estimator.
Since the different observations of the consecutive days are (time)-dependent, we first fit the popular and simple {\em heterogeneous autoregressive} (HAR) model of Corsi09 to each of the individual log-transformed realized variance series. Graphical displays of the residual series of these 31 HAR models suggest that almost all autocorrelation in the series is captured. We then apply the tag-lasso to the residual series to learn the conditional dependency structure among stock market indices.
\paragraph{Estimated Graphical Model.} We fit the tag-lasso estimator, with 5-fold cross validation to select tuning parameters, to the full data set, with the matrix ${\bf A}$ encoding the tree structure in Figure (ref). The tag-lasso returns a solution with $K=6$ aggregated blocks; the sparsity pattern of the full estimated precision matrix is shown in the top left panel of Figure (ref). The coloring of the row labels and the numbering of columns convey the memberships of each variable to aggregated blocks (to avoid clutter, only the first column of each block is labeled).
Dimension reduction mainly occurs through node aggregation, as can be seen from the aggregated precision matrix in the bottom right panel of Figure (ref). The resulting aggregated graphical model is rather dense with only about half of the off-diagonal entries being non-zero in the estimated aggregated precision matrix, thereby suggesting strong volatility connectedness. The solution returned by the tag-lasso estimator consists of one single-market block (block 5: Canada) and five multi-market blocks, which vary in size. The Australian, South-America, and all Asian stock markets form one aggregated block (block 6). Note that the tag-lasso has “aggregated" these merely because they have the same {\em non}-dependence structure (i.e.\ all of these markets are estimated to be conditionally independent of each other and all other markets). The remaining aggregated nodes concern the US market (block 4) and three European markets, which are divided into North-Europe (block 1), Central-, South-Europe & STOXX50E (block 2), and West-Europe (block 3). In the aggregated network, the latter two and the US play a central role as they are the most strongly connected nodes: These three nodes are connected to each other, the US node is additionally connected to Canada, whereas these European nodes are additionally connected with North-Europe.
\paragraph{Out-of-sample Performance.} We conduct an out-of-sample exercise to compare the tag-lasso estimator to the glasso estimator. We take a random $n=203$ observations (80% of the full data set) to form a “training sample" covariance matrix and use the remaining data to form a “test sample" covariance matrix ${\bf S}^{\text{test}}$, and repeat this procedure ten times. We fit both the tag-lasso and glasso estimator to the training covariance matrix, with 5-fold cross-validation on the training data to select tuning parameters. Next, we compute their corresponding out-of-sample errors on the test data, as in (ref).
The top right panel of Figure (ref) shows each of these ten test errors for both the tag-lasso (x-axis) and the glasso estimator (y-axis). The fact that in all ten replicates the points are well above the 45-degree line indicates that the tag-lasso estimator has better estimation error than the glasso. Tag-lasso has a lower test error than glasso in all ten replicates, resulting in a substantial reduction in glasso's test errors. This indicates that jointly exploiting edge and node dimension reduction is useful for precision matrix estimation in this context.
We next turn to a data set of gut microbial amplicon data in HIV patients Rivera-Pinto2018, where our goal is to estimate an interpretable graphical model, capturing the interplay between different taxonomic groups of the microbiome. bien2020tree recently showed that tree-based aggregation in a supervised setting leads to parsimonious predictive models. The data set has $n=152$ HIV patients, and we apply the tag-lasso estimator to all $p=104$ bacterial operational taxonomic units (OTUs) that have non-zero counts in over half of the samples. We use the taxonomic tree that arranges the OTUs into natural hierarchical groupings of taxa: with 17 genera, 11 families, five orders, five classes, three phyla, and one kingdom (the root node). We employ a standard data transformation from the field of compositional data analysis (see e.g., aitchison1982statistical) called the centered log-ratio (clr) transformation that is commonly used in microbiome graphical modeling kurtz2015sparse, lo2018pglasso, kurtz2019disentangling. After transformation, kurtz2015sparse apply the glasso, lo2018pglasso incorporate phylogenetic information into glasso's optimization problem through weights within the $\ell_1$-penalty, and kurtz2019disentangling estimate a latent graphical model which combines sparsity with a low-rank structure. We instead, use the tag-lasso to learn a sparse aggregated network from the clr-transformed microbiome compositions. While the clr-transform induces dependence between otherwise independent components, Proposition 1 in cao2019large provides intuition that as long as the underlying graphical model is sparse and $p$ is large, these induced dependencies may have minimal effect on the covariance matrix. Future work could more carefully account for the induced dependence, incorporating ideas from cao2019large or kurtz2019disentangling.
\paragraph{Estimated Graphical Model.} We fit the tag-lasso to the full data set and use 5-fold cross-validation to select the tuning parameters. The tag-lasso estimator provides a sparse aggregated graphical model with $K=28$ aggregated blocks (a substantial reduction in nodes from the original $p=104$ OTUs). The top panel of Figure (ref) shows the sparsity pattern of the $p\times p$ estimated precision matrix (top left) and of the $K \times K$ estimated aggregated precision matrix (top right). A notable feature of the tag-lasso solution is that it returns a wide range of aggregation levels: The aggregated network consists of 17 OTUs, 7 nodes aggregated to the genus level (these nodes start with “g_"), 3 to the family level (these nodes start with“f_"), and 1 node to the kingdom level (this node starts with “k_"). Some aggregated nodes, such as the “g_Blautia" node (block 19), contain all OTUs within their taxa; some other aggregated nodes, indicated with an asterisk like the “k_Bacteria*" node (block 28), have some of their OTUs missing. This latter “block" consists of 18 OTUs from across the phylogenetic tree that are estimated to be conditionally independent with all other OTUs in the data set.
While the tag-lasso determines the aggregation level in a data-driven way through cross validation, practitioners or researchers may also sometimes wish to restrict the number of blocks $K$ to a pre-determined level when such prior knowledge is available or if this is desirable for interpretability. As an illustration, we consider a {\em constrained cross-validation} scheme in which we restrict the number of blocks $K$ to maximally ten and select the sparsity parameters with the best cross validated error among those solutions with $K\leq 10$. The bottom panel of Figure (ref) shows the sparsity pattern of the full and aggregated precision matrices estimated by this constrained version of the tag-lasso.
The resulting network consists of $K=8$ aggregated nodes. The “k_Bacteria*" node now aggregates 78 OTUs that are estimated to be conditionally independent with each other and all others. The interactions among the remaining nodes are shown in the left panel of Figure (ref), which consists of three OTUs (OTU134, OTU156, and OTU161, in pink), three genera (Prevotella, Bacteroides, and Alistipes in orange) and one family (Porphyromonadaceae in blue). The resulting network is much simpler than the one estimated by the glasso, shown in the middle panel of Figure (ref). The glasso finds 58 OTUs to be conditionally independent with all others, but the interactions among the remaining 46 OTUs are much more difficult to interpret. The glasso is limited to working at the OTU-level, which prevents it from providing insights about interactions that span different levels of the taxonomy.
\paragraph{Out-of-sample Performance.} We conduct the same out-of-sample exercise as described in Section (ref). The right panel of Figure (ref) presents the ten test errors (black dots) for the unconstrained CV tag-lasso and glasso. In all but one case, the tag-lasso leads to a better fit than the glasso, suggesting that it is better suited for modeling the conditional dependencies among the OTUs. The unfilled blue dots show the same but for the constrained CV tag-lasso. In all ten cases, it underperforms the unconstrained CV tag-lasso (see shift to the right on the horizontal axis); however, its performance is on a par with the glasso, with test errors close to the 45 degree line. Thus, there does not appear to be a cost in out-of-sample-performance to the interpretability gains of the constrained tag-lasso over the glasso.
Detecting conditional dependencies between variables, as represented in a graphical model, forms a cornerstone of multivariate data analysis. However, graphical models, characterized by a set of nodes and edges, can quickly explode in dimensionality due to ever-increasing fine-grained levels of resolution at which data are measured. In many applications, a tree is available that organizes the measured variables into various meaningful levels of resolution. In this work, we introduce the tag-lasso, a novel estimation procedure for graphical models that curbs this curse of dimensionality through joint node and edge dimension reduction by leveraging this tree as side information. Node dimension reduction is achieved by a penalty that allows nodes to be aggregated according to the tree structure; edge dimension reduction is achieved through a standard sparsity-inducing penalty. As such, the tag-lasso generalizes the popular glasso approach to sparse graphical modelling. An \verb|R| package called \verb|taglasso| implements the proposed method and is available on the GitHub page of the first author.
\paragraph{Acknowledgments} We thank Christian M{\"u}ller for useful discussions. Jacob Bien was supported in part by NSF CAREER Award DMS-1653017 and NIH Grant R01GM123993.