Table Of ContentSparse Nonparametric Graphical Models
John Lafferty, Han Liu and Larry Wasserman
January 1, 2012
Abstract: We present some nonparametric methods for graphical modeling. In the
discrete case, where the data are binary or drawn from a finite alphabet, Markov
random fields are already essentially nonparametric, since the cliques can take only a
2 finite number of values. Continuous data are different. The Gaussian graphical model
1 is the standard parametric model for continuous data, but it makes distributional
0
assumptions that are often unrealistic. We discuss two approaches to building more
2
flexible graphical models. One allows arbitrary graphs and a nonparametric extension
n
of the Gaussian; the other uses kernel density estimation and restricts the graphs to
a
J trees and forests. Examples of both methods are presented. We also discuss possible
4 future research directions for nonparametric graphical modeling.
]
L
Contents
M
1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
.
t
a 2 Two Families of Nonparametric Graphical Models . . . . . . . . . . . . . . . . . . 2
t
s 3 The Nonparanormal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
[
3.1 Connection to Copula . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1
3.2 Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
v
4 3.3 Statistical Properties of S (f(cid:101)) . . . . . . . . . . . . . . . . . . . . . . . . . . 9
9 n
7 4 Forest Density Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
0
4.1 A Two-Step Procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
.
1
4.2 Statistical Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
0
2 5 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1
: 5.1 Gene-Gene Interaction Graphs . . . . . . . . . . . . . . . . . . . . . . . . . . 17
v
i 5.2 Graphs for Equities Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
X
6 Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
r
a 7 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
1. Introduction
This paper presents two methods for constructing nonparametric graphical models for con-
tinuous data. In the discrete case, where the data are binary or drawn from a finite alphabet,
Markov random fields or log-linear models are already essentially nonparametric, since the
cliques can take only a finite number of values. Continuous data are different. The Gaus-
sian graphical model is the standard parametric model for continuous data, but it makes
1
distributional assumptions that are typically unrealistic. Yet few practical alternatives to
the Gaussian graphical model exist, particularly for high dimensional data. We discuss two
approaches to building more flexible graphical models that exploit sparsity. These two ap-
proaches are at different extremes in the array of choices available. One allows arbitrary
graphs, but makes a distributional restriction through the use of copulas; this is a semi-
parametric extension of the Gaussian. The other approach uses kernel density estimation
and restricts the graphs to trees and forests; in this case the model is fully nonparametric,
at the expense of structural restrictions. We describe two-step estimation methods for both
approaches. We also outline some statistical theory for the methods, and compare them in
some examples. This article is in part a digest of two recent research articles where these
methods first appeared, Liu et al. (2009) and Liu et al. (2011).
The methods we present here are relatively simple, and many more possibilities remain for
nonparametric graphical modeling. But as we hope to demonstrate, a little nonparametricity
can go a long way.
2. Two Families of Nonparametric Graphical Models
The graph of a random vector is a useful way of exploring the underlying distribution.
If X = (X ,...,X ) is a random vector with distribution P, then the undirected graph
1 d
G = (V,E) corresponding to P consists of a vertex set V and an edge set E where V has d
elements, one for each variable X . The edge between (i,j) is excluded from E if and only
i
if X is independent of X given the other variables X ≡ (X : 1 ≤ s ≤ d, s (cid:54)= i,j),
i j \{i,j} s
written
(cid:12)
X (cid:113)X (cid:12) X . (2.1)
i j (cid:12) \{i,j}
The general form for a (strictly positive) probability density encoded by an undirected
graph G is
1 (cid:88)
p(x) = exp fC(xC), (2.2)
Z(f)
C∈Cliques(G)
wherethesumisoverallcliques,orfullyconnectedsubsetsofverticesofthegraph.Ingeneral,
this is what we mean by a nonparametric graphical model. It is the graphical model analogue
of the general nonparametric regression model. Model (2.2) has two main ingredients, the
graph G and the functions {f }. However, without further assumptions, it is much too
C
general to be practical. The main difficulty in working with such a model is the normalizing
constant Z(f), which cannot, in general, be efficiently computed or approximated.
In the spirit of nonparametric estimation, we can seek to impose structure on either
the graph or the functions f in order to get a flexible and useful family of models. One
C
approach parallels the ideas behind sparse additive models for regression. Specifically, we
replace the random variable X = (X ,...,X ) by the transformed random variable f(X) =
1 d
2
nonparanormal forest densities
univariate marginals nonparametric nonparametric
bivariate marginals determined by Gaussian copula nonparametric
graph unrestricted acyclic
Fig 1. Comparison of properties of the nonparanormal and forest-structured densities.
(f (X ),...,f (X )), and assume that f(X) is multivariate Gaussian. This results in a non-
1 1 d d
parametric extension of the Normal that we call the nonparanormal distribution. The non-
paranormal depends on the univariate functions {f }, and a mean µ and covariance matrix
j
Σ, all of which are to be estimated from data. While the resulting family of distributions
is much richer than the standard parametric Normal (the paranormal), the independence
relations among the variables are still encoded in the precision matrix Ω = Σ−1, as we show
below.
The second approach is to force the graphical structure to be a tree or forest, where each
pair of vertices is connected by at most one path. Thus, we relax the distributional assump-
tion of normality but we restrict the allowed family of undirected graphs. The complexity of
the model is then regulated by selecting the edges to include, using cross validation.
Figure 1 summarizes the tradeoffs made by these two families of models. The nonpara-
normal can be thought of as an extension of additive models for regression to graphical
modeling. This requires estimating the univariate marginals; in the copula approach, this is
done by estimating the functions f (x) = µ + σ Φ−1(F (x)), where F is the distribution
j j j j j
function for variable X . After estimating each f , we transform to (assumed) jointly Normal
j j
via Z = (f (X ),...,f (X )) and then apply methods for Gaussian graphical models to esti-
1 1 d d
mate the graph. In this approach, the univariate marginals are fully nonparametric, and the
sparsity of the model is regulated through the inverse covariance matrix, as for the graphical
lasso, or “glasso” (Banerjee et al., 2008; Friedman et al., 2007)1 The model is estimated in a
two-stage procedure; first the functions f are estimated, and then inverse covariance matrix
j
Ωisestimated.Thehighlevelrelationshipbetweenlinearregressionmodels,Gaussiangraph-
ical models, and their extensions to additive and high dimensional models is summarized in
Figure 2.
In the forest graph approach, we restrict the graph to be acyclic, and estimate the bi-
variate marginals p(x ,x ) nonparametrically. In light of equation (4.1), this yields the full
i j
1Throughoutthepaperweusethetermgraphicallasso,orglasso,coinedbyFriedmanetal.(2007)torefer
tothesolutionobtainedby(cid:96) -regularizedlog-likelihoodundertheGaussiangraphicalmodel.Thisestimator
1
goes back at least to Yuan and Lin (2007), and an iterative lasso algorithm for doing the optimization was
first proposed by Banerjee et al. (2008). In our experiments we use the R packages glasso (Friedman et al.,
2007) and huge to implement this algorithm.
3
assumptions dimension regression graphical Models
low linear model multivariate Normal
parametric
high lasso graphical lasso
low additive model nonparanormal
nonparametric
high sparse additive model sparse nonparanormal
Fig 2. Comparison of regression and graphical models. The nonparanormal extends additive models to the
graphical model setting. Regularizing the inverse covariance leads to an extension to high dimensions, which
parallels sparse additive models for regression.
nonparametric family of graphical models having acyclic graphs. Here again, the estimation
procedure is two-stage; first the marginals are estimated, and then the graph is estimated.
Sparsity is regulated through the edges (i,j) that are included in the forest.
Clearly these are just two tractable families within the very large space of possible non-
parametric graphical models specified by equation (2.2). Many interesting research possibil-
ities remain for novel nonparametric graphical models that make different assumptions; we
discuss some possibilities in a concluding section. We now discuss details of these two model
families, beginning with the nonparanormal.
3. The Nonparanormal
We say that a random vector X = (X ,...,X )T has a nonparanormal distribution and write
1 d
X ∼ NPN(µ,Σ,f)
in case there exist functions {f }d such that Z ≡ f(X) ∼ N(µ,Σ), where f(X) =
j j=1
(f (X ),...,f (X )). When the f ’s are monotone and differentiable, the joint probability
1 1 d d j
density function of X is given by
(cid:26) (cid:27) d
1 1 (cid:89)
p (x) = exp − (f(x)−µ)T Σ−1(f(x)−µ) |f(cid:48)(x )|, (3.1)
X (2π)d/2|Σ|1/2 2 j j
j=1
where the product term is a Jacobian.
Note that the density in (3.1) is not identifiable—we could scale each function by a
constant, and scale the diagonal of Σ in the same way, and not change the density. To make
the family identifiable we demand that f preserves marginal means and variances:
j
µ = E(Z ) = E(X ) and σ2 ≡ Σ = Var(Z ) = Var(X ). (3.2)
j j j j jj j j
These conditions only depend on diag(Σ) but not the full covariance matrix.
4
Fig 3.Densitiesofthree2-dimensionalnonparanormals.Theleftplotshavecomponentfunctionsoftheform
f (x) = sign(x)|x|α, with α = 0.9, and α = 0.8. The center plots have component functions of the form
α 1 2
g (x)=(cid:98)x(cid:99)+1/(1+exp(−α(x−(cid:98)x(cid:99)−1/2))) with α =10 and α =5, where x−(cid:98)x(cid:99) is the fractional part.
α 1 2
The right plots have component functions of the form h (x)=x+sin(αx)/α, with α =5 and α =10. In
α 1 2
each case µ=(0,0) and Σ=(cid:0)1 .5(cid:1).
.5 1
Now, let F (x) denote the marginal distribution function of X . Since the component
j j
f (X ) is Gaussian, we have that
j j
(cid:18) (cid:19)
f (x)−µ
F (x) = P(X ≤ x) = P(Z ≤ f (x)) = Φ j j
j j j j
σ
j
which implies that
f (x) = µ +σ Φ−1(F (x)). (3.3)
j j j j
The form of the density in (3.1) implies that the conditional independence graph of the
nonparanormal is encoded in Ω = Σ−1, as for the parametric Normal, since the density
factors with respect to the graph of Ω, and therefore obeys the global Markov property of
the graph.
In fact, this is true for any choice of identification restrictions; thus, it is not necessary to
estimate µ or σ to estimate the graph, as the following result shows.
Lemma 3.1. Define
h (x) = Φ−1(F (x)) (3.4)
j j
and let Λ be the covariance matrix of h(X). Then X (cid:113)X |X if and only if Λ−1 = 0.
j k \{j,k} jk
5
Proof. We can rewrite the covariance matrix as
Σ = Cov(Z ,Z ) = σ σ Cov(h (X ),h (X )).
jk j k j k j j k k
Hence Σ = DΛD and
Σ−1 = D−1Λ−1D−1,
where D is the diagonal matrix with diag(D) = σ. The zero pattern of Λ−1 is therefore
identical to the zero pattern of Σ−1.
Figure 3 shows three examples of 2-dimensional nonparanormal densities. The component
functions are taken to be from three different families of monotonic functions—one using
power transforms, one using logistic transforms, and another using sinusoids:
f (x) = sign(x)|x|α
α
1
g (x) = (cid:98)x(cid:99)+
α 1+exp(cid:8)−α(x−(cid:98)x(cid:99)− 1)(cid:9)
2
sin(αx)
h (x) = x+ .
α
α
The covariance in each case is Σ = (cid:0)1 .5(cid:1) and the mean is µ = (0,0). It can be seen how
.5 1
the concavity and number of modes of the density can change with different nonlinearities.
Clearly the nonparanormal family is much richer than the Normal family.
The assumption that f(X) = (f (X ),...,f (X ) is Normal leads to a semiparametric
1 1 d d
model where only one dimensional functions need to be estimated. But the monotonicity of
thefunctionsf ,whichmapontoR,enablescomputationaltractabilityofthenonparanormal.
j
For more general functions f, the normalizing constant for the density
(cid:26) (cid:27)
1
p (x) ∝ exp − (f(x)−µ)T Σ−1(f(x)−µ)
X
2
cannot be computed in closed form.
3.1. Connection to Copula
If F is the distribution of X , then U = F (X ) is uniformly distributed on (0,1). Let C
j j j j j
denote the joint distribution function of U = (U ,...,U ), and let F denote the distribution
1 d
function of X. Then we have that
F(x ,...,x ) = P(X ≤ x ,...,X ≤ x ) (3.5)
1 d 1 1 d d
= P(F (X ) ≤ F (x ),...,F (X ) ≤ F (x )) (3.6)
1 1 1 1 d d d d
= P(U ≤ F (x ),...,U ≤ F (x )) (3.7)
1 1 1 d d d
= C(F (x ),...,F (x )). (3.8)
1 1 d d
6
This is known as Sklar’s theorem (Sklar, 1959), and C is called a copula. If c is the density
function of C then
d
(cid:89)
p(x ,...,x ) = c(F (x ),...,F (x )) p(x ) (3.9)
1 d 1 1 d d j
j=1
where p(x ) is the marginal density of X . For the nonparanormal we have
j j
F(x ,...,x ) = Φ (Φ−1(F (x )),...,Φ−1(F (x ))) (3.10)
1 d µ,Σ 1 1 d d
where Φ is the multivariate Gaussian cdf and Φ is the univariate standard Gaussian cdf.
µ,Σ
The Gaussian copula is usually expressed in terms of the correlation matrix, which is
given by R = diag(σ)−1Σdiag(σ)−1. Note that the univariate marginal density for a Normal
can be written as p(x ) = 1 φ(u ) where u = (x −µ )/σ . The multivariate Normal density
j σj j j j j j
can thus be expressed as
(cid:18) (cid:19)
1 1
p (x ,...,x ) = exp − uTR−1u (3.11)
µ,Σ 1 d (2π)d/2|R|1/2(cid:81)d σ 2
j=1 j
(cid:18) (cid:19) d
1 1 (cid:89) φ(u )
= exp − uT(R−1 −I)u j . (3.12)
|R|1/2 2 σ
j
j=1
Since the distribution F of the jth variable satisfies F (x ) = Φ((x −µ )/σ ) = Φ(u ), we
j j j j j j j
have that (X −µ )/σ =d Φ−1(F (X )). The Gaussian copula density is thus
j j j j j
1 (cid:110) 1 (cid:111)
c(F (x ),...,F (x )) = exp − Φ−1(F(x))T(R−1 −I)Φ−1(F(x))
1 1 d d |R|1/2 2
(3.13)
where Φ−1(F(x)) = (Φ−1(F (x )),...,Φ−1(F (x ))). This is seen to be equivalent to (3.1)
1 1 d d
using the chain rule and the identity
1
(Φ−1)(cid:48)(η) = . (3.14)
φ(Φ−1(η))
3.2. Estimation
Let X(1),...,X(n) be a sample of size n where X(i) = (X(i),...,X(i))T ∈ Rd. We’ll design
1 d
a two-step estimation procedure where first the functions f are estimated, and then the
j
inverse covariance matrix Ω is estimated, after transforming to approximately Normal.
In light of (3.4) we define
(cid:98)h (x) = Φ−1(F(cid:101) (x))
j j
7
where F(cid:101) is an estimator of F . A natural candidate for F(cid:101) is the marginal empirical distri-
j j j
bution function
n
1 (cid:88)
F(cid:98)j(t) ≡ n 1(cid:110)X(i)≤t(cid:111).
j
i=1
However, in this case (cid:98)h (x) blows up at the largest and smallest values of X(i). For the
j j
high dimensional setting where n is small relative to d, an attractive alternative is to use a
truncated or Winsorized2 estimator:
δn if F(cid:98)j(x) < δn
F(cid:101)j(x) = F(cid:98)j(x) if δn ≤ F(cid:98)j(x) ≤ 1−δn (3.15)
(1−δ ) if F(cid:98) (x) > 1−δ ,
n j n
whereδ isatruncationparameter.Thereisabias-variancetradeoffinchoosingδ ;increasing
n n
δ increases the bias while it decreases the variance.
n
GiventhisestimateofthedistributionofvariableX ,wethenestimatethetransformation
j
function f by
j
f(cid:101)(x) ≡ µ +σ (cid:101)h (x) (3.16)
j (cid:98)j (cid:98)j j
where
(cid:16) (cid:17)
(cid:101)h (x) = Φ−1 F(cid:101) (x)
j j
and µ and σ are the sample mean and standard deviation:
(cid:98)j (cid:98)j
(cid:118)
n (cid:117) n
µ ≡ 1 (cid:88)X(i) and σ = (cid:117)(cid:116)1 (cid:88)(cid:16)X(i) −µ (cid:17)2.
(cid:98)j n j (cid:98)j n j (cid:98)j
i=1 i=1
Now, let S (f(cid:101)) be the sample covariance matrix of f(cid:101)(X(1)),...,f(cid:101)(X(n)); that is,
n
n
1 (cid:88)(cid:16) (cid:17)(cid:16) (cid:17)T
S (f(cid:101)) ≡ f(cid:101)(X(i))−µ (f(cid:101)) f(cid:101)(X(i))−µ (f(cid:101)) (3.17)
n n n
n
i=1
n
1 (cid:88)
µ (f(cid:101)) ≡ f(cid:101)(X(i)).
n
n
i=1
We then estimate Ω using S (f(cid:101)). For instance, the maximum likelihood estimator is Ω(cid:98)MLE =
n n
S (f(cid:101))−1.
n
The (cid:96) -regularized estimator is
1
(cid:110) (cid:16) (cid:17) (cid:111)
Ω(cid:98) = argmin tr ΩS (f(cid:101)) −log|Ω|+λ(cid:107)Ω(cid:107) (3.18)
n n 1
Ω
2After Charles P. Winsor, the statistician whom John Tukey credited with his conversion from topology
to statistics (Mallows, 1990).
8
where λ is a regularization parameter, and (cid:107)Ω(cid:107) = (cid:80)d (cid:80)d |Ω |. The estimated graph
1 j=1 k=1 jk
is then E(cid:98) = {(j,k) : Ω(cid:98) (cid:54)= 0}.
n jk
Thus, we use a two-step procedure to estimate the graph.
1. Replace the observations, for each variable, by their respective Normal scores, subject
to a Winsorized truncation.
2. Apply the graphical lasso to the transformed data to estimate the undirected graph.
The first step is non-iterative and computationally efficient. The truncation parameter δ
n
is chosen to be
1
δ = √ (3.19)
n
4n1/4 πlogn
and does not need to be tuned. As will be shown in Theorem 3.1, such a choice makes the
nonparanormal amenable to theoretical analysis.
3.3. Statistical Properties of S (f(cid:101))
n
The main technical result is an analysis of the covariance of the Winsorized estimator above.
In particular, we show that under appropriate conditions,
(cid:115)
(cid:12) (cid:12) logd+log2n
mj,akx(cid:12)(cid:12)Sn(f(cid:101))jk −Sn(f)jk(cid:12)(cid:12) = OP n1/2
where S (f(cid:101)) denotes the (j,k) entry of the matrix S (f(cid:101)). This result allows us to leverage
n jk n
the significant body of theory on the graphical lasso (Ravikumar et al., 2009; Rothman et al.,
2008) which we apply in step two.
Theorem 3.1. Suppose that d = nξ and let f(cid:101)be the Winsorized estimator defined in (3.16)
1
with δ = √ . Define
n
4n1/4 πlogn
48 (cid:16)√ (cid:17)
C(M,ξ) ≡ √ 2M −1 (M +2)
πξ
(cid:115)
logd+log2n
for M,ξ > 0. Then for any (cid:15) ≥ C(M,ξ) and sufficiently large n, we have
n1/2
(cid:18) (cid:12) (cid:12) (cid:19) c d c d (cid:18) c n1/2(cid:15)2 (cid:19)
P max(cid:12)S (f(cid:101)) −S (f) (cid:12) > (cid:15) ≤ 1 + 2 +c exp − 4 ,
jk (cid:12) n jk n jk(cid:12) (n(cid:15)2)2ξ nMξ−1 3 logd+log2n
where c ,c ,c ,c are positive constants.
1 2 3 4
9
The proof of this result involves a detailed Gaussian tail analysis, and is given in Liu et al.
(2009).
Using Theorem 3.1 and the results of Rothman et al. (2008) it can then be shown that the
precisionmatrixisestimatedatthefollowingratesintheFrobeniusnormandthe(cid:96) -operator
2
norm.
(cid:115)
(s+d)logd+log2n
(cid:107)Ω(cid:98)n −Ω0(cid:107)F = OP n1/2
and
(cid:115)
slogd+log2n
(cid:107)Ω(cid:98)n −Ω0(cid:107)2 = OP n1/2 ,
where
s ≡ Card({(i,j) ∈ {1,...,d}×{1,...,d}|Ω (i,j) (cid:54)= 0, i (cid:54)= j})
0
is the number of nonzero off-diagonal elements of the true precision matrix.
Using the results of Ravikumar et al. (2009), it can also be shown, under appropriate
conditions, that the sparsity pattern of the precision matrix is estimated accurately with
high probability. In particular, the nonparanormal estimator Ω(cid:98) satisfies
n
(cid:16) (cid:16) (cid:17)(cid:17)
P G Ω(cid:98) ,Ω ≥ 1−o(1)
n 0
where G(Ω(cid:98) ,Ω ) is the event
n 0
(cid:110) (cid:16) (cid:17) (cid:111)
(cid:0) (cid:1)
sign Ω(cid:98) (j,k) = sign Ω−1(j,k) , ∀j,k ∈ {1,...,d} .
n 0
We refer to Liu et al. (2009) for the details of the conditions and proofs. These O(cid:101) (n−1/4)
P
ratesareslowerthantheO(cid:101) (n−1/2)ratesobtainableforthegraphicallasso.However,inmore
P
recent work (Liu et al., 2012) we use estimators based on Spearman’s rho and Kendall’s tau
statistics to obtain the parametric rate.
4. Forest Density Estimation
We now describe a very different, but equally flexible and useful approach. Rather than
assuming a transformation to normality and an arbitrary undirected graph, we restrict the
graph to be a tree or forest, but allow arbitrary nonparametric distributions.
Let p∗(x) be a probability density with respect to Lebesgue measure µ(·) on Rd and let
X(1),...,X(n) be n independent identically distributed Rd-valued data vectors sampled from
p∗(x)whereX(i) = (X(i),...,X(i)).LetX denotetherangeofX(i) andletX = X ×···×X .
1 d j j 1 d
10
