In this paper we provide a short tutorial illustrating the new functions in the package ggm that deal with ancestral, summary and ribbonless graphs. These are mixed graphs (containing three types of edges) that are important because they capture the modified independence structure after marginalisation over, and conditioning on, nodes of directed acyclic graphs. We provide functions to verify whether a mixed graph implies that
Graphical Markov models have become a part of the mainstream of statistical theory and application in recent years. These models use graphs to represent conditional independencies among sets of random variables. Nodes of the graph correspond to random variables and edges to some type of conditional dependency.
In the literature on graphical models the two most used classes of graphs are directed acyclic graphs (DAGs) and undirected graphs. DAGs have proven useful, among other things, to specify the data generating processes when the variables satisfy an underlying partial ordering.
For instance, suppose that we have four observed variables:
This model can be represented by the DAG in Figure 1(a) with
nodes associated with the variables and edges indicating the
dependencies represented by the regression coefficients (
From the graph it is seen, for instance, that the ratio of the two blood
pressures (
A remarkable result is that the independencies can be deduced from the
graph alone, without reference to the equations, by using a criterion
called
The model has four observed variables but includes an unobserved
variable, that is, the genetic factor
The induced model is said to be obtained after marginalisation over
A mixed graph with arrows and arcs, as shown in Figure 1(b),
can be used to represent the induced independence model after
marginalisation over
The graph of Figure 1(b) belongs to a class of models called regression chain graph models. This class generalises the recursive generating process of DAGs by permitting joint responses, coupled in the graph by arcs, and thus appears to be an essential extension for applications; see Cox and N. Wermuth (1996). Regression chain graphs can be used as a conceptual framework for understanding multivariate dependencies, for example in longitudinal studies. The variables are arranged in a sequence of blocks, such that (a) all variables in one block are of equal standing and any dependence between them is represented by an arc, and (b) all variables in one block are responses to variables in all blocks to their right, so that any dependencies between them are directed, represented by an arrow pointing from right to left. The graph shows how the data analysis can be broken down into a series of regressions and informs about which variables should or should not be controlled for in each regression.
The class of regression chain graphs is not, however, stable under
marginalisation. For instance, suppose that the generating process for
the blood pressure data is defined by the more general regression chain
graph of Figure 2(a) where
Then, after marginalisation over
The previous illustrations show that when there are unobserved variables, DAG or regression chain graph models are no longer appropriate. The discussion could be extended to situations where there are some selection variables that are hidden variables that are conditioned on.
This motivates the introduction of a more general class of mixed graphs,
which contains three types of edges, denoted by lines,
There are at least three known classes of mixed graphs without self loops that remain in the same class, i.e. that are stable under marginalisation and conditioning. The largest one is that of ribbonless graphs (RGs) (Sadeghi 2012b), defined as a modification of MC-graphs (Koster 2002). Then, there is the subclass of summary graphs (SGs) (Wermuth 2011), and finally the smallest class of the ancestral graphs (AGs) (Richardson and P. Spirtes 2002).
In this paper, we focus on the implementation of four important tasks performed on the class of mixed graphs in R:
Generating different types of stable mixed graphs after marginalisation and conditioning.
Verifying whether an independency of the form
Generating a graph that induces the same independence structure as an input mixed graph such that the generated graph is maximal, i.e. each missing edge of the generated graph implies at least an independence statement.
Verifying whether two graphs are Markov equivalent, i.e. they induce the same independencies, and whether, given a graph of a specific type, there is a graph of a different type that is Markov equivalent to it.
The tasks above are illustrated by using a set of new functions introduced into the R package ggm (Marchetti, M. Drton, and K. Sadeghi 2012). In the next section we give the details of how general mixed graphs are defined. The following four sections deal with the four tasks respectively. For each task we give a brief introduction at the beginning of its corresponding section.
Some of the functions generalise previous contributions of ggm
discussed in Marchetti (2006). The ggm package has been improved and it is now
more integrated with other contributed packages related to graph theory,
such as graph (Gentleman, E. Whalen, W. Huber, and S. Falcon 2012),
igraph (Csardi and T. Nepusz 2006), and
gRbase (Dethlefsen and S. Højsgaard 2005), which
are now required for representing and plotting graphs. Specifically, in
addition to adjacency matrices, all the functions in the package now
accept graphNEL
and igraph
objects as input, as well as a new
character string representation. A more detailed list of available
packages for graphical models can be found at the CRAN Task View
gRaphical Models in R at
http://cran.r-project.org/web/views/gR.html.
For a comprehensive discussion on the ways of defining a directed
acyclic graph, see Højsgaard, D. Edwards, and S. Lauritzen (2012). A mixed graph is a more general graph type
with at most three types of edge: directed, undirected and bi-directed,
with possibly multiple edges of different types connecting two nodes. In
ggm we provide some special tools for mixed graphs that are not
present in other packages. Here we briefly illustrate some methods to
define mixed graphs and we plot them with a new function, plotGraph
,
which uses a Tk GUI for basic interactive graph manipulation.
The first method is based on a generalisation of the adjacency matrix.
The second uses a descriptive vector and is easy to use for small
graphs. The third uses a special function makeMG
that allows the
directed, undirected, and bi-directed components of a mixed graph to be
combined.
In the adjacency matrix of a mixed graph we code the three different edges with a binary indicator: 1 for directed, 10 for undirected and 100 for bi-directed edges. When there are multiple edges the codes are added.
Thus the adjacency matrix of a mixed graph
For instance consider the following general mixed graph.:
Notice that this graph is not of much interest per se, because it is not a stable graph, but it is introduced just to illustrate the structure of the adjacency matrix.
This graph can be defined by the commands
> mg <- matrix(c( 0, 101, 0, 0, 110,
100, 0, 100, 0, 1,
0, 110, 0, 1, 0,
0, 0, 1, 0, 100,
110, 0, 0, 100, 0),
5, 5, byrow = TRUE)
> N <- c("X","Y","Z","W","Q")
> dimnames(mg) <- list(N, N)
> mg
X Y Z W Q
X 0 101 0 0 110
Y 100 0 100 0 1
Z 0 110 0 1 0
W 0 0 1 0 100
Q 110 0 0 100 0
and plotted with plotGraph(mg)
.
A more convenient way of defining small mixed graphs is based on a
simple vector coding as follows. The graph is defined by a character
vector of length type, label1, label2
type
is the edge
type and label1
and label2
are the labels of the two nodes. The edge
type accepts "a"
for a directed arrow , "b"
for an arc and "l"
for
a line. Notice that isolated nodes may not be created by this method.
For example, the vector representation of the previous mixed graph is
> mgv <- c("b","X","Y","a","X","Y","l","X","Q",
"b","Q","X","a","Y","Q","b","Y","Z",
"a","Z","W","a","W","Z","b","W","Q")
Once again as in the DAG case we can use plotGraph(mgv)
to plot the
defined graph.
makeMG
Finally the adjacency matrix of a mixed graph may be built up with the
function makeMG
. This function requires three arguments dg
, ug
and
bg
, corresponding respectively to the three adjacency matrices DG
and UG
of ggm for directed and undirected
graphs respectively. Thus for the previous mixed graph we can issue the
command
> mg <- makeMG(dg = DG(Y ~ X, Z ~ W, W ~ Z),
ug = UG(~ X*Q),
bg = UG(~ Y*X + X*Q + Q*W + Y*Z))
obtaining the same adjacency matrix (up to a permutation).
There are four general classes of stable mixed graphs.
The more general class is that of ribbonless graphs: these are mixed graphs without a specific set of subgraphs called ribbons. Figure 3 below shows two examples of ribbons. The exact definition of ribbons is given in Sadeghi (2012b).
The lack of ribbons ensures that, for any RG, there is a DAG whose
independence structure, i.e. the set of all conditional independence
statements that it induces after marginalisation over, and conditioning
on, two disjoint subsets of its node set can be represented by the given
RG. This is essential, as it shows that the independence structures
corresponding to RGs are probabilistic, that is, there exists a
probability distribution
The other classes of stable graphs are further simplification of the
class of ribbonless graphs. Summary graphs have the additional property
that there are neither arrowheads pointing to lines (i.e.
Ancestral graphs have the same constraints as summary graphs plus the additional prohibition of bows, i.e. arcs with one endpoint that is an ancestor of the other endpoint; see Richardson and P. Spirtes (2002).
However, for some ribbonless and summary graphs the corresponding parametrisation is sometimes not available even in the case of a standard joint Gaussian distribution.
If we suppose that stable mixed graphs are only used to represent the independence structure after marginalisation and conditioning, we can consider all types as equally appropriate. However, each of the three types has been used in different contexts and for different purposes. RGs have been introduced in order to straightforwardly deal with the problem of finding a class of graphs that is closed under marginalisation and conditioning by a simple process of deriving them from DAGs. SGs are used when the generating DAG is known, to trace the effects in the sets of regressions as described earlier. AGs are simple graphs, meaning that they do not contain multiple edges and the lack of bows ensures that they satisfy many desirable statistical properties.
In addition, when one traces the effects in regression models with latent and selection variables (as described in the introduction) ribbonless graphs are more alerting to possible distortions (due to indirect effects) than summary graphs, and summary graphs are more alerting than ancestral graphs; see also Wermuth and D. R. Cox (2008). For the exact definition and a thorough discussion of all such graphs, see Sadeghi (2012b).
Sadeghi (2012b) also defines the algorithms for generating stable mixed graphs of a specific type for a given DAG or for a stable mixed graph of the same type after marginalisation and conditioning such that they induce the marginal and conditional DAG-independence structure. We implement these algorithms in this paper.
By “generating graphs” we mean applying the defined algorithms, e.g. those for generating stable mixed graphs to graphs, in order to generate new graphs.
Three main functions RG
, SG
, and AG
are available to generate and
plot ribbonless, summary, and ancestral graphs from DAGs, using the
algorithms in Sadeghi (2012b). These algorithms look for the paths with three
nodes and two edges in the graph whose inner nodes are being
marginalised over or conditioned on, and generate appropriate edges
between the endpoints. These have two important properties: (a) they are
well-defined in the sense that the process can be performed in any order
and will always produce the same final graph, and (b) the generated
graphs induce the modified independence structure after marginalisation
and conditioning; see Sadeghi (2012b) for more details.
The functions RG
, SG
, and AG
all have three arguments: a
, the
given input graph, M
, the marginalisation set and C
, the
conditioning set. The graph may be of class "graphNEL"
or of class
"igraph"
or may be represented by a character vector, or by an
adjacency matrix, as explained in the previous sections. The sets M
and C
(default c()
) must be disjoint vectors of node labels, and
they may possibly be empty sets. The output is always the adjacency
matrix of the generated graph. There are two additional logical
arguments showmat
and plot
to specify whether the adjacency matrix
must be explicitly printed (default TRUE
) and the graph must be
plotted (default FALSE
).
We start from a DAG defined in two ways, as an adjacency matrix and as a character vector:
> ex <- matrix(c(0,1,0,0,0,0,0,0,
0,0,1,0,0,0,0,0,
0,0,0,0,0,0,0,0,
0,0,1,0,1,0,1,0,
0,0,0,0,0,1,0,0,
0,0,0,0,0,0,0,0,
0,0,0,0,0,1,0,0,
0,0,0,0,0,1,1,0),
8, 8, byrow = TRUE)
>
> exvec <- c("a",1,2,"a",2,3,"a",4,3,
"a",4,5,"a",4,7,"a",5,6,
"a",7,6,"a",8,6,"a",8,7)
> plotGraph(ex)
Then we define two disjoint sets
> M <- c(5,8)
> C <- 3
and we generate the ribbonless, summary and ancestral graphs from the DAG with the associated plot.
> RG(ex, M, C, plot = TRUE)
1 2 4 6 7
1 0 1 0 0 0
2 0 0 10 0 0
4 0 10 0 1 1
6 0 0 0 0 100
7 0 0 0 101 0
The summary graph is also plotted:
> plotGraph(SG(ex,M,C))
1 2 4 6 7
1 0 10 0 0 0
2 10 0 10 0 0
4 0 10 0 1 1
6 0 0 0 0 100
7 0 0 0 101 0
The induced ancestral graph is obtained from the DAG defined as a vector.
> AG(exvec, M, C, showmat = FALSE, plot = TRUE)
To globally verify whether an independence statement of the form
In general, a path is said to be
The msep
. Note that there is still a
function dSep
in ggm for msep
.
The function has four arguments, where the first is the graph a
, in
one of the forms discussed before, and the other three are the disjoint
sets A
, B
, and C
.
For example, consider the DAG of Figure 1(a):
> a <- DAG(Y ~ U + Z, X ~ U + W, Z ~ W)
We see that
> msep(a, "Y", "W", "Z")
[1] TRUE
and the same statement holds for the induced ancestral graph after
marginalisation over U
:
> b <- AG(a, M = "U")
> msep(b, "Y", "W", "Z")
[1] TRUE
This was expected because the induced ancestral graph respects all the
independence statements induced by U
.
As a more complex example, consider the following summary graph,
> a <- makeMG(dg= DG(W ~ Z, Z ~ Y + X),
bg= UG(~ Y*Z))
> plotGraph(a)
Then, the two following statements verify whether X
is Y
given Z
, and whether X
is Y
(given the
empty set):
> msep(a, "X", "Y", "Z")
[1] FALSE
> msep(a, "X", "Y")
[1] TRUE
For many subclasses of graphs a missing edge corresponds to some independence statement, but for the more complex classes of mixed graphs this is not necessarily true. A graph where each of its missing edges is related to an independence statement is called a maximal graph. For a more detailed discussion on maximality of graphs and graph-theoretical conditions for maximal graphs, see Richardson and P. Spirtes (2002) and Sadeghi and S. L. Lauritzen (2012). Sadeghi and S. L. Lauritzen (2012) also gave an algorithm for generating maximal ribbonless graphs that induces the same independence structure as an input non-maximal ribbonless graph. This algorithm has been implemented in ggm as illustrated below.
Given a non-maximal graph, we can obtain the adjacency matrix of a
maximal graph that induces the same independence statements with the
function Max
. This function uses the algorithm by Sadeghi (2012a), which is an
extension of the implicit algorithm presented in Richardson and P. Spirtes (2002). The related
functions MAG
, MSG
, and MRG
, are just handy wrappers to obtain
maximal AGs, SGs and RGs, respectively. For example,
> H <- matrix(c(0 ,100, 1, 0,
100,0 ,100, 0,
0 ,100, 0,100,
0, 1 ,100, 0), 4, 4)
> plotGraph(H)
is a non-maximal ancestral graph, with the missing edge between nodes 1 and 4 that is not associated with any independence statement. Its associated maximal graph is obtained by
> Max(H)
1 2 3 4
1 0 100 0 100
2 100 0 100 1
3 1 100 0 100
4 100 0 100 0
> plotGraph(Max(H))
As the graph H
is an ancestral graph (as may be verified by the
function isAG
), we obtain the same result with
> MAG(H)
1 2 3 4
1 0 100 0 100
2 100 0 100 1
3 1 100 0 100
4 100 0 100 0
Two graphical models are said to be Markov equivalent when their
associated graphs, although nonidentical, imply the same independence
structure, that is the same set of independence statements. Thus two
Markov equivalent models cannot be distinguished on the basis of
statistical tests of independence, even for arbitrary large samples. For
instance, it is easy to verify that the two directed acyclic graphs
models
Sometimes, we can check whether graphs of different types are Markov
equivalent. For instance the DAG
Markov equivalent models may be useful in applications because (a) they may suggest alternative interpretations of a given well-fitting model or (b) on the basis of the equivalence one can choose a simpler fitting algorithm. For instance, the previous bi-directed graph model may be fitted, using the Markov equivalent DAG, in terms of a sequence of univariate regressions.
In the literature several problems related to Markov equivalences have been discussed. These include (a) verifying the Markov equivalence of given graphs, (b) presenting conditions under which a graph of a specific type can be Markov equivalent to a graph of another type, and (c) providing algorithms for generating Markov equivalent graphs of a certain type from a given graph.
The function MarkEqRcg
tests whether two regression chain graphs are
Markov equivalent. This function simply finds the skeleton and all
unshielded collider V-configurations in both graphs and tests whether
they are identical, see Wermuth and K. Sadeghi (2012). The arguments of this function are the
two graphs a
and b
in one of the allowed forms. For example,
> H1 <- makeMG(dg = DAG(W ~ X, Q ~ Z),
bg = UG(~ X*Y + Y*Z + W*Q))
> H2 <- makeMG(dg = DAG(W ~ X, Q ~ Z, Y ~ X + Z),
bg = UG(~ W*Q))
> H3 <- DAG(W ~ X, Q ~ Z + W, Y ~ X + Z)
> plotGraph(H1); plotGraph(H2); plotGraph(H3)
We can now verify Markov equivalence as follows
> MarkEqRcg(H1,H2)
[1] TRUE
> MarkEqRcg(H1,H3)
[1] FALSE
> MarkEqRcg(H2,H3)
[1] FALSE
To test Markov equivalence for maximal ancestral graphs the algorithm is
much more computationally demanding (see Ali and T. Richardson (2004)) and, for this purpose,
the function MarkEqMag
has been provided. Of course, one can use this
function for Markov equivalence of regression chain graphs (which are a
subclass of maximal ancestral graphs). For example,
> A1 <- makeMG(dg = DG(W ~ Y),
bg = UG(~ X*Y + Y*Z + Z*W))
> A2 <- makeMG(dg = DG(W ~ Y, Y ~ X),
bg = UG(~ Y*Z + Z*W))
> A3 <- makeMG(dg = DG(W ~ Y, Y ~ X, Z ~ Y),
bg = UG(~ Z*W))
> plotGraph(A1); plotGraph(A2); plotGraph(A3)
> MarkEqMag(H1,H2)
[1] TRUE
> MarkEqMag(H1,H3)
[1] FALSE
> MarkEqMag(H2,H3)
[1] FALSE
To obtain an alternative interpretation of an independence structure by
using different graphical models, it is important to verify if a given
graph is capable of being Markov equivalent to a graph of a specific
class of graphs (such as DAGs, undirected graphs, or bidirected graphs),
and if so, to obtain as a result such a graph. The functions
RepMarDAG
, RepMarUG
, and RepMarBG
do this for DAGs, undirected
graphs, and bidirected graphs, respectively. For associated conditions
and algorithms, see Sadeghi (2012a). For example, given the following graph
> H <- matrix(c( 0,10, 0, 0,
10, 0, 0, 0,
0, 1, 0,100,
0, 0,100, 0), 4, 4)
> plotGraph(H)
we can see that it is Markov equivalent to a DAG, by
> RepMarDAG(H)
$verify
[1] TRUE
$amat
1 2 3 4
1 0 1 0 0
2 0 0 1 0
3 0 0 0 0
4 0 0 1 0
> plotGraph(RepMarDAG(H))
On the other hand it is not Markov equivalent to an undirected graph or to a bidirected graph.
> RepMarUG(H)
$verify
[1] FALSE
$amat
[1] NA
> RepMarBG(H)
$verify
[1] FALSE
$amat
[1] NA
The authors are grateful to Steffen Lauritzen for helpful suggestions on codes and comments on an earlier version of the paper and to Nanny Wermuth, the editor, and referees for their insightful comments.
Bayesian, GraphicalModels, Optimization
This article is converted from a Legacy LaTeX article using the texor package. The pdf version is the official version. To report a problem with the html, refer to CONTRIBUTE on the R Journal homepage.
Text and figures are licensed under Creative Commons Attribution CC BY 4.0. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".
For attribution, please cite this work as
Sadeghi & M. Marchetti, "Graphical Markov Models with Mixed Graphs in R", The R Journal, 2012
BibTeX citation
@article{RJ-2012-015, author = {Sadeghi, Kayvan and M. Marchetti, Giovanni}, title = {Graphical Markov Models with Mixed Graphs in R}, journal = {The R Journal}, year = {2012}, note = {https://rjournal.github.io/}, volume = {4}, issue = {2}, issn = {2073-4859}, pages = {65-73} }