Model selection and estimation in the Gaussian graphical model

Abstract

1 Introduction

Let X = (X⁽¹⁾, …, X^(p)) be a p-dimensional random vector following a multivariate normal distribution with unknown mean μ and nonsingular covariance matrix Σ. Given a random sample X₁, …, X_n of X, we wish to estimate the concentration matrix C = Σ⁻¹. Of particular interest is the identification of zero entries in the concentration matrix C = (c_ij), since a zero entry c_ij = 0 indicates the conditional independence between the two random variables X⁽ⁱ⁾ and X^(j) given all other variables. This is the covariance selection problem (Dempster, 1972) or the model-selection problem in the Gaussian concentration graph model (Cox & Wermuth, 1996).

A Gaussian concentration graph model for the Gaussian random vector X is represented by an undirected graph G = (V, E), where V contains p vertices corresponding to the p coordinates and the edges describe the conditional independence relationships among X⁽¹⁾, …, X^(p). The edge between X⁽ⁱ⁾ and X^(j) is absent if and only if X⁽ⁱ⁾ and X^(j) are independent conditional on the other variables, and corresponds to c_ij = 0. Thus parameter estimation and model selection in the Gaussian concentration graph model are equivalent to estimating parameters and identifying zeros in the concentration matrix; see Whittaker (1990), Lauritzen (1996) and Edwards (2000) for statistical properties of Gaussian concentration graph models and commonly used model selection and parameter estimation methods in such models.

The standard approach to model selection in Gaussian graphical models is greedy stepwise forward-selection or backward-deletion, and parameter estimation is based on the selected model. In each step the edge selection or deletion is typically done through hypothesis testing at some level α. It has long been recognized that this procedure does not correctly take account of the multiple comparisons involved (Edwards, 2000). Another drawback of the common stepwise procedure is its computational complexity. To remedy these problems, Drton & Perlman (2004) proposed a method that produces conservative simultaneous 1 − α confidence intervals, and uses these confidence intervals to do model selection in a single step. The method is based on asymptotic considerations. Meinshausen & Bühlmann (2006) proposed a computationally attractive method for covariance selection that can be used for very large Gaussian graphs. They perform neighbourhood selection for each node in the graph and combine the results to learn the structure of a Gaussian concentration graph model. They showed that their method is consistent for sparse high-dimensional graphs. In all of the above mentioned methods, model selection and parameter estimation are done separately. The parameters in the concentration matrix are typically estimated based on the model selected. As demonstrated by Breiman (1996), the discrete nature of such procedures often leads to instability of the estimator: small changes in the data may result in very different estimates. Other recent advances include a Duke University discussion paper by A. Dobra and M. West, who presented a novel Bayesian framework for building Gaussian graphical models and illustrated their approach in a large scale gene expression study, and Li & Gui (2006), who adopted gradient-directed regularization, which is described in a technical report by J. Friedman and B. Popescu, available at http://www-stat.stanford.edu/∼jhf, to construct sparse Gaussian graphical models.

In this paper, we propose a penalized-likelihood method that does model selection and parameter estimation simultaneously in the Gaussian concentration graph model. We employ an ℓ₁ penalty on the off-diagonal elements of the concentration matrix. This is similar to the idea of the lasso in linear regression (Tibshirani, 1996). The ℓ₁ penalty encourages sparsity and at the same time gives shrinkage estimates. In addition, we explicitly ensure that the estimator of the concentration matrix is positive definite. We also introduce a ‘nonnegative garrote’ type method that is closely related to the aforementioned approach.

There is a connection between the neighbourhood-selection method in Meinshausen & Bühlmann (2006) and our penalized-likelihood approach, which we illustrate in 5. The neighbourhood-selection method can be cast as a penalized M-estimation without incorporating the positive definiteness or symmetry constraint. The loss function in the penalized M-estimation is a particular quadratic form. The neighbourhood selection method is computationally faster because of its simpler form and because it does not consider the positive definite constraint. Our method is more efficient because of the incorporation of the positive definite constraint and the use of likelihood.

Throughout the paper we assume that the observations are suitably centred and scaled. The sample mean is centred to be zero. One may scale to have the diagonal elements of the sample covariance matrix equal to one or to have the diagonal elements of the sample concentration matrix equal to one. In our experience these two scalings give very similar performance, and in this paper we assume the latter since it seems to be more natural for estimating the concentration matrix.

2 Methodology

2·1 Lasso-type estimator

The loglikelihood for μ and C = Σ⁻¹ based on a random sample X₁, …, X_n of X is

up to a constant not depending on μ and C. The maximum likelihood estimator of (μ, Σ) is (X, A), where

The commonly used sample covariance matrix is S = nA/(n − 1). The concentration matrix C can be naturally estimated by A⁻¹ or S⁻¹. However, because of the large number, p(p + 1)/2, of unknown parameters to be estimated, S is not a stable estimator of Σ for moderate or large p. In general, the matrix S⁻¹ is positive definite when ⁠, but does not lead to ‘sparse’ graph structure since the matrix typically contains no zero entry.

To achieve ‘sparse’ graph structure and to give a better estimator of the concentration matrix, we adapt the lasso idea and seek the minimizer of

1

over the set of positive definite matrices. Here is a tuning parameter. When t = ∞, the solution to (1) is the maximum likelihood estimator A⁻¹ provided that the inverse exists. On the other hand, if t = 0, then the constraint forces C to be diagonal, which implies that X⁽¹⁾, …, X^(p) are mutually independent. It is clear that regardless of t. Since the observations are centred, we have ⁠. Therefore, Ĉ is the positive definite matrix that minimizes

2

We can further rewrite (2) as

3

Since both the objective function and feasible region of (3) are convex, we can equivalently use the Lagrangian form

4

with being the tuning parameter.

2·2 Nonnegative garrote-type estimator

If a relatively reliable estimator of C is available, a shrinkage estimator can be defined through ⁠, where the symmetric matrix D = (d_ij) is the minimizer of

5

and with C positive definite. Equivalently, this can be written as

6

subject to and with C positive definite. For a relatively large sample size, A⁻¹ is an obvious choice for the preliminary estimator. This procedure is similar in spirit to the nonnegative garrote estimator proposed by Breiman (1995) for linear regression.

2·3 Illustration

Consider the special case in which p = 2. Denote the maximum likelihood estimator of the concentration matrix by

where the diagonal elements are 1 because of the scaling. Therefore,

Substitution in (4) gives

where we used the fact that C is symmetric.

LEMMA 1.

In the case of the bivariate normal, the proposed penalized likelihood estimator given by the solution to (4) is

where (x)₊ = max(x, 0) and

7

Similarly, the garrote type estimator can also be obtained in an explicit form in this case.

LEMMA 2.

With⁠, the minimizer of (6) is

and ĉ₁₁ = ĉ₂₂are given by (7).

The estimators are illustrated Fig. 1. If the true value of c₁₂ is zero, r will tend to be small in magnitude. With an appropriately chosen λ, both estimates of c₁₂ can be shrunk to zero, so that model selection for the graphical model can be achieved. Note that the proposed estimators are continuous functions of r and consequently of the data. Such continuity, not shared by the existing methods that perform maximum likelihood estimation on a selected graph structure, ensures the stability of our estimators. The garrote-type estimator penalizes large r’s less heavily than small r’s. As will be demonstrated in the next section, this can be advantageous for model-fitting. However, the disadvantage of the garrote-type estimator is that it can only be applied when a good initial estimator is available.

Fig. 1

Open in new tabDownload slide

(a) Lasso-type estimator, (b) garrote type estimator for the case of p = 2.

3 Asymptotic Theory

In this section, we derive asymptotic properties of the proposed estimators that are analogous to those for the lasso (Knight & Fu, 2000). For simplicity, we assume that p is held fixed as the sample size ⁠. Although it might be more realistic to consider the case when as ⁠, the following results nevertheless provide an adequate illustration of the mechanism of the proposed estimators.

THEOREM 1.

Ifas⁠, the lasso-type estimator is such that

in distribution, where

in which W is a random symmetric p × p matrix such that⁠, and Λ is such that

As an illustration, consider an example where p = 3 and

Note that, for i ≠ j ≠ k ≠ l,

After some tedious algebraic manipulation, we obtain that

We simulated 1000 observations from the distribution of arg min V. Figure 2 gives the scatterplot of the off-diagonal elements for λ₀ = 0, 0·5 and 1. When λ₀ = 0, our estimator is asymptotically equivalent to the maximum likelihood estimator, and the asymptotic distribution for the elements of C is multivariate normal; see Fig. 2(a). If λ₀ > 0, the proposed estimator will have a positive probability of estimating c₁₃ by its true value 0, and this probability increases as λ₀ increases. From Theorem 1 pr(ĉ₁₃ = 0) tends to 0·30 if λ₀ = 0·5 and to 0·45 when λ₀ = 1.

Fig. 2

Open in new tabDownload slide

Example with p = 3. Asymptotic distribution of our estimator as estimated by 1000 simulations, for different values of λ₀, (a) λ₀ = 0, (b) λ₀ = 0·5, (c) λ₀ = 1.

Similarly to Theorem 1, we can derive the asymptotic properties of the nonnegative garrote-type estimator.

Theorem 2 indicates that the garrote-type estimator enjoys the so-called oracle property: it selects the right graph with probability tending to one and at the same time gives a root-n consistent estimator of the concentration matrix.

4 Computation

4·1 The maxdet problem

The nonlinearity of the objective function and the positive-definiteness constraint make the optimization problem (3) nontrivial. We take advantage of the connection between (3) and the determinant-maximization problem, the maxdet problem (Vandenberghe et al., 1998), which can be solved very efficiently with the interior point algorithm.

The maxdet problem is of the form

where b ∈ R^m, G(x) is positive definite, F(x) is positive semidefinite, and the functions and are affine:

where F_i and G_i are symmetric matrices. To use the algorithm of Vandenberghe et al. (1998), it is also required that F_i, i = 1, …, m, be linearly independent and the same be true of G_i, i = 1, …, m. It is not hard to see that the garrote-type estimator (5) solves a maxdet problem.

4·2 Algorithm for lasso-type estimator

If the signs of the c_ij’s are known, (3) can be expressed as the following maxdet problem:

subject to ∑_ic_iiI⁽ⁱ⁾ + ∑_i<jc_ijI^(ij) being positive definite,

8

where C = (c_ij), S = (s_ij), A = (a_ij), I⁽ⁱ⁾ is an n × n matrix with the (i, i)th entry being 1 and all other entries being 0, I^(ij) is an n × n matrix with the (i, j)th and the (j, i)th entries being 1 and all other entries being 0, and s_ij is the sign of c_ij. Since the signs of the c_ij’s are not known in advance, we propose to update the s_ij’s and c_ij’s iteratively using the following steps. In our experience the algorithm usually converges within a small number of iterations. Clearly, other initial values for s can also be used.

• Step 1. Let and for all i ≠ j.

• Step 2. Let Ĉ_new solve (8) over the set of positive definite matrices.

• Step 3. If Ĉ_new = Ĉ_old, then stop and let Ĉ = Ĉ_new. Otherwise, set Ĉ_old = Ĉ_new and s_ij = −s_ij for all pairs (i, j) such that ĉ_ij = 0 and go back to Step 2.

4·3 Tuning

So far we have concentrated on the calculation of the minimizer of (3) for any fixed tuning parameter t. In practice, we need to choose a tuning parameter so as to minimize a score measuring the goodness-of-fit. A commonly used such score is the multifold crossvalidation score, but a computationally more efficient alternative is the BIC for model selection and estimation. To evaluate the BIC for the current setting, one must first obtain an estimate of the degrees of freedom, which is defined as the number of unknown parameters in the case of the maximum likelihood estimate.

Let  = Ĉ⁻¹ and Ω = {(i, j):ĉ_ij ≠ 0}. From the Karush-Kuhn-Tucker conditions, it is not hard to see that the lasso-type estimator satisfies â_ij = a_ij + λsign(ĉ_ij) for all pairs (i, j) ∈ Ω. The remaining card(Ω^c)/2 unique entries of  can be obtained by solving card(Ω^c)/2 equations, Ĉ_ij = 0, where i < j and card(·) is the cardinality of a set. Therefore, Ĉ relies on the observations only through a_ij, (i, j) ∈ Ω. Note that the number of parameters in {a_ij:(i, j) ∈ Ω} is ⁠. Since the a_ij’s are maximum likelihood estimates, we can define, for a given tuning parameter t,

where ê_ij = 0 if ĉ_ij = 0, and ê_ij = 1 otherwise.

5 Quadratic Approximation

Provided that A is nonsingular, a second-order approximation to the objective function of (3) around A⁻¹ can be written as (Boyd & Vandenberghe, 2003)

Therefore, the solution to (4) can be approximated by the solution to

9

This second-order approximation is closely connected to the approach proposed by Meinshausen & Bühlmann (2006). In their approach, for each i = 1, …, p, we seek the minimizer of

10

where X^[−i] is the n × (p − 1) matrix resulting from the deletion of the ith column from the data matrix X. A vertex j is taken to be a neighbour of vertex i if and only if ⁠. The two vertices are connected by an edge in the graphical model if either vertex is the neighbour of the other one.

Note that θ_ii, i = 1, …, p, are not determined. For notational purposes, we write θ_ii = 1 for i = 1, …, p. Recall that we scale each component of X so that all the diagonal elements of the sample concentration matrix are unity. The following lemma reveals a close connection between the approach of Meinshausen & Bühlmann (2006) and the second-order approximation (9).

LEMMA 4.

The matrix Θ = (θ_ij) defined by (10) is the unconstrained solution to

11

over all p × p matrices with diagonal elements fixed at 1.

Lemma 4 shows that the approach of Meinshausen & Bühlmann (2006) seeks a sparse C close to A⁻¹. However, it does not incorporate the symmetry and positive-definiteness constraint in the estimation of the concentration matrix, and therefore an additional step is needed to estimate either the covariance matrix or the concentration matrix. Also, the loss function used by Meinshausen & Bühlmann is different from the quadratic approximation to the loglikelihood, and therefore the approach is expected to be less efficient than our penalized likelihood method or the corresponding quadratic approximation (9).

6 Simulation

We consider eight different models in our simulation.

• Model 1. Heterogeneous model with Σ = diag(1, 2, …, n).

• Model 2. An AR(1) model with c_ii = 1 and c_{i, i−1} = c_{i−1, i} = 0·5.

• Model 3. An AR(2) model with c_ii = 1, c_{i, i−1} = c_{i−1, i} = 0·5 and c_{i, i−2} = c_{i−2, i} = 0·25.

• Model 4. An AR(3) model with c_ii = 1, c_{i, i−1} = c_{i−1, i} = 0·4 and c_{i, i−2} = c_{i−2, i} = c_{i, i−3} = c_{i−3, i} = 0·2.

• Model 5. An AR(4) model with c_ii = 1, c_{i, i−1} = c_{i−1, i} = 0·4, c_{i, i−2} = c_{i−2, i} = c_{i, i−3} = c_{i−3, i} = 0·2 and c_{i, i−4} = c_{i−4, i} = 0·1.

• Model 6. Full model with c_ij = 2 if i = j and c_ij = 1 otherwise.

• Model 7. Star model with every node connected to the first node, with c_ii = 1, c_{1, i} = c_{i, 1} = 0·2 and c_ij = 0 otherwise.

• Model 8. Circle model with c_ii = 1, c_{i, i−1} = c_{i−1, i} = 0·5 and c_1n = c_n1 = 0·4.

For each model, we simulated samples with size 25 and dimension p = 5, or size 50 and dimension 10. We compare our methods with the approach of Meinshausen & Bühlmann (2006) and the method proposed by Drton & Perlman (2004) in terms of the Kullback–Leibler loss,

the number of false positives (FP; incorrectly identified edges) and the number of false negatives (FN; incorrectly missed edges). The approach of Meinshausen & Bühlmann (2006) was implemented using the LARS package from R and the method of Drton & Perlman (2004) has also been implemented in the SIN package of R. Their method gives each edge of the full graph a p-value and two different cut-off values, 5% and 25%, were suggested in their original paper. Both of these methods focus on model selection and do not consider the problem of estimating the covariance matrix or the concentration matrix. For comparison, we estimate the concentration matrix by the maximum likelihood estimate after the graph structure is selected using their methods. Table 1 documents the means and standard errors, in parentheses, from 100 runs for each combination. Our penalized likelihood method is referred to as Lasso in the table because of its connection to the idea of the lasso in linear regression. Similarly, the extension described in 2·2 is referred to as Garrote in the table.

As shown in Table 1, the proposed penalized likelihood methods enjoy better performance than the other methods. The method of Meinshausen & Bühlmann (2006) and both versions of the Drton–Perlman method tend to have larger FN, which may partly explain their relatively poor performance. However, the results suggest that the proposed penalized likelihood approach combined with BIC may have relatively larger FP. The solution path of (3) may therefore be more informative in determining the graph structure.

All four methods under comparison require the selection of tuning parameters that control the trade-off between sensitivity and specificity. To appreciate better the merits of different methods independently of the tuning parameter selection, we plotted the receiver operating characteristic curves for different models and methods, each averaged over the 100 simulated datasets. The AR(4) and full models are not included because the specificity in these cases is not well defined when p = 5. From the plot, not shown here because of space restrictions, the proposed methods outperform the other approaches for all models when p = 5 and n = 25. In the cases when p = 10 and n = 50, the performance of all methods improves but Lasso and Garrote still enjoy competitive performance when compared with the other approaches.

7 Real World Examples

We first consider three real-world examples. The cork borings data are presented in Whittaker (1990, p. 267) and were originally used by Rao (1948). The p = 4 measurements are the weights of cork borings on n = 28 trees in the four directions, north, east, south and west. Fret's heads dataset contains head measurements on the first and the second adult son in a sample of 25 families. The 4 variables are the head length of the first son, the head breadth of the first son, the head length of the second son and the head breadth of the second son. The data are also presented in Whittaker (1990, p. 255). The Mathematics marks dataset (Mardia et al., 1979, p. 3) contains the marks of n = 88 students in the p = 5 examinations in mechanics, vectors, algebra, analysis and statistics. The data also appear in Whittaker (1990, p. 1).

Figures 3–5 depict the solution paths of (3) for each of the three datasets.

Fig. 3

Open in new tabDownload slide

Cork borings dataset. The Meinshausen–Bühlmann method selects (d); both methods based on SIN, with cut-off 0·05 and 0·25, select (b); both Lasso and Garrote with BIC select (e).

Fig. 4

Open in new tabDownload slide

Fret's heads dataset. The Meinshausen–Bühlmann method selects (f); SIN with cut-off 0·05 selects (a); SIN with cut-off 0·25 selects (b); both Lasso and Garrote with BIC select (f).

Fig. 5

Open in new tabDownload slide

Mathematics marks dataset. The Meinshausen–Bühlmann method selects (f); SIN with cut-off 0·05 selects (b) with an additional edge between Mechanics and Vectors; SIN with cut-off 0·25 selects (f); Lasso with BIC selects (j); and Garrote with BIC selects (f).

To compare the accuracy of different methods, fivefold crossvalidation was applied on the datasets. Table 2 documents the average values of KL distances for each method, where now

in which Ĉ is the concentration matrix estimated on the training set and is the sample covariance matrix evaluated on the test set.

Next we considered a larger problem. The opening prices of 35 stocks were collected for the years 2003 and 2004. Different methods were applied to estimate the covariance matrix using the data from 2003. The KL loss of the estimates are then evaluated using the data from 2004, and Table 3 reports the improved KL loss over the sample covariance matrix.

As shown in Tables 2 and 3, the proposed penalized likelihood methods enjoy very competitive performance.

The authors wish to thank the editor and two anonymous referees for comments that greatly improved the paper. This research was supported in part by grants from the U.S. National Science Foundation.

Appendix

Proofs

Proof of Lemma 1

The first-order condition leads to

A1

A2

A3

where sign(c₁₂) = 1 if c₁₂ > 0, sign(c₁₂) = −1 if c₁₂ < −1, and sign(c₁₂) is anywhere between −1 and 1 if c₁₂ = 0. Equation (7) can be easily obtained from (A1) and (A2). Together with (A3), we conclude that

A4

The sign of the left-hand side of the above equation, sign(c₁₂), should therefore be equal to the sign of the right-hand side, sign(r). It follows that (A4) implies that

A5

The proof can be completed by the solution of (A5).  □

Proof of Theorem 1

Define V_n(U) as

Note that

where σ_i(·) denotes the ith-largest eigenvalue of a matrix. Since

we conclude that

On the other hand,

Together with the fact that

nV_n(U) can be re-written as

where ⁠. Therefore, ⁠, in distribution. Since both V(U) and nV_n(U) are convex and V(U) has a unique minimum, it follows that,

  □

Proof of Theorem 2

The proof proceeds in the same fashion as that of Theorem 1. Define V_n(U) as

As before,

Note that if c_ij = 0, in probability and ⁠. Therefore, the above expression can be rewritten as

Since ⁠, we conclude that the minimizer of nV_n(U) satisfies u_ij = 0 if c_ij = 0 with probability tending to one. The proof is now completed if we note that the maximum likelihood estimator Ĉ^true for the true graph (V, E = (c_ij ≠ 0)) is such that

in distribution, where the minimum is taken over all symmetric matrices U such that u_ij = 0 if c_ij = 0.  □

Proof of Lemma 2

Simple matrix calculus shows that the matrix of second derivatives of the objective function in (8) is positive definite and therefore the objective function is strictly convex. Since the feasible region is compact, Ĉ_new is always well defined. We now show that the algorithm will terminate in a finite number of iterations. Note that, at each iteration, Ĉ_old lies in the feasible region of Step 2. If the algorithm does not terminate, that is, at each step Ĉ_new ≠ Ĉ_old, then the minimum attained at Step 2 is strictly smaller than that from the previous iteration. The minima attained in the iterations form a strictly decreasing sequence, which in turn implies that the sign matrix in (8) must be different for all iterations. However, this contradicts the fact that there are only a finite number, 2^p(p−1)/2, of possible choices for the sign matrix S. Therefore the algorithm has to terminate.

Now we show that the algorithm converges to the solution to (3). Denote the solution at convergence of the algorithm by Ĉ. By the algorithm we see there exist two sign matrices Ŝ and ⁠, with ⁠, ⁠, and for any ĉ_ij = 0, such that Ĉ solves (8) with both Ŝ and ⁠. Let ⁠. Then, by the Karush-Kuhn-Tucker conditions (Boyd & Vandenberghe, 2003), there exist λ₁ > 0 and λ₂ > 0 such that

A6

A7

and

A8

A9

Together with the fact that for any ĉ_ij ≠ 0, (A6) and (A8) imply that λ₁ = λ₂ ≡ λ. Combining this with (A7) and (A9), we conclude that

which implies that Ĉ is also the solution to (4), and equivalently (3), again by the Karush-Kuhn-Tucker conditions.  □

Proof of Lemma 3

Let B = A⁻¹. Then b_ii = 1 according to our scaling and b_i,−i is the least-squares estimator corresponding to regressing X⁽ⁱ⁾ on the other elements; see Lauritzen (1996) and Meinshausen & Bühlmann (2006). Using this fact, we may write (11) as

To minimize this function, we have θ_ii = 1 and θ_i,−i as the minimizer of (10).  □

References