Graphical models and networks models are important topics in the modern statistical and machine learning research. Both methods use the graphs a lot but they are quite different ideas. In graphical models, graphs are not our data but are tools to determine relationship among entries of a random vector. In network models, the graph are the data (i.e., we observe the network) and we want to make inference with this type of data (known as network data).
A graphical model uses a graph to represent the conditional independence between a set of RVs. We start with the concepts of graphical models and later we will discuss how this model is constructed. Suppose that \(X{\perp\mspace{-10.5mu}\perp}Y| Z\) then we have \[p_{XYZ}(x,y,z) = p(x,y|z)p(z) = p(x|z) p(y|z) p(z)= g(x,z) h(y,z)\] for some functions \(g\) and \(h\). We then use the following graph to represent their relation:
The edge \(X-Z\) is drawn because the density factorization has a factor, namely \(g(x,z)\), that depends on both \(x\) and \(z\). Similarly, the edge \(Z-Y\) is drawn because of factor \(h(y,z)\).
Note that there is no edge between \(X-Y\). The only path from \(X\) to \(Y\) passes through \(Z\). Later we will see that in the graphical model, this implies conditional independence of \(X\) and \(Y\) given \(Z\).
The above is the basic definition of a graphical model. We now discuss how this model is constructed. The graphical model relies on two properties: graph factorization (how the distribution of a random variable is associated with a graph) and Markov properties (how the graph represents conditional independence).
A graph \(G\) and a random vector \(X\) may or may not have any relationship. The notion of graph factorization connects the joint PDF/PMF of \(X\) using a graph \(G\).
Formally, a graph \(G=(V,E)\) is a pair consisting of a (finite) vertex set \(V\) and an edge set \(E\subset V\times V\). Here, we consider undirected graphs where an edge \(v-w\) is represented by the fact that \((v,w)\) and \((w,v)\) are both in \(E\). We assume no self-loops, so \((v,v)\notin E\) for all \(v\in V\).
Example 1: If \(V=\{1,2,3,4\}\) and \[E=\{(1,2),(2,1), (2,3),(3,2), (2,4),(4,2), (3,4),(4,3)\}\] then the picture is
A non-empty subset of nodes \(A\subseteq V\) is complete if there is an edge \(v-w\) between any pair of nodes \(v,w\in A\). Complete sets are also called cliques. Sometimes, clique refers to an inclusion-maximal complete set. In this case, we often call it a maximal clique. We denote the family of all complete sets/maximal cliques as \(\mathcal{C}(G)\).
In the above example, complete sets/cliques are \[\begin{gather*} \{1\},\;\{2\},\;\{3\},\;\{4\}, \{1,2\},\;\{2,3\},\;\{2,4\},\;\{3,4\}, \{2,3,4\}. \end{gather*}\] And maximal cliques are \(\{1,2\}\), \(\{2,3,4\}\).
Definition 1. Let \(X = (X_1,\cdots, X_d)\) be a random vector and \(G= (V,E)\) be a graph where \(V = \{V_1,\cdots,V_d\}\) is the node set. We say that \(X\) factorizes over/with respect to a graph \(G\) if there exist (potential) functions \(\{\psi_C\geq 0: C\in \mathcal{C}(G)\}\) such that \[p(x_1,\cdots, x_d) = \frac{1}{Z}\prod_{C\in \mathcal{C}(G)} \psi_C(x_C)\] and \(Z = \int \prod_{C\in \mathcal{C}(G)} \psi_C(x_C) dx_1,\cdots dx_d\) is known as the partition function.
Note that we call the distribution of \(X\) a Gibbs distribution with respect to \(G\) if \[p(x_1,\cdots, x_d) = \frac{1}{Z}\prod_{C\in \mathcal{C}(G)} \psi_C(x_C) = \frac{1}{Z}\exp\left(\sum_{C\in \mathcal{C}(G)}\log \psi_C(x_C)\right)\] for some positive functions \(\{\psi_C> 0: C\in \mathcal{C}(G)\}\).
Example 2: If the following graph is a graphical model of random variables \(X=(X_1,X_2,X_3,X_4)\):
then \[% \begin{multline*} p_X(x_1,x_2,x_3,x_4) = \psi_{12}(x_1,x_2) \psi_{23}(x_2,x_3)\psi_{34}(x_3,x_4)\psi_{14}(x_1,x_4). % \end{multline*}\]
Definition 1 defines the meaning of graph factorization that connects the distribution of a random vector \(X\) to a graph \(G\). However, it does not imply anything about the conditional independence. The graph factorization and conditional independence are associated via the Markov properties of graphs.
A Markov random field is a random variable \(X\) satisfying Markov properties with respect to a graph. However, this definition is not very precise because we have not yet specified what we mean by Markov properties.
It turns out that there are three common Markov properties that associate the graph factorization to the notion of conditional independence. We start with the most common type of Markov properties–global Markov property.
The global Markov property relies on the notion of path and separation of a graph. A path in \(G\) is a sequence of distinct nodes \(v_0,v_1,\dots,v_d\) s.t. there is an edge between any two consecutive nodes, \(v_{i-1}-v_i\) for \(i=1,\dots,n\). Let \(A,B,C\subset V\) be subsets of nodes. Then \(C\) separates \(A\) and \(B\) if every path from a node \(v\in A\) to a node \(w\in B\) intersects \(C\). For instance, in example 1, \(X_2\) separates \(X_1\) and \((X_3,X_4)\) and in example 2, \((X_2,X_4)\) separates \(X_1\) and \(X_3\).
Definition 2 (Global Markov Property). A probability distribution \(P\) for a random vector \(X = (X_1,\cdots,X_d)\) satisfies the global Markov property with respect to a graph \(G\) if for any disjoint vertex subsets \(A, B\), and \(C\) such that \(C\) separates \(A\) and \(B\), then the random variables \(X_A\) are conditionally independent of \(X_B\) given \(X_C\).
It is very easy to see that
Graph Factorization (definition \(\ref{def::factor}\)) \(\Rightarrow\) Global Markov Property
as stated in the following theorem.
Theorem 3 (Global Markov theory). Suppose the distribution of \(X=(X_v:v\in V)\) factorizes over \(G=(V,E)\). Let \(A,B,C\subset V\) be subsets of nodes. Then \[C \ \; \text{separates}\ \; A \ \text{and} \ B \;\implies\; X_A{\perp\mspace{-10.5mu}\perp} X_B\mid X_C.\]
A distribution that satisfies the global Markov property is said to be a Markov random field or Markov network with respect to the graph. A more general type of Markov property is the local Markov property, which is defined as follows.
Definition 4 (Local Markov Property). A probability distribution \(P\) for a random vector \(X = (X_1,\cdots,X_d)\) satisfies the local Markov property with respect to a graph \(G\) if the conditional distribution of a variable given all its neighbor is independent of any other vertices. Namely, let \(N(j) = \{X_i: E_{ij} =1\}\) be the neighbors of \(X_j\). Then the local Markov property means that \[P(X_j|X_{- j}) = P(X_j|X_{N(j)}),\] where \(X_{-j} = \{X_i: i\neq j\}\).
A more general definition is the pairwise Markov property.
Definition 5 (Pairwise Markov Property). A probability distribution \(P\) for a random vector \(X = (X_1,\cdots,X_d)\) satisfies the pairwise Markov property with respect to a graph \(G\) if for any two non-adjacent vertices \(X_i\) and \(X_j\) (i.e., \(E_{ij}=0\)), \[X_i{\perp\mspace{-10.5mu}\perp}X_j|X_{V\backslash \{i,j\}}.\]
Proposition 6 (Equivalence of Markov properties). For any undirected graph \(G\) and any distribution \(P\), we have \[\mbox{Global Markov Property}\Rightarrow \mbox{Local Markov Property}\Rightarrow \mbox{Pairwise Markov Property}.\]
The proof is very straightforward so we omit it.
Example: local Markov property but no global Markov property. Define binary random variables \(X_1,\cdots, X_5\) such that \(P(X_1 =1)=P(X_5=1) = \frac{1}{2}\) and \(X_2=X_1\) and \(X_4=X_5\) and \(X_3 = X_2X_4\). You can easily verify that the random vector satisfies the local Markov property with respect to the chain graph \(G_0\) that is \((X_1-X_2-X_3-X_4-X_5).\) In particular, the PMF of \(X_3\) is conditionally independent of \(X_1\) and \(X_5\) given \(X_2\) and \(X_4\). However, the global Markov property is violated. To see this, consider the case of \(X_3 = 0\) and it is easy to see that \(P(x_2,x_4|X_3=0) = \frac{1}{3}\) when \((x_2,x_4) = (1,0), (0,1), (0,0)\). However, the marginal probability \(P(X_2=0|X_3=0)=P(X_4=0|X_3=0) = \frac{2}{3}\). Thus, \[P(X_2=0,X_4=0|X_3=0) = \frac{1}{3}\neq P(X_2=0|X_3=0)\times P(X_4=0|X_3=0) = \frac{4}{9}\] so the global Markov property does not hold for this graph \(G_0\).
Example: pairwise Markov property but no local Markov property. Define binary random variables \(X_1,X_2,X_3\) and \(X_1=X_2=X_3\) with \(P(X_1=1) = \frac{1}{2}\). The random vector \(X = (X_1,X_2,X_3)\) has a very degenerated PMF. Consider a graph \(G\) such that there is only one edge \(E_{23} =1\). Then you can easily verify that \(X\) satisfies the pairwise Markov property with respect to \(G\) but not the local Markov property (specifically, \(P(X_1=1|X_2=0,X_3=0) = 0\neq P(X_1=1)=\frac{1}{2}\)). This example also shows a fact about the Markov properties– the same distribution may satisfy different Markov properties on different graphs! In the above example, the same pairwise Markov property holds for another graph \(G'\) with only a single edge \(E'_{12}=1\) or a graph \(G''\) with only a single edge \(E''_{13}=1\).
From the above examples, we see that the term Markov random field is not a well-defined term by itself. We have to specifically state which Markov property a random variable \(X\) satisfies with respect to a graph.
The good news is that when the PDF/PMF is positive, the three Markov properties are equivalent.
Proposition 7. Let \(p\) be the PDF/PMF of random variable \(X\in\mathbb{R}^d\) such that \(\mathcal{X}_j\) is the support of \(X_j\) for \(j=1,\cdots, d\). If \(p(x)>0\) for all \(x\in \mathcal{X}_1\times\mathcal{X}_2\times \cdots\times\mathcal{X}_d,\) then the three Markov properties are equivalent.
The condition \[\begin{equation} p(x)>0\mbox{ for all }x\in \mathcal{X}_1\times\mathcal{X}_2\times \cdots\times\mathcal{X}_d \label{eq::positivity} \end{equation}\] is also known as the positivity condition (for graphical models).
Proposition 7 is a powerful result because it shows that under the positivity condition, all three Markov properties are equivalent. So the term Markov random fields is well-defined.
The above proposition relies on the intersection lemma from
Pearl, J., & Paz, A. (1985). Graphoids: A graph-based logic for reasoning about relevance relations. University of California (Los Angeles). Computer Science Department.
Lemma 8 (Intersection lemma; Pearl, J., & Paz (1985)). Suppose that for any subsets \(A,B,C,D\subset V\) we have \[X_A{\perp\mspace{-10.5mu}\perp}X_B| X_{C\cup D},\quad X_A{\perp\mspace{-10.5mu}\perp}X_C| X_{B\cup D}\Rightarrow X_A{\perp\mspace{-10.5mu}\perp}X_{B\cup C}| X_{D}.\] Then the three Markov properties are equivalent.
[Proof of Proposition 7]
Without loss of generality, we consider the three-variable case: \(X= (X_1,X_2,X_3)\). To use Lemma 8, we need to show that \[X_1{\perp\mspace{-10.5mu}\perp}X_2|X_3,\quad X_1{\perp\mspace{-10.5mu}\perp}X_3|X_2\Rightarrow X_1 {\perp\mspace{-10.5mu}\perp}\{X_2,X_3\}.\]
Assume the two conditional independences on the left-hand side of the above equation. Then we have \[p(x_1,x_2,x_3) = f_{13}(x_1,x_3) f_{23}(x_2,x_3) = g_{12}(x_1,x_2) g_{23}(x_2,x_3)\] for some functions \(f_{13}, f_{23}, g_{12}, g_{23}\). Thus, \[g_{12}(x_1,x_2) = \frac{f_{13}(x_1,x_3) f_{23}(x_2,x_3)}{g_{23}(x_2,x_3)} = f_{13}(x_1,x_3)\frac{ f_{23}(x_2,x_3)}{g_{23}(x_2,x_3)}.\] An interesting implication from the above equation is that the left-hand side does not depend on \(x_3\) so this holds for any \(x_3\). WLOG, we choose \(x_3=0\) and this leads to \[g_{12}(x_1,x_2) = f_{13}(x_1,0)\frac{ f_{23}(x_2,0)}{g_{23}(x_2,0)} = h(x_1) k(x_2).\] Putting this back to the joint PDF/PMF, we obtain \[p(x_1,x_2,x_3) = g_{12}(x_1,x_2) g_{23}(x_2,x_3) = h(x_1) k(x_2) g_{23}(x_2,x_3),\] which implies \(X_1 {\perp\mspace{-10.5mu}\perp}\{X_2,X_3\}\). So by Lemma 8, the three Markov properties are equivalent.
Theorem 3 shows that if the distribution of a random vector \(X\) factorizes over a graph, then it satisfies the global Markov property. However, the reverse direction is unclear to us. Specifically, we want to know
if a random vector satisfies the global/local/pairwise Markov property with respect to a graph, can it always be factorized with respect to the graph?
If we can provide a positive answer to the above question, then we can use graph factorization for a Markov random field. The good news is: the following theorem, known as the Hammersley-Clifford (or Hammersley-Clifford-Besag) theorem, provides a positive answer to this question.
Theorem 9 (Hammersley-Clifford (1971)). Suppose that \(G = (V, E)\) is a graph and \(X_1,\cdots, X_d\) are random variables that take on a finite number of values. Let \(\mathcal{X}_j\) be the support of \(X_j\) for \(j=1,\cdots, d\). If \(p(x)>0\) for all \(x\in \mathcal{X}_1\times\mathcal{X}_2\times \cdots\times\mathcal{X}_d,\) and satisfies the local Markov property with respect to \(G\), then it factors with respect to \(G\).
The Hammersley-Clifford theorem shows that under the positivity condition in equation \(\eqref{eq::positivity}\), the graph factorization and all three Markov properties are equivalent.
The following paper is the original paper that states this theorem:
Hammersley, J. M., & Clifford, P. (1971). Markov fields on finite graphs and lattices.
Note that they did not publish this paper in a journal article but you can still find the original manuscript online.
A formal paper that includes this theorem (and improves the proof and mentions the generalization to continuous random vectors) is the following paper:
Besag, J. (1974). Spatial interaction and the statistical analysis of lattice systems. Journal of the Royal Statistical Society: Series B (Methodological), 36(2), 192-225.
Although the Hammersley-Clifford theorem only proves the case for discrete random variables, the result can be generalized to continuous random variables as well. The Hammersley-Clifford theorem together with Proposition 7 imply the following conclusion:
For a random vector \(X\) with a positive PDF/PMF (satisfying the positivity condition), then \[\mbox{satisfying Markov Properties}\Leftrightarrow \mbox{factorizing with respect to $G$}.\]
Thus, Theorem 3 together with the Hammersley-Clifford theorem provide the foundation of graphical models that we can interchangeably use graph factorization and conditional independence. This is why the Hammersley-Clifford theorem is sometimes referred to as the fundamental theorem of graphical models.
The Hammersley-Clifford theorem not only informs us that satisfying the three Markov properties is equivalent to graph factorization, its required condition, the positivity condition in equation \(\eqref{eq::positivity}\), is the same as the condition for showing equivalence among the three Markov properties. Thus, once a random variable \(X\) satisfies the positivity condition, all three Markov properties are equivalent (Proposition 7) and it factorizes with respect to a graph \(G\) (Theorem 9), so the Markov random field is well-defined and factorizes with respect to a graph \(G\).
Consider the problem of a Gaussian random vector \(X = (X_1,X_2,\cdots, X_p)\in\mathbb{R}^p\) with a mean vector \(\mu\) and a covariance matrix \(\Sigma\). Assume that \(\Sigma\) is positive definite, then the joint PDF can be written as \[p_X(x) = \frac{1}{\sqrt{(2\pi)^p {\sf det}(\Sigma)}}\exp\left\{-\frac{1}{2}(x-\mu)^T \Sigma^{-1}(x-\mu)\right\},\] where \(x=(x_1,\cdots, x_p)\).
In this model, there are two parameters \(\mu\) and \(\Sigma\). When \(\Sigma\) is invertible, the PDF satisfies the positivity condition, so \(X\) is a well-defined Markov random field with respect to a particular graph \(G\) and we can use graph factorization and all Markov definitions.
Now we consider the case where \(G\) has no edge between node \(X_1\) and \(X_2\). This implies the conditional independence \(X_1{\perp\mspace{-10.5mu}\perp}X_2|X_3,\cdots, X_p\). What does this conditional independence tell us about the underlying parameters?
Using the graph factorization, we can factorize \(p_X\) into \[p_X(x) = g(x_1,x_3,x_4,\cdots, x_p) h(x_2,x_3,\cdots, x_p).\] Therefore, \[\log p_X(x) = \tilde{g}(x_1,x_3,x_4,\cdots, x_p)+ \tilde{h}(x_2,x_3,\cdots, x_p) = -\frac{1}{2}(x-\mu)^T \Sigma^{-1}(x-\mu)+ C_0,\] where \(C_0\) is a constant with respect to \(x\).
Because \[(x-\mu)^T \Sigma^{-1}(x-\mu) = \sum_{i,j=1}^p (x_i-\mu_i)(x_j-\mu_j) \left(\Sigma^{-1}\right)_{ij},\] we conclude that \(\left(\Sigma^{-1}\right)_{12}=0\). Namely, for a Gaussian random vector, if we see the \((i,j)\)-th element of the inverse covariance matrix (also known as the precision matrix) is \(0\), we have the conditional independence of \(X_i\) and \(X_j\) given the other elements.
The log-linear model is a parametrization for the PMF of multinomials. Suppose that each \(X_j \in \{0,1,2,\cdots, m_j-1\}\) for each \(j=1,\cdots, d\) and \(X=(X_1,\cdots, X_{d})\) is the random vector of interest.
Log-linear model. The log-linear model expands the log PMF of \(X\) as \[\begin{equation} \log p(x) = \sum_{A\subset V}\psi_A (x_A), \label{eq::LL} \end{equation}\] with the constraint that if a variable \(j\in A\) with \(x_j=0\), \(\psi_A(x_A) = 0\). Equation \(\eqref{eq::LL}\) is known as the log-linear expansion of \(p(x)\). Although \(\psi_A (x_A)\) behaves like a function, it is a set of several parameters since the variable(s) \(x_A\) only takes discrete values. In fact, there are only \(\prod_{j\in A} (m_j-1)\) number of possible values of \(\psi_A\) so it is often referred to as the parameter of a log-linear model. You can interpret the parameter/function \(\psi_A\) as the (joint) interaction effect of variables in \(A\).
Hierarchical model. A hierarchical log-linear model is a log-linear model such that if \(\psi_A(x_A) =0\) implies \(\psi_B(x_B) = 0\) for all \(B\supset A\). Namely, a hierarchical log-linear model has a nested structure that if a parameter \(\psi_A = 0\), any parameter that is a superset of \(A\) must be \(0\). You can interpret a hierarchical log-linear model as the model that any higher-order interaction exists only if all lower-order interactions exist. However, even if all lower-order interactions exist, the higher-order interaction needs not to exist.
Graphical log-linear model. A graphical log-linear model with respect to a graph \(G\) is the log-linear model such that \(\psi_A(x_A)\) is not zero if and only if \(A\) is a clique (not necessarily a maximal clique). The graphical model requires that if all lower-order interactions exist, the higher-order interaction MUST exist. Thus, one can see that the graphical log-linear model is a sub-class of the hierarchical model, as stated in the following lemma.
Lemma 10. A graphical log-linear model is hierarchical log-linear model but not vice versa.
Suppose that for a graphical model of \(G\) with \(\psi_A = 0\), this implies that \(A\) is not a clique in \(G\). Thus, any set \(B\supset A\) will not be a clique in \(G\) so the model is hierarchical.
Now consider a three variable log-linear model with \[\log p(x) = \psi_1(x_1) +\psi_2(x_2) + \psi_3(x_3) + \psi_{12}(x_1,x_2)+\psi_{13}(x_1,x_3)+\psi_{23}(x_2,x_3).\] Clearly, this is a hierarchical model but not a graphical model (it will requires \(\psi_{123}(x_1,x_2,x_3)\neq 0\)).
With the above lemma, we conclude that
graphical log-linear model \(\Rightarrow\) hierarchical model \(\Rightarrow\) log-linear (multinomial) model
and a graphical log-linear model is both a graphical and a log-linear model. However,
log-linear model + graphical model \(\nRightarrow\) graphical log-linear model.
While the above result might sound strange, a concrete example is the famous Ising model (Remark 12).
Example 11 (Ising model). The Ising model is a special case of hierarchical log-linear models and a graphical model as well. However, it is not a graphical log-linear model. It is a hierarchical model with binary variables with only pairwise interactions. Specifically, the Ising model is the case where \[\begin{equation} \log p(x) = \sum_{i=1}^d \theta_i x_i + \sum_{(j,k)\in E} \theta_{j,k} x_jx_k. \label{eq::ising} \end{equation}\] Since the Ising model only contains pairwise interaction, it can be viewed as a discrete analogue of the Gaussian graphical model. The Ising model is related to the logistic regression. By the local Markov property, a random variable \(X_i\) only depends on its neighborhoods so the conditional probability \[P(X_i=1|X_{-i}) = P(X_i=1|X_{j}, (i,j)\in E) = \frac{\exp(\theta_i + \sum_{(i,j)\in E} \theta_{i,j}x_j)}{1+\exp(\theta_i + \sum_{(i,j)\in E} \theta_{i,j}x_j)},\] where \(X_{-i}\) is the collection of all variables except \(X_i\). While Ising model is a graphical model plus a hierarchical log-linear models, it is NOT a graphical log-linear models.
Remark 12 (Why graphical models + log-linear model is NOT graphical log-linear?). Now we consider the Ising model with three variables \(X_1,X_2,X_3\) that are fully connected. First, since Ising model is a hierarchical log-linear model, it is a log-linear model. Also, you can easily show that the joint PMF does satisfy the graph factorization, so it is a graphical model. However, since Ising model lacks the three-way interaction, it is NOT a graphical log-linear model!
Example 13 (Potts model). The Potts model is a generalized Ising model that allows variables to have \(m\) distinct outcomes, i.e., \(X_i \in\{0,1,2,\cdots, m-1\}\) and the pairwise interaction contributes only if the two variables are in the same ‘state’. Specifically, the joint PMF in the Potts model can be factorized as \[\begin{equation} \log p(x) = \sum_{i=1}^d \theta_i x_i + \sum_{(j,k)\in E} \theta_{j,k} \delta(x_j, x_k), \label{eq::potts} \end{equation}\] where \(\delta(a,b) = I(a=b).\) The Potts model is motivated by statistical mechanics in which each variable \(X_i\) is a particle and a particle has \(m\) different states. In a stable scenario, two adjacent particles (particles are variables \(X_i\)’s) will avoid being in the same state. So the distribution can be modeled using the Potts model with a negative \(\theta_{j,k}\).
A graph where the edges are directional is called a directed graph. In statistics and machine learning, we often focus on one particular directed graph called directed acyclic graphs (DAGs). A DAG is a directed graph that has no directed cycles (i.e., arrows do not form a cycle). Directed graphical models are often viewed as a generative model. To illustrate the idea, consider \(5\) random variables \(X_1,\cdots, X_5\) with the following generative models: \[\begin{align*} X_1&\sim p_1(x_1)\\ X_2&\sim p_2(x_2)\\ X_3&\sim p_3(x_3)\\ X_4|X_1,X_2&\sim p_4(x_4|X_1,X_2)\\ X_5|X_1,X_3,X_4&\sim p_5(x_5|X_1,X_3,X_4). \end{align*}\] Then we can summarize this model using the left panel of Figure 1.
Because of the popularity of DAG in the probability generative model, a DAG is also called a Bayesian network. Note that a Bayesian network has nothing to do with Bayesian inference or Bayesian statistics; it is just a graphical model that relied on Bayes’ rule to describe a probability distribution.
The DAG in the left panel of Figure 1 implies that the joint density can be written as \[p(x_1,\cdots,x_5) = p(x_5|x_1,x_3,x_4) p(x_4|x_1,x_2)p(x_3)p(x_2)p(x_1) = \psi_{1,3,4,5}(x_1,x_3,x_4,x_5) \psi_{1,2,4}(x_1,x_2,x_4)\] so the corresponding undirected graphical model is the right panel of Figure 1 that has two maximal cliques \((1,3,4,5)\) and \((1,2,4)\).
More generally, we can always convert a DAG into an UG using the idea of moralizing. If there is an arrow from node \(X_i\) to node \(X_j\), we call \(X_i\) a parent (node) of \(X_j\) and \(X_j\) a child (node) of \(X_i\). Note that every node may have multiple parents and children.
Definition 14. The moral graph \(M\) of a DAG \(G\) is an undirected graph where there is an edge between two vertices \(X_i\) and \(X_j\) if one of the following conditions met:
There is an edge between \(X_i\) and \(X_j\) in \(G\).
\(X_i\) and \(X_j\) are the parents of the same child node.
Informally, the moralized graph can be constructed by ‘marrying the parents’–we connect all parents of each child node (and remove arrows) to form the corresponding undirected graph. Two different DAGs may have the same moralized graph, as illustrated in Figure 2. Also, the arrow direction matters in the construction of moral graph; Figure 3 shows an example that we only reverse one arrow’s direction in the DAG of the left panel in Figure 2 and the resulting moral graph is different.
Important note on learning DAG. The data alone CANNOT give you a unique DAG! This is because the data can only give you associations (e.g., correlations) among variables. We can have multiple DAGs leading to the same associations. Figure 2 is an example of this – we can only identify the dependency as the undirected graph in the right but we are unable to tell which DAG (left or middle) is the actual DAG that generates the data. In actual scientific problem, we have additional knowledge on the variables, which may allow us to choose the correct DAG. But data alone CANNOT tell you which DAG is the correct one.
The d-separation is common criterion to determine independence or conditional independence in a DAG. Roughly speaking, d-separation implies independence. To explain the concept of d-separation, we first introduce a few concepts.
For a pair of vertices \(X_a\) and \(X_b\), an undirected path is a collection of vertices \((V_1,\cdots, V_k)\) such that \(V_1 = X_a\), and \(V_k = X_b\), and each pair \((V_j, V_{j+1})\) contains an arrow (ignored the orientation). A collider of a path is the triplet \((V_{j-1}, V_j,V_{j+1})\) of the path such that the arrows are \(V_{j-1}\rightarrow V_j \leftarrow V_{j+1}\). A vertex \(X_a\) is a descendent of \(X_b\) if there exists a directed path \(X_b = V_1\rightarrow V_2\rightarrow \cdots \rightarrow V_k = X_a\).
Rule 1: a collider of a path creates a block on the path. If a path has one or more blocks, then the path is d-separated.
For two sets of variables \(A,B \subset \{1,2,\cdots,n\}\), they are d-separated if all paths between \(A\) and \(B\) are d-separated.
If \(A\) and \(B\) are d-separated, then \(X_A \perp X_B\).
Rule 1 is often used to verify independence between two sets of variables. Now we consider the case of conditional independence, for which we will introduce two more rules.
Rule 2: Conditioned on a non-collider will create a block.
Rule 3: Conditioned on a collider or its descendent will remove the block caused by the collider.
For three sets of variables \(A,B,C \subset \{1,2,\cdots,n\}\), \(A\) and \(B\) are d-separated given \(C\) if all paths between \(A\) and \(B\) are d-separated after applying the above rules (Rules 2 and 3 with conditioning on \(C\)).
If \(A\) and \(B\) are d-separated given \(C\), then \(X_A\perp X_B|X_C\).
Example 1. Consider the graph that \[X_1\rightarrow X_2\rightarrow X_3\leftarrow X_4\leftarrow X_5.\] Since \(X_3\) is a collider, we have \(X_1\perp X_5\). However, if we conditioned on \(X_3\), this will remove the block at \(X_3\), making the path between \(X_1\) and \(X_5\) not d-separated. So \(X_1\) and \(X_5\) are not conditional independent given \(X_3\). If we further condition on \(X_2\), this will add an additional block at \(X_2\), making the path d-separated. So we have \(X_1\perp X_5|X_2,X_3\).
Example 2. Now we consider another graph \[Y_1\rightarrow Y_2\leftarrow Y_3\rightarrow Y_4\leftarrow Y_5.\] There are two colliders in the above DAG (\(Y_2\) and \(Y_4\)). We now consider the relation between \(Y_1\) and \(Y_5\). Clearly, \(Y_1\perp Y_5\). When we conditioned on \(Y_2\), while this removes the block at \(Y_2\), we still have another collider \(Y_4\). So we still have \(Y_1\perp Y_5 |Y_2\) and similarly \(Y_1\perp Y_5|Y_4\) Note that \(Y_1\) and \(Y_5\) will be dependent if we conditioned on \(Y_2\) and \(Y_4\). However, since \(Y_3\) is not a collider, conditioned on \(Y_3\) will create a block, making the path d-separated. So we have \(Y_1\perp Y_5|Y_3\) and even \(Y_1\perp Y_5|Y_2,Y_3,Y_4\). As long as there is one block on the path, the path will be d-separated.
A Bayesian hierarchical model is a scenario where DAGs are often applied. To illustrate the idea, we consider the following example. Suppose that we have \(n\) individuals participating in an exam and their scores can be summarized using univariate random variables \(X_1,\cdots, X_n\). We all know that scores are measurements (with noises) of the individual’s capability so we can view each random variable as \[X_i|\mu_i \sim N(\mu_i, \sigma^2),\] where \(\mu_i\) can be interpreted as the individual’s actual performance on the exam. Suppose that these \(n\) individuals are randomly chosen from a population. To model the randomness of the selection, we assume that \[\mu_1,\cdots,\mu_n\sim N(\theta,\tau^2),\] where \(\theta\) reflects the average performance of the sampled population. To account for our uncertainty about \(\theta\), we may introduce a prior \(\pi(\theta)\) over it. Under this model specification, all random quantities can be written as the DAG in Figure 4.
The DAG is also used frequently in causal inference. The arrow is interpreted as a causal relation. For instance, if we have a DAG \(X_1\rightarrow X_2\rightarrow X_3\), then we mean that \(X_1\) causes \(X_2\) and \(X_2\) causes \(X_3\). The above graph also implies that conditioning on \(X_2\), \(X_1\) and \(X_3\) are independent. In the causal relation, this means that if we control \(X_2\), then \(X_1\) does not cause any change in \(X_3\). So conditional independence becomes an elegant mathematical tool to discuss causal relations.
Here is another example to illustrate how DAGs provide useful insights on causal relation. Consider the left DAG in Figure 5–
Causal interpretation: \(X_2\) causes both \(X_1\) and \(X_3\). Thus, if \(X_2\) is unobserved, then \(X_1\) and \(X_3\) are associated (in this case, \(X_2\) is a confounder). On the other hand, if \(X_2\) is controlled, then \(X_1\) and \(X_3\) are independent.
Graphical model interpretation: The generative model is \[p(x_1,x_2,x_3) = p(x_1|x_2) p(x_3|x_2) p(x_2).\] Thus, the marginal density \[p(x_1,x_3) = \int p(x_1,x_2,x_3)dx_2 = \int p(x_1|x_2) p(x_3|x_2) p(x_2)dx_2 = g(x_1,x_3)\] for some function \(g\). Thus, \(X_1\) and \(X_3\) are marginally dependent. However, \(p(x_1,x_3|x_2) = p(x_1|x_2) p(x_3|x_2)\) so \(X_1\) and \(X_3\) are conditionally independent.
Now we consider the right DAG in Figure 5–
Causal interpretation: Both \(X_1\) and \(X_3\) causes \(X_2\) but they are independent causes. However, if \(X_2\) is observed, then \(X_1\) and \(X_3\) will be associated. Note that \(X_2\) in this case will be called a collider.
Graphical model interpretation: The generative model is \[p(x_1,x_2,x_3) = p(x_2|x_1,x_3) p(x_1) p(x_3)\Rightarrow p(x_1,x_3) = p(x_1)p(x_3)\] so \(X_1\) and \(X_3\) are marginally independent. And the conditional density \[p(x_1,x_3|x_2) = \frac{p(x_1,x_2,x_3)}{p(x_2)} = \frac{p(x_2|x_1,x_3) p(x_1) p(x_3)}{p(x_2)}\] cannot be factorized into the product of \(g_1(x_1,x_2)\) and \(g_2(x_2,x_3)\) so \(X_1\) and \(X_3\) are conditionally dependent given \(X_2\).
Therefore, the probabilistic structure implied by a DAG and the causal interpretation of variables have an elegant correspondence. This is why DAGs are very popular in causal inference.
The structural equation model (SEM) is a popular model that associates a DAG to a probability distribution via equations. For each variable \(X_i\), we let \({\sf PA}(i)\subset V\) denote the parents of \(X_i\). The structural equation model is a set of generative equations such that for each \(i\), \[\begin{equation} X_i = f_i(X_{{\sf PA}(i)}, \epsilon_i), \label{eq::NPSEM} \end{equation}\] where \(\epsilon_i\) are IID random variables with mean \(0\) and \(f_i\) are functions belonging to certain classes. Equation \(\eqref{eq::NPSEM}\) is sometimes called a nonparametric SEM, as it does not specify any parametric form.
Linear SEM. Equation \(\eqref{eq::NPSEM}\) may run into an identification problem if the class of functions is not specified properly. To resolve this issue, people often assume a specific parametric class for each equation. The most popular class is the linear function class, which leads to the famous linear SEM. The linear SEM formulates \[X_i = \gamma_{i, {\sf PA}(i)}^T X_{\sf PA(i)} +\epsilon_i.\] Very often people would further assume that \(\epsilon_i\) are IID \(N(0,\sigma^2)\)1. Under the linear SEM, you can show that the whole random vector \(X\in\mathbb{R}^n\) can be written as \[X = (\mathbf{I}_n-\Gamma)^{-1} \epsilon,\] where \(\epsilon^T = (\epsilon_1,\cdots,\epsilon_n)\) and \(\Gamma\) is constructed from \(\gamma_{i,{\sf PA}(i)}\) with \({\sf diag}(\Gamma) = (0,0,\cdots, 0)\) and the feature that if there is no arrow from \(X_j\) to \(X_i\), then \(\Gamma_{ij}=0\). Under linear SEM, if \(\epsilon\sim N(0, \sigma^2\mathbf{I}_n)\), we further have \[X\sim N(\mu, \sigma^2 (\mathbf{I}_n-\Gamma)^{-1}[(\mathbf{I}_n-\Gamma)^{-1}]^T).\] The parameters \(\sigma^2\) and \(\Gamma\) can be estimated by either the MLE or eigen-analysis on the covariance matrix of \(X\).
Example: the DAG in Figure 1. We now use the DAG in Figure 1 as an example of the linear SEM with Gaussian noises. This DAG implies the following five equations: \[\begin{align*} X_1 &= \epsilon_1\\ X_2&= \epsilon_2\\ X_3 &= \epsilon_3\\ X_4 & = \gamma_{4,1}X_1 + \gamma_{4,2}X_2+\epsilon_4\\ X_5 & = \gamma_{5,1}X_1 + \gamma_{5,3}X_3+\gamma_{5,4}X_4+\epsilon_5. \end{align*}\] This can be written as \[\begin{align*} X = \underbrace{\begin{pmatrix} 0&0&0&0&0\\ 0&0&0&0&0\\ 0&0&0&0&0\\ \gamma_{4,1}&\gamma_{4,2}&0&0&0\\ \gamma_{5,1}&0&\gamma_{5,3}&\gamma_{5,4}&0 \end{pmatrix}}_{=\Gamma} X + \epsilon, \end{align*}\] which after rearrangements, becomes \[(\mathbf{I}_5-\Gamma)X = \epsilon\Rightarrow X = (\mathbf{I}_5-\Gamma)^{-1}\epsilon.\]
Factor analysis is a special type of linear SEM that it includes latent/hidden/unobserved variables called factors. Now we consider the simplest factor model that we have \(n\) observed variables \(X = (X_1,\cdots, X_n)^T\) and \(k\) latent factors \(G = (G_1,\cdots, G_k)^T\). In this simplest case, the factor model assumes a graphical model with arrows from every \(G_j\) to every \(X_i\), which under the linear SEM, leads to the following equation \[\begin{equation} X_i = \sum_{j=1}^k \gamma_{i,j} G_j + \epsilon_i \label{eq::FM1} \end{equation}\] for each \(i\) and \(G_j,\epsilon_i\) are independent for every pair \((i,j)\). We can write equation \(\ref{eq::FM1}\) as \[X = \Gamma G + \epsilon\] with \(\mathbb{E}(\epsilon) = 0,{\sf Cov}(\epsilon) = \sigma^2 \mathbf{I}_n\). Since the latent variables \(G\) are unobserved, we often assume that \(\mathbb{E}(G) = 0, {\sf Cov}(G) = \mathbf{I}_n\) for identifiability. Under this assumption, we have \[{\sf Cov}(X) = \Gamma \Gamma^T + \sigma^2 \mathbf{I}_n.\] Thus, the problem reduces to an eigenanalysis on the covariance matrix, so the principal component analysis (PCA) is a popular method to estimate \(\Gamma\). Alternatively, we may assume that \(\epsilon,G\) are both multivariate normal, and then use the MLE (maximum likelihood estimator) to estimate \(\Gamma\).
There are many variants of the above factor models. For instances, you may use two sets of latent factors (known as the bi-factor model), or pre-specify some \(0\)’s in \(\Gamma\) so that not all factors are influencing every observables. You may even include some arrow among latent factors. The factor analysis is very popular in social sciences such as economics, psychology, and education because factors are often from prior scientific knowledge and we can use our knowledge to add/remove arrows among these variables (latent or observed).
Having explored graphs as a tool to model dependencies within a single observation (graphical models), we now shift our perspective. In the following sections, we will treat the graph itself as the observed data and discuss methods for statistical inference on networks
I would recommend the following lecture notes if you are interested in learning more about network models in Statistics:
CMU 36-720, Statistical Network Models (by C. Shalizi): https://www.stat.cmu.edu/~cshalizi/networks/16-1/
UW CSSS-STAT 567 Statistical Analysis of Social Networks (by P. Hoff): https://www.stat.washington.edu/people/pdhoff/courses/567/
In statistics, network models are often used to model network data. Unlike graphical model problems, in handling network data, we directly observe a network. Studies on statistical network models attempt to use network data to make scientific inference. There are several scenarios where network data can be used in statistical inference, for instance:
Random networks. We may view the network as random quantities (called random networks) and study the distribution that generates a random network.
Community detection. We want to find communities (nodes that are highly interconnected) within a network – these communities often represent certain groups of nodes.
Networks as covariates. In some scenarios, we may use the network as a covariate in a regression/classification task.
Sampling network data. In many realistic situations, we may not observe the complete network data but only a fraction of it. Different sampling schemes in this case lead to different estimators of the properties of the entire network.
A statistical network model is a probability model that describes the generating process of a random graph. Often the model describes the probability structure of a random undirected and unweighted graph, although many models can be generalized to directed graphs as well. In a network model, the nodes are often assumed to be fixed and non-random, and the edges are randomly formed (although this is not strict – there are network models where nodes can be randomly generated). One of the most famous network models is the Erdős-Rényi graph, which states that any pair of nodes has an equal probability of forming an edge. A generalized version of the Erdős-Rényi graph is the stochastic block model where we assume that all nodes can be partitioned into \(K\) unknown groups and pairs of nodes have different probabilities of forming an edge depending on whether they belong to the same group or not. We will briefly review and discuss some famous network models. For undirected and unweighted random networks, the probability model is equivalent to an \(n\times n\) random matrix with Bernoulli random variables in every entry. So the random network is a special case of a random matrix. Let \(G = (V,E)\) be a random graph and \(\|V\|=n\) is the number of vertices and we may use the edge/adjacency matrix \(E\in \{0,1\}^{n\times n}\) to denote the edges with \(E_{ij}=1\) means that there exists an edge between node \(i\) and node \(j\).
Erdös-Rényi model. The Erdös-Rényi model is a very simple stochastic model for generating a random graph. There are two variants of the Erdös-Rényi model.
\(ER(n,p)\) model. This variant is the model to which the Erdös-Rényi model is most commonly referred. It generates a random graph where every possible edge has an independent probability \(p\) of forming. Namely, \(P(E_{ij}=1) = p\) and \(\{E_{ij}: i\geq j\}\) are IID. Essentially, its randomness can be described by \({n\choose 2}\) independent Bernoulli random variables. This model has several interesting properties on the asymptotic behavior of \(n\) and \(p\), for instance:
If \(np<1\), then the graph will almost surely have no connected components with a size larger than \(O(\log n)\).
If \(np=1\), then the graph will almost surely have the largest connected component with a size on the order of \(O( n^{2/3})\).
If \(np\rightarrow c>1\), then the graph will almost surely have a giant connected component and no other connected component has a size larger than \(O(\log n)\).
If \(p<\frac{(1-\epsilon)\log n}{n}\) for some fixed number \(\epsilon>0\), then the graph will almost surely be disconnected.
If \(p>\frac{(1+\epsilon)\log n}{n}\) for some fixed number \(\epsilon>0\), then the graph will almost surely be (fully) connected.
The above results are summarized from the following famous paper by Erdös and Rényi:
Erdös, Paul, and Alfréd Rényi . “On the evolution of random graphs." Publ. Math. Inst. Hung. Acad. Sci 5, no. 1 (1960): 17-60.
\(ER(n,m)\) model. This variant creates a random graph with a fixed amount of edges. Here, \(n\) stands for the number of vertices and \(m\) stands for the total number of edges. The \(ER(n,m)\) model generates a random graph such that any graph with \(m\) edges has an equal probability being selected.
In most cases, the Erdös-Rényi model refers to the first variant. A simple statistical estimator of \(p\) is \[\widehat p = \frac{1}{{n\choose 2}}\sum_{i\neq j} E_{ij},\] the fraction of existing edges. However, there is a very sad news about this model–most of the observed networks are not from Erdös-Rényi model. The first three properties partition the possible range of \(n,p\) into three categories and from this model, it is unlikely that the graph will have multiple ‘stars’ (stars refer to the vertices that have high degrees, i.e., many other nodes connecting to them). Many realistic networks such as the social networks often contain several stars.
Stochastic block model (SBM). A popular alternative to the Erdös-Rényi model is the stochastic block model. The idea is very simple. Suppose that there is a partition of the vertices that forms \(K\) groups of vertices. The SBM places an equal probability for forming a within-group edge and another equal probability for forming a between-group edge. And every edge is formed independently from each other. This probability model has \(K(K+1)/2\) parameters and the parameters form a symmetric \(K\times K\) matrix \(\theta\) such that the diagonal describes the probability of forming a within-group edge in each group and the off-diagonal parts are the probability of forming a between-group edge that corresponds to the two groups. Let \(g_i\in \{1,2,\cdots, K\}\) denotes the group that the node \(i\) belongs to. Then the stochastic block model can be written as \[P(E_{ij} = 1) = \theta_{g_i, g_j}.\] If we rearrange the vertices such that vertices are ordered with respect to the index of the nodes, you can easily see that the probability matrix \(P(E_{ij}=1)\) forms a block-diagonal structure. Often the label of which group that a vertex belongs to is unknown and the parameters are also unknown. To estimate the parameter, we can estimate the parameters by the MLE but the likelihood function is often non-convex so finding the MLE is computationally challenging. Often people use approximation approach to find a surrogate of the MLE. One common approximation is the variational approximation (also known as the variational inference/variational Bayes). See the following papers for more details:
1. Bickel, Peter, et al. “Asymptotic normality of maximum likelihood and its variational approximation for stochastic blockmodels." The Annals of Statistics 41.4 (2013): 1922-1943.
2. Celisse, Alain, Jean-Jacques Daudin, and Laurent Pierre. “Consistency of maximum-likelihood and variational estimators in the stochastic block model." Electronic Journal of Statistics 6 (2012): 1847-1899.
Another approach to estimate the parameter is via the spectral clustering algorithm. The idea is due to the fact that spectral clustering is a relaxation of the optimal graph cut problem and the optimal graph cut is a good approximation to the partition that forms the group structures in SBM. Thus, the spectral clustering gives a partition of the graph and after forming the partition, we can simply use the average within/between-group edge proportion as an estimator of the parameter. See the following for more details:
Rohe, Karl, Sourav Chatterjee, and Bin Yu. “Spectral clustering and the high-dimensional stochastic blockmodel." The Annals of Statistics 39, no. 4 (2011): 1878-1915.
Random dot product graph (RDPG). Random dot product graph assumes \[P(E_{ij}=1) = X_i^T X_j\] with \(X_i,X_j\) are supported on the ball \(\mathbb{S}_d = \{x\in\mathbb{R}^d: \|x\| = 1\}\) and \(X_1,\cdots, X_n\) are assumed to be IID from an unknown distribution \(F\) over \(\mathbb{S}_d\) and denote \(X \in \mathbb{R}^{n\times d}\) be the matrix of \(X_1,\cdots, X_n\). Interestingly, under RDPG, there exists a simple approach to recover the latent position up to rotations called the Adjacency spectral embedding (ASE). Given an observed edge matrix \(E\), let \(\Omega_d\in\mathbb{R}^{d\times d}\) be the diagonal matrix consists of the top \(d\) eigenvalues and \(U_d\in\mathbb{R}^{n\times d}\) be the corresponding eigenvector matrix. Define a matrix \(\widehat X= U_d \Omega_d^{1/2} \in \mathbb{R}^{n\times d}\). Then \(\widehat X\) is a consistent estimator of the latent position matrix \(X\) in the sense that there exists an orthogonal matrix \(Q\in\mathbb{R}^{n\times n}\) such that \[\max_i \|Q\widehat X_i - X_i \|\leq \frac{C \log n}{n}\] with a probability of at least \(1-Cn^{-2}\). This result is from the following paper:
Lyzinski, Vince, et al. “Perfect clustering for stochastic blockmodel graphs via adjacency spectral embedding." Electronic journal of statistics 8.2 (2014): 2905-2922.
For a recent review/survey on RDPG, please see
Athreya, Avanti, et al. “Statistical inference on random dot product graphs: a survey." The Journal of Machine Learning Research 18.1 (2017): 8393-8484.
Latent space model. Both SBM and RDPG are latent space models. The latent space model assumes that there is a latent space \(\mathbb{S}\subset\mathbb{R}^d\) such that every node \(V_i\) has a latent position \(X_i\in\mathbb{S}\). And the probability of forming an edge between node \(i\) and \(j\) depends on their relative position in the latent space. Namely, \[P(E_{ij}=1) = \mu(X_i, X_j).\] A simple choice is \(\mu(x,y) = \log\mbox{odds}(\alpha+\|x-y\|)\), namely, \[P(E_{ij}=1) = \log\mbox{odds}(\alpha +\|X_i-X_j\|) = \log\mbox{odds}(\alpha +A_{ij}),\] where \(\alpha\) is a parameter. The estimation of the parameter \(\alpha\) and the matrix \(A_{ij}\) are often done by the ML procedure but this could be computationally challenging. Note that we can only recover \(A_{ij}\) rather than the exact location \(X_i\) and \(X_j\) because the model will be translational and rotational invariant with respect to \(X_i\)’s. For latent space model for networks, I would recommend the first paper on this topic:
Hoff, Peter D., Adrian E. Raftery, and Mark S. Handcock. “Latent space approaches to social network analysis." Journal of the american Statistical association 97.460 (2002): 1090-1098.
Exponential family Random Graph Model (ERGM). The ERGM utilizes the exponential family in statistics to model the generating probability of a specific graph. It does not assume any independence between pairs of edges so it is a very flexible and powerful model. Recall that to describe the probability model of the network, we only need to specify the randomness of edges or the edge matrix \(E\). The ERGM models the probability of the random matrix \(E\) as \[P(E=e;\theta) \propto \exp\left(\sum_{\ell=1}^d \theta_\ell T_\ell(e)\right)=\exp\left(\theta^T T(e)\right),\] where \(e\in \{0,1\}^{n\times n}\) is a realization of the edge matrix and \(T(e) = (T_1(e),\cdots, T_d(e))\) are the sufficient statistics of the model and \(\theta = (\theta_1,\cdots,\theta_d)\) is the parameter. Sometimes, we will introduce the partition function \(Z(\theta)\) and write the above model as \[P(E=e;\theta) = \frac{1}{Z(\theta)}\exp\left(\theta^T T(e)\right),\] where \[Z(\theta) = \sum_{e} \exp\left(\theta^T T(e)\right).\]
The Erdös-Rényi model is a special case of ERGM. To see this, note that from Erdös-Rényi model, \[\begin{align*} P(E=e) = \prod_{i,j=1}^n p^{e_{ij}} (1-p)^{1-e_{ij}} &=\exp\left(\sum_{i,j} e_{ij} \log p + \sum_{i,j}(1-e_{ij})\log (1-p)\right)\\ & = \exp\left( n^2\log (1-p) + \sum_{i,j}e_{ij}\log \left(\frac{p}{1-p}\right)\right). \end{align*}\] So a sufficient statistic is \(T(e) = \sum_{i,j}e_{ij}\), the total number of edges.
In ERGM, all we need is to estimate the parameter \(\theta\). In a usual exponential family, this is often done by the MLE. Given a random adjacency matrix \(E\), you can actually show that the MLE of \(\theta_\ell\) satisfies \[T_\ell(E) = \mathbb{E}_{E'\sim P_{\widehat\theta }}(T_\ell(E')).\] The left-hand-side is a fixed sufficient statistic and all we need to do is to find the parameter such that the expected value of the sufficient statistic happens to be the same as the observed sufficient statistic. Although this seems to be not difficult (in the regular exponential family problem), it is actually a computationally challenging problem. For any parameter \(\theta\), we do not have a closed-form of \(\mathbb{E}_{E'\sim P_{\theta}}(T_\ell(E'))\) so the only way to compute the expectation is to compute the probability \(P(E';\theta)\) for each possible \(E'\) and then use the fact \[\mathbb{E}_{E'\sim P_{\theta}}(T_\ell(E')) = \sum_{e} T_\ell(e) P(E=e;\theta).\] The problem of the above summation is that there are too many terms in the summation of \(e\). For a network with \(n\) vertices, there are \(2^{{n\choose 2}}\) possible edges. So it is almost impossible to compute the sum when \(n\) is not small. There are some stochastic approximations or MCMC methods for approximating the MLE. For instance, the statnet2 is an R package developed at our department to address this issue.
Statistical inference. Recently, there is more and more attentions on making statistical inference with network models. But in general, it is hard to construct a confidence interval for a parameter of interest under a network model. And it is not so easy to construct a resampling method such as the bootstrap in the network models. However, there are some progress on the resampling inference from a network model. One possible approach is a parametric bootstrap: we assume a parametric model (such as SBMs) and resample from the fitted parametric model; see, e.g.,
Bickel, Peter, et al. “Asymptotic normality of maximum likelihood and its variational approximation for stochastic blockmodels." The Annals of Statistics 41.4 (2013): 1922-1943.
Another approach is to use the latent space model. We resample the fitted latent space positions and regenerate a new network from the resampled latent space positions. See the following paper for an example:
Levin, Keith, and Elizaveta Levina. “Bootstrapping Networks with Latent Space Structure." arXiv preprint arXiv:1907.10821 (2019).
Incorporating covariates. It is possible to incorporate the covariates into the network model. The ERGMs can easily have some effects from the covariates. For instance, if we have a network representing by an adjacency matrix \(E\) and covariates \(X\in \mathbb{R}^{n\times p}\) for each node, we may use the logistic regression model \[P(E=e|X;\beta) = \prod_{i\neq j} \left(\frac{\exp(E_{ij}(\beta_0 + \beta_1^T (X_{i} + X_j)))}{1+\exp(\beta_0 + \beta_1^T (X_{i} + X_j))}\right).\] See the following lecture note for more details:
https://www.stat.washington.edu/people/pdhoff/courses/567/Notes/l13_ergmcov.pdf
Community detection. The goal of community detection is to find communities within a network data. While there are several definitions of communities, often a community is a collection of nodes such that nodes share many inter-community connections. In a sense, a community is like a cluster of observations. The stochastic block model is a common model for modeling communities– nodes within the same block are in the same community. A recent survey on this topic can be found in
Abbe, Emmanuel. "Community detection and stochastic block models: recent developments." The Journal of Machine Learning Research 18, no. 1 (2017): 6446-6531.
Networks as covariates. The network may be used as a covariate that informs us the dependency among the response variable. A good news of this use of network is that we no longer have to worry about the probability model that generates the observed networks. For some work along this direction, please see:
1. Forastiere, Laura, Edoardo M. Airoldi, and Fabrizia Mealli. “Identification and estimation of treatment and interference effects in observational studies on networks." arXiv preprint arXiv:1609.06245 (2016).
2. Basse, Guillaume W., and Edoardo M. Airoldi. “Model-assisted design of experiments in the presence of network-correlated outcomes." Biometrika 105, no. 4 (2018): 849-858.
3. Basse, Guillaume W., and Edoardo M. Airoldi. “Limitations of design-based causal inference and A/B testing under arbitrary and network interference." Sociological Methodology 48, no. 1 (2018): 136-151.
Sampling a network data In many realistic scenarios, we do not observe the entire network but only a fraction of it. For instance, medical researchers often use ‘coupons’ for participants to recruit other participants. This generates samples from a network structure (assuming that a participant only gives the coupon to his/her friends).
Note that sometimes we can design how the network is sampled but sometimes we cannot design the sampling scheme–we already observed the network. However, even we cannot design the sampling scheme, if we have information about how the network is sampled, we can construct a corresponding probabilistic model that helps us understand properties of an estimator.
There are four common types of sampling a network data:
Node-induced subgraph sampling. We random choose \(m\) out of \(n\) nodes from the graph and examine if there are edges within these \(m\) sampled nodes. This generates a subgraph \(G'\subset G\). This is the common scenario that we recruit several participants to join a study and examine their relationships.
Edge-induced subgraph sampling. Randomly choose a set of edges from \(G\) and construct the corresponding subgraph. Note that a node will be included in the sampled subgraph if any of its edge is sampled. Again, this generates a subgraph \(G'\subset G\).
Egocentric sampling. Similar to the node-induced approach but whenever we observe a node, we also observe all its edges along with adjacent nodes. In addition, if any of the adjacent nodes are linked, we also observe this information. This occurs in questionnaires such as the ones asking
Who are your friends?
Among your friends, who are friends with each other?
Link-tracing sampling (also known as snowball sampling). Link-tracing sampling is similar to the egocentric sampling but we repeat this process several times. We perform egocentric sampling on those who are already recruited and expand the sampled network gradually.
For more details, I would recommend the following lecture note
https://www.stat.washington.edu/people/pdhoff/courses/567/Notes/l18_sampling.pdf
Recall that two random variables \(X\) and \(Y\) are independent if \[P(X\leq x, Y\leq y) = P(X\leq x) P(Y\leq y).\] In this case, we write it as \(X{\perp\mspace{-10.5mu}\perp}Y\). Let \(p_X\) and \(p_Y\) denote the PDF or PMF of \(X\) and \(Y\), respectively. Then independence also implies \[p_{XY}(x,y) = p_X(x) p_Y(y) \Leftrightarrow p_{X|Y}(x|y) =p_X(x).\]
Consider a special case where both \(X\) and \(Y\) are categorical variables such that \(X\in \{1,2,\cdots, m\}\) and \(Y\in \{1,2,\cdots, n\}\). We further define \[q_{ij} = P(X=i, Y=j)\quad q_{i+} = P(X=i)\quad q_{+j} = P(Y=j).\] Then \(X{\perp\mspace{-10.5mu}\perp}Y\) if and only if \[q_{ij} = q_{i+}\cdot q_{+j}\quad \mbox{for all $i,j$.}\]
Lemma 15. Let \(Q\) be an \(m\times n\) matrix such that \(Q_{ij} = q_{ij}\). Then \(X{\perp\mspace{-10.5mu}\perp}Y\) if and only if the matrix \(Q\) has rank \(1\).
\(\Rightarrow:\)
This direction is easy to see because \(q_{ij} = q_{i+}\cdot q_{+j}\) implies that \(Q = uv^T\), where \(u = (q_{1+},q_{2+},\cdots, q_{n+})\) and \(v=(q_{+1},q_{+2},\cdots,q_{+m})\).
\(\Leftarrow:\)
If \(Q\) has rank \(1\), there exists vectors \(u\in\mathbb{R}^n\) and \(b\in\mathbb{R}^m\) such that \(Q = uv^T\). Because \(q_{ij}\geq 0\), we may choose every elements of \(u\) and \(v\) to be non-negative, i.e., \(u_j\geq 0\) and \(v_j\geq0\) for every \(i\) and \(j\).
Since \(Q_{ij} = p_{ij}= u_i v_j\), \[p_{i+} = \sum_{j=1,\cdots, m} p_{ij} = \sum_{j=1}^mu_i v_j = u_i v_+,\] where \(v_+ = \sum_{j=1}^m v_j>0\). Similarly, \[p_{+j} = u_+ v_j,\quad u_+ = \sum_{i=1}^n u_i .\] Therefore, we obtain \[u_i = \frac{p_{i+}}{v_+},\quad v_j = \frac{p_{+j}}{u_+}\] and \[p_{ij}= u_i v_j = \frac{p_{i+}p_{+j}}{v_+u_+}= p_{i+}p_{+j}\] because \(v_+u_+ = \sum_{j=1}^m v_j\sum_{i=1}^n u_i = \sum_{i,j} u_iv_j = \sum_{i,j} p_{ij} = 1\).
For three RVs \(X,Y,\) and \(Z\), we say \(X,Y\) are conditionally independent given \(Z\) if \[P(X\leq x,Y\leq y|Z= z) = P(X\leq x|Z=z) P(Y\leq y|Z=z)\] for every \(x\) and \(y\) and \(P_Z\)-almost everywhere of \(z\). \(P_Z\)-almost everywhere of \(z\) means that the above equality holds for all \(z\) except for a set of values that has \(0\) probability. It is a slightly weaker notion than ‘for every \(z\)’. We use the notation \[X{\perp\mspace{-10.5mu}\perp}Y|Z\] to denote the case where \(X,Y\) are conditionally independent given \(Z\).
Note that \(X{\perp\mspace{-10.5mu}\perp}Y|Z\) also implies \[P(X\leq x|Y=y,Z= z) = P(X\leq x|Z=z)\] for every \(x\) and \(P_{Y,Z}\)-almost everywhere of \((y,z)\).
Theorem 16. Let \(p_{XYZ}\) be the joint PDF/PMF of \(X,Y,\) and \(Z\). Then the following are equivalent:
\(X{\perp\mspace{-10.5mu}\perp}Y|Z\).
\(p_{XY|Z}(x,y|z)= p_{X|Z}(x|z) p_{Y|Z}(y|z)\) a.e.
\(p_{X|YZ}(x|y,z)= p_{X|Z}(x|z)\) a.e.
\(p_{XYZ}(x,y,z)= \frac{p_{XZ}(x,z)p_{YZ}(y,z)}{p_Z(z)}\) a.e.
\(p_{XYZ}(x,y,z)= g(x,z) h(y,z)\), where \(g\) and \(h\) are some (measurable) functions.
\(p_{X|YZ}(x|y,z)= w(x,z)\), where \(w\) is some (measurable) function.
The equivalence between (i), (ii), (iii), and (iv) are trivial so we focus on case (v) and (vi).
(ii) \(\Rightarrow\) (v):
Because \[p_{XY|Z}(x,y|z)= p_{X|Z}(x|z)
p_{Y|Z}(y|z),\] we have \[\frac{p_{XYZ}(x,y,z)}{p_Z(z)}=
\frac{p_{XZ}(x,z)}{p_Z(z)} \frac{p_{YZ}(y,z)}{p_Z(z)}\] so \[p_{XYZ}(x,y,z) =
\frac{p_{XZ}(x,z)p_{YZ}(y,z)}{p_Z(z)} = h(x,z) g(y,z),\] which
proves (v).
(v) \(\Rightarrow\) (vi):
Based on (v), we have \[p_{YZ}(y,z) = \int
p_{XYZ}(x,y,z)dx = h(y,z) \int g(x,z)dx= h(y,z) q(z).\] Thus,
\[p_{X|YZ}(x|y,z) =
\frac{p_{XYZ}(x,y,z)}{p_{YZ}(y,z)} = \frac{g(x,z) h(y,z)}{h(y,z) q(z)} =
\frac{g(x,z)}{q(z)} = w(x,z).\]
Finally, we show that (vi) \(\Rightarrow\) (iii):
\[\begin{align*}
p_{X|Z} (x|z) = \int p_{XY|Z} (x,y |z)dy &= \int p_{X|YZ} (x|y ,z)
p_{Y|Z}(y|z)dy\\
&= w(x,z) \int p_{Y|Z}(y|z)dy = w(x,z) = p_{X|YZ} (x|y ,z).
\end{align*}\]
Here are five important properties of conditional independence. Let \(X,Y,Z,W\) be RVs.
(symmetry) \(X{\perp\mspace{-10.5mu}\perp}Y|Z\Longleftrightarrow Y{\perp\mspace{-10.5mu}\perp}X|Z\).
(decomposition) \(X
{\perp\mspace{-10.5mu}\perp}Y|Z\Longrightarrow
h(X){\perp\mspace{-10.5mu}\perp}Y|Z\) for any (measurable)
function \(h\).
A special case is: \((X,W)
{\perp\mspace{-10.5mu}\perp}Y|Z\Longrightarrow
X{\perp\mspace{-10.5mu}\perp}Y|Z\).
(weak union) \(X{\perp\mspace{-10.5mu}\perp}Y|Z\Longrightarrow
X{\perp\mspace{-10.5mu}\perp}Y|Z, h(X)\) for any (measurable)
function \(h\).
A special case is: \((X,W){\perp\mspace{-10.5mu}\perp}Y|Z
\Longrightarrow X{\perp\mspace{-10.5mu}\perp}Y|(Z,W)\)
(contraction) \[X{\perp\mspace{-10.5mu}\perp}Y|Z\mbox{ and } X{\perp\mspace{-10.5mu}\perp}W|(Y,Z)\Longleftrightarrow X{\perp\mspace{-10.5mu}\perp}(W,Y)|Z.\]
If the joint PDF \(p_{XYZW}(x,y,z,w)\) satisfies \(f_{YZW}(y,z,w)>0\) almost everywhere. Then \[X{\perp\mspace{-10.5mu}\perp}Y|(W,Z) \mbox{ and } X{\perp\mspace{-10.5mu}\perp}W|(Y,Z) \Longleftrightarrow X{\perp\mspace{-10.5mu}\perp}(W,Y)|Z.\]
You may add the intercept term into the linear SEM. Here we ignore it to obtain an elegant result.↩︎