Sparse estimation of huge networks with a block‐wise structure

Summary

1. Introduction

Estimation of large covariance matrices and their inverse has several applications in various areas, from economics and finance to health, biology, computer science and engineering. One important technique developed by the statistical and computer science literature is the graphical modelling approach, which aims at exploring the relationships among a set of random variables through their joint distribution. Under this framework, the Gaussian distribution is often assumed and, in this case, the dependence structure is completely determined by the covariance matrix, or, equivalently, by its inverse, where the off‐diagonal elements are proportional to partial correlations Lauritzen, (1996). Specifically, variables i and j are conditionally independent given all other variables, if and only if the $(i,j)$th element of the inverse covariance matrix, referred to as the precision matrix, is zero. One result in the Gaussian graphical modelling literature is that there is a one‐to‐one correspondence between the joint Gaussian distribution of a vector of random variables and its conditional Gaussian distribution. Under the latter, the distribution of a variable observed in a certain node, given values observed in all other nodes, depends only on the observations in its neighbourhood; see, e.g. Mardia (1988) and Meinshausen and Buhlmann (2006). Hence, the problem of estimating the (inverse) covariance matrix is equivalent to a neighbourhood selection problem. This observation has led to efficient node‐wise LASSO approaches for sparse high‐dimensional graphs; see, e.g. Meinshausen and Buhlmann (2006) and Peng et al. (2009). In contrast to these approaches, Friedman et al. (2008) have developed the Graphical LASSO (GLASSO) approach, where the inverse covariance matrix is directly estimated via penalized likelihood.

Conditional Gaussian models are known in the spatial econometrics literature as conditional autoregressive (CAR) models, representing data from a given spatial location as a function of data in neighbouring locations; see, e.g. Cressie (1993) and Anselin (2010). In a CAR model, the neighbourhood structure is represented by means of the so‐called spatial weights matrix, usually assumed to be known a priori using information on distance between units, such as the geographic, economic, policy or social distance. It is interesting to observe that the problem of estimating the spatial weights matrix in a CAR model is equivalent to a neighbourhood selection problem in a graphical model; for more details, see Section 5. Hence, the spatial weights matrix for CAR models can be estimated by using methods from the Gaussian graphical modelling literature for estimating inverse covariance matrices. While the spatial econometrics literature has been largely immune to the developments in Gaussian graphical modelling, these methods may be useful for a large number of applications in the social sciences.

In this paper, we consider the case of networks with a very large number of nodes and we focus on the estimation of Gaussian graphical models when the dependency between variables has a block‐wise structure. We assume that units can be split into a set of non‐overlapping groups, or blocks, in such a way that the dependence between units only varies across blocks, instead of individual observations. Hence, rather than estimating the links between each pair of units in the sample, we propose to estimate the dependence (links) between groups of cross‐sectional units. Our approach consists of applying the GLASSO methodology of Friedman et al. (2008) to block‐level averages of observations rather than to single observations. When the size of the group is unity, our method collapses to the conventional GLASSO. A major advantage of this method is that its computational cost is greatly reduced and hence it can be adopted for estimation and prediction using very large, or huge, networks. Our approach is also particularly useful when there is a problem of missing values and outliers or when the focus of the analysis is out‐of‐sample prediction.

There exist several examples where it is reasonable to assume a block‐wise dependence structure between units. In economics, preferences for consumer goods of individuals belonging to the same household may react similarly in response to consumption decisions of neighbouring households. Companies belonging to the same sector of economic activity and located within the same geographical area (e.g. the postcode, the region or the country) tend to behave similarly because they have similar characteristics or face similar opportunities and constraints. Thus, it is reasonable to assume that the way they interact with companies from other sectors and/or geographical areas is similar. A block‐wise dependence structure is also a realistic assumption when the variable of interest displays an explicit hierarchical or group membership structure, namely, clustering of units in an organized fashion, such as students within classrooms, members of a household, General Practitioners in a clinic, etc. This is common, for example, when dealing with large, individual‐level, microeconomic or health data sets. Other examples are in neuroscience, where the networks used to represent brain activity have a hierarchical structure, with billions of neurons connected to each other through hub nodes, called voxels, and with connected voxels forming areas that are again connected with each other Luo, (2015). In biology, regulatory networks are thought to have a hub‐type structure, with groups of genes having a similar dependency structure and regulated by a small number of unobserved proteins (Hao et al., 2012). When the grouping is not fully known a priori, we could use methods that allow us to determine endogenously the optimal grouping of cross‐sectional units, by employing techniques from the clustering literature; see, e.g. Lin and Ng (2012), Bonhomme and Manresa (2015) and Ando and Bai (2016).

Exploitation of a priori information on the group structure of variables is not new in the social interaction literature and in the statistical and graphical modelling literature. Empirical works from the social interaction literature typically assume that an individual reacts to the average of others in a predefined group; see Durlauf and Young (2001) and Blume et al. (2013) for a review. Such an assumption implies that the spatial weights matrix has a group‐membership structure, where the weights are identical for all units belonging to the same group, while they are set to zero for the interaction between units belonging to different groups. Lee and Yu (2007) considered the identification and estimation of interaction effects in the context of a spatial autoregressive model where the spatial weights matrix (and the associated precision matrix) has such a block diagonal structure with equal entries. Note that this is a more restrictive assumption to that used in this paper, as it does not allow for dependences between groups. Nevertheless, this model has been widely adopted in several different areas of the social sciences, such as education (Calvó‐Armengol et al., 2009), labour market outcomes (Bayer et al., 2008), crime Sirakaya, (2006) and welfare participation (Bertrand et al., 2000). Similar models have been proposed by the statistical literature, where mixed effect models are commonly used to represent variables with a hierarchical or known group membership structure Goldstein, (2011). When the random effects are assumed to be correlated, these models lead to a covariance matrix that has a block‐wise structure of the same type that we use in this paper, with equal correlation within groups and equal correlation between any two elements of two specified groups Laird and Ware, (1982). Maximum likelihood approaches are typically used for parameter estimation in these models. In the case of a large number of regressors, penalized approaches based on the L₁ penalty are used for estimation and variable selection (Schelldorfer et al., 2014). However, these methods typically require a small number of random effects (blocks).

A number of authors in the literature on graphical modelling have proposed sparse estimation of graphs with a block structure. These methods exploit a priori information on group membership of observations to propose fast, sparse estimation algorithms. Guo et al. (2011) consider a heterogeneous data set where variables, while independent across groups, have a sparse dependency structure within group. The corresponding precision matrix has a block diagonal structure, and the authors propose joint estimation of various blocks by maximizing the corresponding penalized log‐likelihood functions. A similar approach is taken by Mazumder and Hastie (2012), who propose thresholding estimation of a sparse inverse covariance that is a block diagonal matrix of connected components. Wit and Abbruzzo (2015) impose block equality constraints on the parameters of an undirected graphical model to reduce the number of parameters to be estimated. Vinciotti et al. (2016) discuss various forms of block structures for dynamic networks and propose estimation of the associated precision matrix under sparsity and equality constraints on parameters (also known as parameter tying). The inclusion of equality constraints, while reducing the number of parameters, often increases the computational complexity of the estimation procedures. For example, the general block structures considered by Wit and Abbruzzo (2015) and Vinciotti et al. (2016) imply a computational cost of the estimation procedure that is higher compared to the approaches of Guo et al. (2011) and Mazumder and Hastie (2012), where the assumed block structure allows the large GLASSO problem to be split into many, smaller tractable problems.

In this paper, we use block structures with the intent to achieve computational efficiency, allowing us to infer networks of very large dimensions. Differently from Guo et al. (2011) and Mazumder and Hastie (2012), our approach does not need to impose block‐diagonality of the precision matrix. However, we assume that units can be split into groups in such a way that the covariance (and associated precision matrix) only varies across blocks, rather than individual observations.

The rest of the paper is structured as follows. In Section 2, we describe the main features of our graphical model with block‐wise dependence structure, while in Section 3 we propose our estimator based on GLASSO. In Section 4, we run Monte Carlo experiments to investigate the small‐sample properties of the proposed estimator. In Section 5, we carry out an empirical study on the economic growth of a set of small regions in Europe. Finally, in Section 6, we provide some concluding remarks. The Appendix provides the proofs.

We use $|{\lambda }_1(\mathbf{A})|\ge |{\lambda }_2(\mathbf{A})|\ge \cdots \ge |{\lambda }_n(\mathbf{A})|$ to denote the eigenvalues of a matrix $\mathbf{A}\in {\mathbb{M}}^{n\times n}$⁠, where ${\mathbb{M}}^{n\times n}$ is the space of $n\times n$ matrices. Tr$(\mathbf{A})$ is the trace of $\mathbf{A}\in {\mathbb{M}}^{n\times n}$⁠, while its Frobenius norm is ${\displaystyle {\parallel \mathbf{A}\parallel }_F={(\sum _{i,j=1}^m{a}_{ij}^2)}^{1/2}}$⁠. K is used for a fixed positive constant that does not depend on N; $S^c$ is used to denote the complement of a set S.

2. Block‐Wise Dependence Structure in Huge Networks

3. Block‐Glasso Approach

• Step 1. Estimate ${\mathbf{\Phi }}_G={\mathbf{\Psi }}_G^{-1}$ by applying the GLASSO to ${\overline{\mathbf{y}}}_{G,t}$⁠, $t=1,2,\dots ,T$⁠. This allows us to obtain ${\widehat{\sigma }}_{gh}$ for $g\ne h=1,2,\dots ,G$⁠, and ${\widehat{\psi }}_{\mathit{gg}},g=1,2,\dots ,G$⁠.

• Step 2. Estimate ${\gamma }_g$ by exploiting identity 2.4 and 3.4. Noting that $E[(1/\mathit{MT})$${\displaystyle \sum _{i\in g}\sum _{t=1}^T{y}_{\mathit{it}}^2]={\sigma }_{\mathit{gg}}+{\gamma }_g}$⁠, while ${\displaystyle E[(1/\mathit{MT})\sum _{i\in g}\sum _{t=1}^T{\overline{y}}_{gt}^2]={\sigma }_{\mathit{gg}}+(1/M){\gamma }_g}$⁠, we can consider the following estimator for ${\widehat{\gamma }}_g$⁠:

  ${\displaystyle {\widehat{\gamma }}_g=\frac{M}{M-1}(\frac{1}{\mathit{MT}}\sum _{i\in g}\sum _{t=1}^T{y}_{\mathit{it}}^2-{\widehat{\psi }}_{\mathit{gg}}),g=1,2,\dots ,G.}$

  (3.5)

  Hence, use 2.6 to recover $\widehat{\mathbf{\Theta }}$⁠:

  $\widehat{\mathbf{\Theta }}=(\frac{1}{M}{\widehat{\mathbf{\Phi }}}_G\otimes \frac{1}{M}{\mathbf{1}}_M)+({\widehat{\mathbf{\Gamma }}}_G^{-1}\otimes ({\mathbf{I}}_M-\frac{1}{M}{\mathbf{1}}_M)).$

  (3.6)

The following theorems derive the asymptotic properties of estimator 13 when both N and T go to infinity.  

Theorem 3.2 is a straightforward consequence of the sparsistency theorem of Lam and Fan (2009) applied to ${\widehat{\mathbf{\Phi }}}_G$⁠; see also Rothman et al. (2008) and Guo et al. (2011).

Hence, for $\widehat{\mathbf{\Theta }}$ to be a good proxy of $\mathbf{\Theta }$⁠, G needs to be small (or, equivalently, M large) and ${\mathbf{\Phi }}_G$ needs to be a sparse matrix, as measured by $s_G$⁠. Note, however, that, from 2.6, the off‐diagonal elements of $\mathbf{\Theta }$ are proportional to $1/M^2$⁠. Hence, for fixed G, as M increases the (relative) effect of each individual neighbour on each unit would disappear and, in the limit, the precision matrix would become a diagonal matrix. A similar result has been obtained by Lee (2002) in the context of a spatial autoregressive (SAR) model where each spatial unit is influenced aggregately by a significant portion of other spatial units in the sample. Lee (2002) showed that if each spatial unit in the limit has infinitely many neighbours (which would happen in our case for G fixed and M increasing), then the ordinary least‐squares (OLS) estimator for a SAR model would still be consistent and even asymptotically efficient. In Section 4, we investigate the properties of our estimator for different values of G relative to N.

A major advantage of our proposed estimation procedure is that it is considerably faster than the conventional GLASSO for estimating an $N\times N$ precision matrix. Using the algorithm proposed by Friedman et al. (2008), the computational cost associated with a coordinate descendent update would decrease from $O(N^2)$ to $O(G^2)$⁠. This could decrease further to $O(G)$ using faster algorithms, such as QUIC (Hsieh et al., 2014). Another advantage of our approach is that using block averages rather than single observations greatly helps in the presence of missing values, a common problem in statistical analysis. Exploiting group membership information is also very useful for prediction purposes on a hold‐out sample of units, for which the position in the (individual‐level) network is usually unknown. It is important, however, to remark that our approach requires a priori information on the block structure. If this is not available, then one could exploit methods from the clustering literature that allow us to determine endogenously the optimal grouping of cross‐sectional units, such as the k‐means algorithm Forgy, (1965) extended to allow for covariates in the model; see, in particular, Lin and Ng (2012) and Bonhomme and Manresa (2015), and also Ando and Bai (2016). Our approach also has potential application in the area of spatial econometrics. Given the equivalence between CAR models and the joint Gaussian distribution emphasized by many authors – see, among others, Mardia (1988) and Meinshausen and Buhlmann (2006) – this method provides a means for estimating spatial weights matrices in the context of very large panel data. Later in the paper, we offer a small empirical exercise using CAR models.

Finally, it is important to remark that our approach does not allow us to estimate consistently the precision matrix when this arises from one or more common, pervasive factors. Unobserved common factors occur in time series as a result of global shocks, namely unexpected events that may hit all statistical units, although with different intensities Stock and Watson, (2010). These large‐scale perturbations affect micro‐level population units and are often responsible for observable co‐movements of a large number of time series. We observe that our model is more parsimonious than the common factor specification and may be useful in situations where T is too short to allow for fully unrestricted common effects. However, in a large T setting, in the presence of unobserved common factors, our approach can be applied to de‐factored residuals, after estimating common factors using methods such as principal components Bai, (2003) or the Common Correlated Effects methodology Pesaran, (2006).

3.1. Case of Blocks With Unequal Size

3.2. Allowing for General Intra‐Block Correlation Structure

We can show that, under the condition that ${\displaystyle (1/M)\sum _{k\in g}{\pi }_{\mathit{ik}}\approx 0}$ for all i, the Π matrix can be estimated by the covariance matrix of $y_{\mathit{it}}-{\overline{y}}_{\mathit{gt}}$ for each block. This results in a relatively easy implementation, whereby we first calculate ${\overline{y}}_{gt}$ and apply the block‐GLASSO outlined in Section 3 to compute ${\widehat{\mathbf{\Sigma }}}_G$⁠. Hence, we calculate the deviations of each value $y_{\mathit{it}}$ from its corresponding group‐level average, namely $y_{\mathit{it}}-{\overline{y}}_{\mathit{gt}}$⁠, and we apply the conventional GLASSO to all $y_{\mathit{it}}-{\overline{y}}_{gt}$ for each block, separately. This approach requires that ${\pi }_{\mathit{ij}}$⁠, namely the deviations of ${\sigma }_{\mathit{ij}}$ from ${\sigma }_{\mathit{gg}}$⁠, are not too large, so that the ${\overline{y}}_{\mathit{gt}}$ can be used to consistently estimate ${\sigma }_{\mathit{gg}}$⁠. The computational complexity of this procedure rises to $O(G^2)+O({\mathit{GM}}^2)$⁠, as it is necessary to estimate G blocks of size M. In the rest of the paper, we refer to this approach as the flexible block‐GLASSO.

4. Monte Carlo Experiments

In each experiment, we compute the block‐GLASSO and the conventional GLASSO, for all pairs of N and T with $N=50$ and 100 and $T=10$⁠, 50 and 200. As for the choice of G, we try $G=N/2$⁠, and $N/5$⁠. Each experiment is replicated $R=250$ times. We also carry out another set of experiments with N much larger than T, and set $N=500$⁠, 1,000 and 2,000 and $T=20$⁠. In this set of experiments, given the computational difficulties and poor performance in computing conventional GLASSO for such large networks, we only provide results for the block‐GLASSO. Under the dynamic set‐up 4.3, we only run experiments for large T (i.e. $T=50$ and 200) to avoid incurring bias of the OLS estimator for short panels.²

4.1. Results

The results are summarized in Tables 1–6 and Figures 1–2. The results from Table 1 show that, when data have block‐wise dependence structure, our method greatly outperforms the conventional GLASSO for all combinations of N, T and G. In particular, the F1 score and AUC show that block‐GLASSO has higher true positive rates and substantially lower false positive rates, while the EL and FL are always lower for block‐GLASSO, indicating that the latter provides a better estimation of the precision matrix. However, it is interesting to note that when $T=10$ and $G=N/2$ the block‐GLASSO does not perform well relative to other cases, and its properties are much worse than the case $T=10$ and $G=N/5$⁠. More generally, Tables 1 and 2 show that for the same pair of N and T, the properties of block‐GLASSO deteriorate as G rises, thus confirming our theoretical results that, holding N and T fixed, the estimation error is higher when G is large or, equivalently, M small. This result is also confirmed by Figure 1, showing the ROC curves for the block‐GLASSO for varying N, T and G. As expected, the performance of the estimator improves as N increases (and hence M) for fixed T and G, and as T increases for fixed N and G, while it deteriorates as G rises, holding N and T constant.

Figure 1.

Open in new tabDownload slide

Block‐GLASSO ROC curves: varying values of (a) N, (b) T and (c) G. [Color figure can be viewed at wileyonlinelibrary.com]

Figure 2.

Open in new tabDownload slide

Flexible block‐GLASSO, group LASSO and conventional GLASSO: within‐block variation, $N=100$⁠, $T=200$⁠. [Color figure can be viewed at wileyonlinelibrary.com]

Table 3 reports the small‐sample properties of OLS and GLS estimators as well as of the block‐GLASSO. As expected in the case of cross‐sectionally correlated regression errors, the OLS estimator, while having a bias comparable to that of the GLS, has higher RMSE and is oversized for all combinations of N, T and G. Hence, ignoring the network leads to severe over‐rejection of the null hypothesis. Looking at the GLS estimator, its empirical size is close to the nominal size of 5% in most cases, although some size distortions can be observed when $T=10$ and $G=N/2$⁠, namely, for short panels characterized by the presence of many, small groups. In fact, under this case the block‐GLASSO does not perform well, having small F1 and AUC and large EL and FL, thus confirming our asymptotic results reported in Section 3. Similar results can be observed in Table 4 for the case where the dependent variable is generated by the first‐order autoregressive model 4.3. Under non‐normal errors (Table 5), the block‐GLASSO still performs well in detecting the network, as confirmed by F1 and AUC values similar to those reported in Table 1, although its EL and FL are much higher than in the normal counterpart.

Table 6 shows results when the error covariance matrix displays general intra‐block variation (see 4.4 and 4.5). It is interesting to observe that the empirical size of the GLS estimator of β when ignoring the intra‐block variation (block‐GLASSO) is in some cases still close to the nominal value of 5%. The GLS estimator based on the more general procedure (flexible block‐GLASSO) shows a good performance only for smaller values of G, perhaps because under small G (and hence large M) the covariance of ${\overline{y}}_{gt}$ better approximates the part of the covariance that is block‐wise. We also remark that the more flexible procedure is computationally much slower than the block‐GLASSO. Figure 2 shows the ROC for the flexible block‐GLASSO and the conventional GLASSO, as well as the group LASSO by Yuan and Lin (2006). The use of a group penalty in the group LASSO encourages the recovery of the block structure, although it does not impose it as in the block‐GLASSO. Because the group LASSO has been developed in the context of regression analysis, we apply it to our model as a neighbourhood selection problem for each node of the network. It is interesting to see from Figure 2 that the group LASSO approach performs less well than the block‐GLASSO, but slightly better than the conventional GLASSO, as the latter does not use any a priori information about the blocks.

5. An Empirical Example: Spatial Spillovers in Regional Growth and Convergence in Europe

We use block‐GLASSO for estimating a growth equation in per‐capita gross value‐added and for testing for economic convergence of European regions. The debate on whether there exists convergence in per‐capita input and income across nations is still open, with results obtained that differ depending on the sample period and the regions included, as well as the estimation methods adopted. A number of authors have highlighted the importance of incorporating spatial effects when studying economic growth and regional convergence and have proposed the use of spatial econometric techniques; see, among others, Rey and Montouri (1999), Ertur and Koch (2007) and Cuaresma and Feldkircher (2013). Spatial dependence in regional economic growth is likely to arise from technology spillover across neighbouring regions and from factor mobility, as well as from the presence of spatial heterogeneity Rey and Montouri, (1999). In the presence of spatial dependence in economic growth data, if ignored, estimates of the speed of income convergence across geographical regions will be biased.

We contribute to this literature by estimating a growth equation with spatial spillovers and we use the block‐GLASSO procedure to estimate the spatial weights matrix. We use data on gross value‐added per worker (GVA) for 1,088 NUTS3 observed over the period 1980–2012 in 14 European countries.³ The NUTS classification is a hierarchical system for dividing up the economic territory of the European Union (EU) for the purpose of socio‐economic analysis of the regions and design of EU regional policies. It subdivides the EU territory into regions at the three different levels, NUTS1, NUTS2 and NUTS3, moving from larger to smaller geographical units.

In 5.1, $w_{\mathit{ij}}$ is the (⁠$i,j$⁠)th element of an $N\times N$ matrix, $\mathbf{W}$⁠, known as the spatial weights matrix, such that $w_{ii}=0$⁠. In spatial econometrics, $\mathbf{W}$ is often assumed to be known using a priori information (e.g. from economic theory) on how statistical units potentially interact. Spatial weights based on geographical or travel distance, or contiguity have been used for modelling spatial spillovers in the economic growth equation, although this has been pointed out as being unrealistic Cuaresma and Feldkircher, (2013). In this application, we keep $\mathbf{W}$ as unknown and estimate it using our block‐GLASSO approach. While the unit of analysis is the NUTS3 region, we take as groups larger geographical areas, given by 80 NUTS1 and then 211 NUTS2 European regions. Other grouping criteria may undoubtedly be suggested, for example by looking at the literature on club convergence – see, among others, Corrado et al. (2005) – or using methods for identifying communities in social networks from the graph modelling literature Freeman, (1979).