Package 'pcalg'

Description

Add background knowledge x -> y to an adjacency matrix and complete the orientation rules from Meek (1995).

Usage

addBgKnowledge(gInput, x = c(), y = c(), verbose = FALSE, checkInput = TRUE)

Arguments

Details

If the input is a graphNEL object, it will be converted into an adjacency matrix of type amat.cpdag. If x and y are given and if amat[y,x] != 0, this function adds orientation x -> y to the adjacency matrix amat and completes the orientation rules from Meek (1995). If x and y are not specified (or empty vectors) this function simply completes the orientation rules from Meek (1995). If x and y are vectors of length k, k>1, this function tries to add x[i] -> y[i] to the adjacency matrix amat and complete the orientation rules from Meek (1995) for every $i \in \{1, \ldots, k\}$ (see Algorithm 1 in Perkovic et. al, 2017).

Value

An adjacency matrix of type amat.cpdag of the maximally oriented pdag with added background knowledge x -> y or NULL, if the backgound knowledge is not consistent with any DAG represented by the PDAG with the adjacency matrix amat.

Author(s)

Emilija Perkovic and Markus Kalisch

References

C. Meek (1995). Causal inference and causal explanation with background knowledge, In Proceedings of UAI 1995, 403-410.

E. Perkovic, M. Kalisch and M.H. Maathuis (2017). Interpreting and using CPDAGs with background knowledge. In Proceedings of UAI 2017.

Examples

## a -- b -- c
amat <- matrix(c(0,1,0, 1,0,1, 0,1,0), 3,3)
colnames(amat) <- rownames(amat) <- letters[1:3]
## plot(as(t(amat), "graphNEL"))             
addBgKnowledge(gInput = amat) ## amat is a valid CPDAG
## b -> c is directed; a -- b is not directed by applying
## Meek's orientation rules
bg1 <- addBgKnowledge(gInput = amat, x = "b", y = "c") 
## plot(as(t(bg1), "graphNEL"))
## b -> c and b -> a are directed
bg2 <- addBgKnowledge(gInput = amat, x = c("b","b"), y = c("c","a")) 
## plot(as(t(bg2), "graphNEL")) 

## c -> b is directed; as a consequence of Meek's orientation rules,
## b -> a is directed as well
bg3 <- addBgKnowledge(gInput = amat, x = "c", y = "b") 
## plot(as(t(bg3), "graphNEL")) 

amat2 <- matrix(c(0,1,0, 1,0,1, 0,1,0), 3,3)
colnames(amat2) <- rownames(amat2) <- letters[1:3]
## new collider is inconsistent with original CPDAG; thus, NULL is returned
addBgKnowledge(gInput = amat2, x = c("c", "a"), y = c("b", "b"))

Description

This function is a wrapper for convenience to the function adjustmentSet from package dagitty.

Usage

adjustment(amat, amat.type, x, y, set.type)

Arguments

Details

If set.type is "minimal", then only minimal sufficient adjustment sets are returned. If set.type is "all", all valid adjustment sets are returned. If set.type is "canonical", a single adjustment set is returned that consists of all (possible) ancestors of x and y, minus (possible) descendants of nodes on proper causal paths. This canonical adjustment set is always valid if any valid set exists at all.

Value

If adjustment sets exist, list of length at least one (list elements might be empty vectors, if the empty set is an adjustment set). If no adjustment set exists, an empty list is returned.

Author(s)

Emilija Perkovic and Markus Kalisch ([email protected])

References

E. Perkovic, J. Textor, M. Kalisch and M.H. Maathuis (2015). A Complete Generalized Adjustment Criterion. In Proceedings of UAI 2015.

E. Perkovic, J. Textor, M. Kalisch and M.H. Maathuis (2017). Complete graphical characterization and construction of adjustment sets in Markov equivalence classes of ancestral graphs. To appear in Journal of Machine Learning Research.

B. van der Zander, M. Liskiewicz and J. Textor (2014). Constructing separators and adjustment sets in ancestral graphs. In Proceedings of UAI 2014.

See Also

gac for testing if a set satisfies the Generalized Adjustment Criterion.

Examples

## Example 4.1 in Perkovic et. al (2015), Example 2 in Perkovic et. al (2017)
mFig1 <- matrix(c(0,1,1,0,0,0, 1,0,1,1,1,0, 0,0,0,0,0,1,
                  0,1,1,0,1,1, 0,1,0,1,0,1, 0,0,0,0,0,0), 6,6)
type <- "cpdag"
x <- 3; y <- 6
## plot(as(t(mFig1), "graphNEL"))

## all
if(requireNamespace("dagitty", quietly = TRUE)) {
adjustment(amat = mFig1, amat.type = type, x = x, y = y, set.type =
"all")
}
## finds adjustment sets: (2,4), (1,2,4), (4,5), (1,4,5), (2,4,5), (1,2,4,5)

## minimal
if(requireNamespace("dagitty", quietly = TRUE)) {
adjustment(amat = mFig1, amat.type = type, x = x, y = y, set.type =
"minimal")
}
## finds adjustment sets: (2,4), (4,5), i.e., the valid sets with the fewest elements

## canonical
if(requireNamespace("dagitty", quietly = TRUE)) {
adjustment(amat = mFig1, amat.type = type, x = x, y = y, set.type =
"canonical")
}
## finds adjustment set: (1,2,4,5)

Description

Estimate an APDAG (a particular PDAG) using the aggregative greedy equivalence search (AGES) algorithm, which uses the solution path of the greedy equivalence search (GES) algorithm of Chickering (2002).

Usage

ages(data, lambda_min = 0.5 * log(nrow(data)), labels = NULL,
     fixedGaps = NULL, adaptive = c("none", "vstructures", "triples"),
     maxDegree = integer(0), verbose = FALSE, ...)

Arguments

Details

This function tries to add orientations to the essential graph (CPDAG) found by ges (ran with lambda=lambda_min). It does it aggregating several CPDAGs present in the solution path of GES. Conceptually, AGES starts with the essential graph found by GES ran with lambda = lambda_min. Then, it checks for further (compatible) orientation information in other essential graphs present in the solution path of GES, i.e., in essential graphs outputted by GES for larger penalty parameters. With compatible we mean that the aggregation process is done such that the final APDAG is still within the Markov equivalence graph represented by the essential graph found by GES in the following sense: an APDAG can always be extended to a DAG without creating new v-structures. This DAG lies in the Markov equivalence class represented by the essential graph found by GES. The algorithm is explained in detail in Eigenmann, Nandy, and Maathuis (2017).

The arguments fixedgaps and adaptive work also with AGES. However, they have not been studied in Eigenmann, Nandy, and Maathuis (2017). Using the argument fixedGaps, one can make sure that certain edges will not be present in the resulting essential graph: if the entry [i, j] of the matrix passed to fixedGaps is TRUE, there will be no edge between nodes $i$ and $j$. The argument adaptive can be used to relax the constraints encoded by fixedGaps according to a modification of GES called ARGES (adaptively restricted greedy equivalence search) which has been presented in Nandy, Hauser and Maathuis (2018):

• When adaptive = "vstructures" and the algorithm introduces a new v-structure $a \longrightarrow b \longleftarrow c$ in the forward phase, then the edge $a - c$ is removed from the list of fixed gaps, meaning that the insertion of an edge between $a$ and $c$ becomes possible even if it was forbidden by the initial matrix passed to fixedGaps.

• When adaptive = "triples" and the algorithm introduces a new unshielded triple in the forward phase (i.e., a subgraph of three nodes $a$, $b$ and $c$ where $a$ and $b$ as well as $b$ and $c$ are adjacent, but $a$ and $c$ are not), then the edge $a - c$ is removed from the list of fixed gaps.

With one of the adaptive modifications, the successive application of a skeleton estimation method and GES restricted to an estimated skeleton still gives a consistent estimator of the DAG, which is not the case without the adaptive modification.

For a detailed explanation of the GES function as well as its related object like essential graphs, we refer to the ges function.

Differences in the arguments with respect to GES: AGES uses data to initialize several scores taken as argument by GES. AGES modifies the forward and backward phases of GES performing single steps in either directions. For this reason, phase, iterate, and turning are not available arguments.

Value

ages returns a list with the following four components:

Author(s)

Marco Eigenmann ([email protected])

References

D.M. Chickering (2002). Optimal structure identification with greedy search. Journal of Machine Learning Research 3, 507–554

M.F. Eigenmann, P. Nandy, and M.H. Maathuis (2017). Structure learning of linear Gaussian structural equation models with weak edges. In Proceedings of UAI 2017

P. Nandy, A. Hauser and M.H. Maathuis (2018). High-dimensional consistency in score-based and hybrid structure learning. Annals of Statistics, to appear.

See Also

ges, EssGraph

Examples

## Example 1: ages adds correct orientations: Bar --> V6 and Bar --> V8

set.seed(77)

p <- 8
n <- 5000
## true DAG:
vars <- c("Author", "Bar", "Ctrl", "Goal", paste0("V",5:8))
gGtrue <- randomDAG(p, prob = 0.3, V = vars)
data = rmvDAG(n, gGtrue)


## Estimate the aggregated PDAG with ages
ages.fit <- ages(data = data)


## Estimate the essential graph with ges
## We specify the phases in order to have a fair comparison of the algorithms
## Without the phases specified it would be easy to find examples
## where each algorithm outperforms the other
score <- new("GaussL0penObsScore", data)
ges.fit <- ges(score, phase = c("forward","backward"), iterate = FALSE)

## Plots
if (require(Rgraphviz)) {
par(mfrow=c(1,3))
plot(ges.fit$essgraph, main="Estimated CPDAG with GES")
plot(ages.fit$essgraph, main="Estimated APDAG with AGES")
plot(gGtrue, main="TrueDAG")
}


## Example 2: ages adds correct orientations: Author --> Goal and Author --> V5

set.seed(50)

p <- 9
n <- 5000
## true DAG:
vars <- c("Author", "Bar", "Ctrl", "Goal", paste0("V",5:9))
gGtrue <- randomDAG(p, prob = 0.5, V = vars)
data = rmvDAG(n, gGtrue)


## Estimate the aggregated PDAG with ages
ages.fit <- ages(data = data)


## Estimate the essential graph with ges
## We specify the phases in order to have a fair comparison of the algorithms
## Without the phases specified it would be easy to find examples
## where each algorithm outperforms the other
score <- new("GaussL0penObsScore", data)
ges.fit <- ges(score, phase = c("forward","backward"), iterate = FALSE)

## Plots
if (require(Rgraphviz)) {
par(mfrow=c(1,3))
plot(ges.fit$essgraph, main="Estimated CPDAG with GES")
plot(ages.fit$essgraph, main="Estimated APDAG with AGES")
plot(gGtrue, main="TrueDAG")
}

## Example 3: ges and ages return the same graph

data(gmG)

data <- gmG8$x

## Estimate the aggregated PDAG with ages
ages.fit <- ages(data = data)


## Estimate the essential graph with ges
score <- new("GaussL0penObsScore", data)
ges.fit <- ges(score)


## Plots
if (require(Rgraphviz)) {
par(mfrow=c(1,3))
plot(ges.fit$essgraph, main="Estimated CPDAG with GES")
plot(ages.fit$essgraph, main="Estimated APDAG with AGES")
plot(gmG8$g, main="TrueDAG")
}

Description

Two types of adjacency matrices are used in package pcalg: Type amat.cpdag for DAGs and CPDAGs and type amat.pag for MAGs and PAGs. The required type of adjacency matrix is documented in the help files of the respective functions or classes. If in some functions more detailed information on the graph type is needed (i.e. DAG or CPDAG; MAG or PAG) this information will be passed in a separate argument (see e.g. gac and the examples below).

Note that you get (‘extract’) such adjacency matrices as (S3) objects of class "amat" via the usual as(., "<class>") coercion,

  as(from, "amat")

Arguments

Details

Adjacency matrices are integer valued square matrices with zeros on the diagonal. They can have row- and columnnames; however, most functions will work on the (integer) column positions in the adjacency matrix.

Coding for type amat.cpdag:

0:

No edge or tail

1:

Arrowhead

Note that the edgemark-code refers to the row index (as opposed adjacency matrices of type mag or pag). E.g.:

    amat[a,b] = 0  and  amat[b,a] = 1   implies a --> b.
    amat[a,b] = 1  and  amat[b,a] = 0   implies a <-- b.
    amat[a,b] = 0  and  amat[b,a] = 0   implies a     b.
    amat[a,b] = 1  and  amat[b,a] = 1   implies a --- b.

Coding for type amat.pag:

0:

No edge

1:

Circle

2:

Arrowhead

3:

Tail

Note that the edgemark-code refers to the column index (as opposed adjacency matrices of type dag or cpdag). E.g.:

  amat[a,b] = 2  and  amat[b,a] = 3   implies   a --> b.
  amat[a,b] = 3  and  amat[b,a] = 2   implies   a <-- b.
  amat[a,b] = 2  and  amat[b,a] = 2   implies   a <-> b.
  amat[a,b] = 1  and  amat[b,a] = 3   implies   a --o b.
  amat[a,b] = 0  and  amat[b,a] = 0   implies   a     b.

See Also

E.g. gac for a function which takes an adjacency matrix as input; fciAlgo for a class which has an adjacency matrix in one slot.

getGraph(x) extracts the graph-class object from x, whereas as(*, "amat") gets the corresponding adjacency matrix.

Examples

##################################################
## Function gac() takes an adjecency matrix of
## any kind as input.  In addition to that, the
## precise type of graph (DAG/CPDAG/MAG/PAG) needs
## to be passed as a different argument
##################################################
## Adjacency matrix of type 'amat.cpdag'
m1 <- matrix(c(0,1,0,1,0,0, 0,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,0,0,0), 6,6)
## more detailed information on the graph type needed by gac()
gac(m1, x=1,y=3, z=NULL, type = "dag")

## Adjacency matrix of type 'amat.cpdag'
m2 <- matrix(c(0,1,1,0,0,0, 1,0,1,1,1,0, 0,0,0,0,0,1,
               0,1,1,0,1,1, 0,1,0,1,0,1, 0,0,0,0,0,0), 6,6)
## more detailed information on the graph type needed by gac()
gac(m2, x=3, y=6, z=c(2,4), type = "cpdag")

## Adjacency matrix of type 'amat.pag'
m3 <- matrix(c(0,2,0,0, 3,0,3,3, 0,2,0,3, 0,2,2,0), 4,4)
## more detailed information on the graph type needed by gac()
mg3 <- gac(m3, x=2, y=4, z=NULL, type = "mag")
pg3 <- gac(m3, x=2, y=4, z=NULL, type = "pag")

############################################################
## as(*, "amat") returns an adjacency matrix incl. its type
############################################################
## Load predefined data
data(gmG)
n <- nrow    (gmG8$x)
V <- colnames(gmG8$x)

## define sufficient statistics
suffStat <- list(C = cor(gmG8$x), n = n)
## estimate CPDAG
skel.fit <- skeleton(suffStat, indepTest = gaussCItest,
                     alpha = 0.01, labels = V)
## Extract the "amat" [and show nicely via  'print()' method]:
as(skel.fit, "amat")

##################################################
## Function fci() returns an adjacency matrix
## of type amat.pag as one slot.
##################################################
set.seed(42)
p <- 7
## generate and draw random DAG :
myDAG <- randomDAG(p, prob = 0.4)

## find skeleton and PAG using the FCI algorithm
suffStat <- list(C = cov2cor(trueCov(myDAG)), n = 10^9)
res <- fci(suffStat, indepTest=gaussCItest,
           alpha = 0.9999, p=p, doPdsep = FALSE)
str(res)
## get the a(djacency) mat(rix)  and nicely print() it:
as(res, "amat")

##################################################
## pcAlgo object
##################################################
## Load predefined data
data(gmG)
n <- nrow    (gmG8$x)
V <- colnames(gmG8$x)

## define sufficient statistics
suffStat <- list(C = cor(gmG8$x), n = n)
## estimate CPDAG
skel.fit <- skeleton(suffStat, indepTest = gaussCItest,
                     alpha = 0.01, labels = V)
## Extract Adjacency Matrix - and print (via method 'print.amat'):
as(skel.fit, "amat")

pc.fit <- pc(suffStat, indepTest = gaussCItest,
             alpha = 0.01, labels = V)
pc.fit # (using its own print() method 'print.pcAlgo')

as(pc.fit, "amat")

Description

This function first checks if the total causal effect of one variable (x) onto another variable (y) is identifiable via the GBC, and if this is the case it explicitly gives a set of variables that satisfies the GBC with respect to x and y in the given graph.

Usage

backdoor(amat, x, y, type = "pag", max.chordal = 10, verbose=FALSE)

Arguments

Details

This function is a generalization of Pearl's backdoor criterion, see Pearl (1993), defined for directed acyclic graphs (DAGs), for single interventions and single outcome variable to more general types of graphs (CPDAGs, MAGs, and PAGs) that describe Markov equivalence classes of DAGs with and without latent variables but without selection variables. For more details see Maathuis and Colombo (2015).

The motivation to find a set W that satisfies the GBC with respect to x and y in the given graph relies on the result of the generalized backdoor adjustment:

If a set of variables W satisfies the GBC relative to x and y in the given graph, then the causal effect of x on y is identifiable and is given by

$P(Y|do(X = x)) = \sum_W P(Y|X,W) \cdot P(W).$

This result allows to write post-intervention densities (the one written using Pearl's do-calculus) using only observational densities estimated from the data.

If the input graph is a DAG (type="dag"), this function reduces to Pearl's backdoor criterion for single interventions and single outcome variable, and the parents of x in the DAG satisfy the backdoor criterion unless y is a parent of x.

If the input graph is a CPDAG C (type="cpdag"), a MAG M (type="mag"), or a PAG P (type="pag") (with both M and P not allowing selection variables), this function first checks if the total causal effect of x on y is identifiable via the GBC (see Maathuis and Colombo, 2015). If the effect is not identifiable in this way, the output is NA. Otherwise, an explicit set W that satisfies the GBC with respect to x and y in the given graph is found.

At this moment this function is not able to work with an RFCI-PAG.

It is important to note that there can be pair of nodes x and y for which there is no set W that satisfies the GBC, but the total causal effect might be identifiable via some other technique.

For the coding of the adjacency matrix see amatType.

Value

Either NA if the total causal effect is not identifiable via the GBC, or a set if the effect is identifiable via the GBC. Note that if the set W is equal to the empty set, the output is NULL.

Author(s)

Diego Colombo and Markus Kalisch ([email protected])

References

M.H. Maathuis and D. Colombo (2015). A generalized backdoor criterion. Annals of Statistics 43 1060-1088.

J. Pearl (1993). Comment: Graphical models, causality and intervention. Statistical Science 8, 266–269.

See Also

gac for the Generalized Adjustment Criterion (GAC), which is a generalization of GBC; pc for estimating a CPDAG, dag2pag and fci for estimating a PAG, and pag2magAM for estimating a MAG.

Examples

#####################################################################
##DAG
#####################################################################
## Simulate the true DAG
suppressWarnings(RNGversion("3.5.0"))
set.seed(123)
p <- 7
myDAG <- randomDAG(p, prob = 0.2) ## true DAG

## Extract the adjacency matrix of the true DAG
true.amat <- (amat <- as(myDAG, "matrix")) != 0 # TRUE/FALSE <==> 1/0
print.table(1*true.amat, zero.=".") # "visualization"

## Compute set satisfying the GBC:
backdoor(true.amat, 5, 7, type="dag")

#####################################################################
##CPDAG
#####################################################################
##################################################
## Example not identifiable
## Maathuis and Colombo (2015), Fig. 3a, p.1072
##################################################
## create the graph
p <- 5
. <- 0
amat <- rbind(c(.,.,1,1,1),
              c(.,.,1,1,1),
              c(.,.,.,1,.),
              c(.,.,.,.,1),
              c(.,.,.,.,.))
colnames(amat) <- rownames(amat) <- as.character(1:5)
V <- as.character(1:5)
edL <- vector("list",length=5)
names(edL) <- V
edL[[1]] <- list(edges=c(3,4,5),weights=c(1,1,1))
edL[[2]] <- list(edges=c(3,4,5),weights=c(1,1,1))
edL[[3]] <- list(edges=4,weights=c(1))
edL[[4]] <- list(edges=5,weights=c(1))
g <- new("graphNEL", nodes=V, edgeL=edL, edgemode="directed")

## estimate the true CPDAG
myCPDAG <- dag2cpdag(g)
## Extract the adjacency matrix of the true CPDAG
true.amat <- (as(myCPDAG, "matrix") != 0) # 1/0 <==> TRUE/FALSE

## The effect is not identifiable, in fact:
backdoor(true.amat, 3, 5, type="cpdag")


##################################################
## Example identifiable
## Maathuis and Colombo (2015), Fig. 3b, p.1072
##################################################

## create the graph
p <- 6
amat <- rbind(c(0,0,1,1,0,1), c(0,0,1,1,0,1), c(0,0,0,0,1,0),
              c(0,0,0,0,1,1), c(0,0,0,0,0,0), c(0,0,0,0,0,0))
colnames(amat) <- rownames(amat) <- as.character(1:6)
V <- as.character(1:6)
edL <- vector("list",length=6)
names(edL) <- V
edL[[1]] <- list(edges=c(3,4,6),weights=c(1,1,1))
edL[[2]] <- list(edges=c(3,4,6),weights=c(1,1,1))
edL[[3]] <- list(edges=5,weights=c(1))
edL[[4]] <- list(edges=c(5,6),weights=c(1,1))
g <- new("graphNEL", nodes=V, edgeL=edL, edgemode="directed")

## estimate the true CPDAG
myCPDAG <- dag2cpdag(g)
## Extract the adjacency matrix of the true CPDAG
true.amat <- as(myCPDAG, "matrix") != 0

## The effect is identifiable and the set satisfying GBC is:
backdoor(true.amat, 6, 3, type="cpdag")


##################################################################
##PAG
##################################################################
##################################################
## Example identifiable
## Maathuis and Colombo (2015), Fig. 5a, p.1075
##################################################

## create the graph
p <- 7
amat <- t(matrix(c(0,0,1,1,0,0,0, 0,0,1,1,0,0,0, 0,0,0,1,0,1,0,
                   0,0,0,0,0,0,1, 0,0,0,0,0,1,1, 0,0,0,0,0,0,0,
                   0,0,0,0,0,0,0),  7, 7))
colnames(amat) <- rownames(amat) <- as.character(1:7)
V <- as.character(1:7)
edL <- vector("list",length=7)
names(edL) <- V
edL[[1]] <- list(edges=c(3,4),weights=c(1,1))
edL[[2]] <- list(edges=c(3,4),weights=c(1,1))
edL[[3]] <- list(edges=c(4,6),weights=c(1,1))
edL[[4]] <- list(edges=7,weights=c(1))
edL[[5]] <- list(edges=c(6,7),weights=c(1,1))
g <- new("graphNEL", nodes=V, edgeL=edL, edgemode="directed")
L <- 5

## compute the true covariance matrix of g
cov.mat <- trueCov(g)

## transform covariance matrix into a correlation matrix
true.corr <- cov2cor(cov.mat)
suffStat <- list(C=true.corr, n=10^9)
indepTest <- gaussCItest

## estimate the true PAG
true.pag <- dag2pag(suffStat, indepTest, g, L, alpha = 0.9999)@amat

## The effect is identifiable  and the backdoor set is:
backdoor(true.pag, 3, 5, type="pag")

Description

This function is DEPRECATED! Use ida instead.

Usage

beta.special(dat=NA, x.pos, y.pos, verbose=0, a=0.01, myDAG=NA,
             myplot=FALSE, perfect=FALSE, method="local", collTest=TRUE,
             pcObj=NA, all.dags=NA, u2pd="rand")

Arguments

Value

estimates of intervention effects

Author(s)

Markus Kalisch ([email protected])

See Also

pcAlgo, dag2cpdag; beta.special.pcObj for a fast version of beta.special(), using a precomputed pc-object.

Description

This function is DEPRECATED! Use ida or idaFast instead.

Usage

beta.special.pcObj(x.pos, y.pos, pcObj, mcov=NA, amat=NA, amatSkel=NA, t.amat=NA)

Arguments

Value

estimates of intervention effects

Author(s)

Markus Kalisch ([email protected])

See Also

pcAlgo, dag2cpdag, beta.special

Description

$G^2$ test for (conditional) independence of binary variables $X$ and $Y$ given the (possibly empty) set of binary variables $S$.

binCItest() is a wrapper of gSquareBin(), to be easily used in skeleton, pc and fci.

Usage

gSquareBin(x, y, S, dm, adaptDF = FALSE, n.min = 10*df, verbose = FALSE)
binCItest (x, y, S, suffStat)

Arguments

Details

The $G^2$ statistic is used to test for (conditional) independence of X and Y given a set S (can be NULL). This function is a specialized version of gSquareDis which is for discrete variables with more than two levels.

Value

The p-value of the test.

Author(s)

Nicoletta Andri and Markus Kalisch ([email protected])

References

R.E. Neapolitan (2004). Learning Bayesian Networks. Prentice Hall Series in Artificial Intelligence. Chapter 10.3.1

See Also

gSquareDis for a (conditional) independence test for discrete variables with more than two levels.

dsepTest, gaussCItest and disCItest for similar functions for a d-separation oracle, a conditional independence test for Gaussian variables and a conditional independence test for discrete variables, respectively.

skeleton, pc or fci which need a testing function such as binCItest.

Examples

n <- 100
set.seed(123)
## Simulate *independent data of {0,1}-variables:
x <- rbinom(n, 1, pr=1/2)
y <- rbinom(n, 1, pr=1/2)
z <- rbinom(n, 1, pr=1/2)
dat <- cbind(x,y,z)

binCItest(1,3,2, list(dm = dat, adaptDF = FALSE)) # 0.36, not signif.
binCItest(1,3,2, list(dm = dat, adaptDF = TRUE )) # the same, here

## Simulate data from a chain of 3 variables: x1 -> x2 -> x3
set.seed(12)
b0 <- 0
b1 <- 1
b2 <- 1
n <- 10000
x1 <- rbinom(n, size=1, prob=1/2) ## = sample(c(0,1), n, replace=TRUE)

## NB:  plogis(u) := "expit(u)" := exp(u) / (1 + exp(u))
p2 <- plogis(b0 + b1*x1) ; x2 <- rbinom(n, 1, prob = p2) # {0,1}
p3 <- plogis(b0 + b2*x2) ; x3 <- rbinom(n, 1, prob = p2) # {0,1}

ftable(xtabs(~ x1+x2+x3))
dat <- cbind(x1,x2,x3)

## Test marginal and conditional independencies
gSquareBin(3,1,NULL,dat, verbose=TRUE)
gSquareBin(3,1, 2,  dat)
gSquareBin(1,3, 2,  dat) # the same
gSquareBin(1,3, 2,  dat, adaptDF=TRUE, verbose = 2)

Description

For each subset of nbrsA and nbrsC where a and c are conditionally independent, it is checked if b is in the conditioning set.

Usage

checkTriple(a, b, c, nbrsA, nbrsC,
            sepsetA, sepsetC,
            suffStat, indepTest, alpha, version.unf = c(NA, NA),
            maj.rule = FALSE, verbose = FALSE)

Arguments

Details

This function is used in the conservative versions of structure learning algorithms.

Value

Author(s)

Markus Kalisch ([email protected]) and Diego Colombo.

References

D. Colombo and M.H. Maathuis (2014).Order-independent constraint-based causal structure learning. Journal of Machine Learning Research 15 3741-3782.

Examples

##################################################
## Using Gaussian Data
##################################################
## Load predefined data
data(gmG)
n <- nrow    (gmG8$x)
V <- colnames(gmG8$x)


## define independence test (partial correlations), and test level
indepTest <- gaussCItest
alpha <- 0.01
## define sufficient statistics
suffStat <- list(C = cor(gmG8$x), n = n)

## estimate CPDAG
pc.fit <- pc(suffStat, indepTest, alpha=alpha, labels = V, verbose = TRUE)

if (require(Rgraphviz)) {
  ## show estimated CPDAG
  par(mfrow=c(1,2))
  plot(pc.fit, main = "Estimated CPDAG")
  plot(gmG8$g, main = "True DAG")
}

a <- 6
b <- 1
c <- 8
checkTriple(a, b, c,
            nbrsA = c(1,5,7),
            nbrsC = c(1,5),
            sepsetA = pc.fit@sepset[[a]][[c]],
            sepsetC = pc.fit@sepset[[c]][[a]],
            suffStat=suffStat, indepTest=indepTest, alpha=alpha,
            version.unf = c(2,2),
            verbose = TRUE) -> ct
str(ct)
## List of 4
## $ decision: int 2
## $ version : int 1
## $ SepsetA : int [1:2] 1 5
## $ SepsetC : int 1

checkTriple(a, b, c,
            nbrsA = c(1,5,7),
            nbrsC = c(1,5),
            sepsetA = pc.fit@sepset[[a]][[c]],
            sepsetC = pc.fit@sepset[[c]][[a]],
            version.unf = c(1,1),
            suffStat=suffStat, indepTest=indepTest, alpha=alpha) -> c2
stopifnot(identical(ct, c2)) ## in this case,  'version.unf' had no effect

Description

Compares the true undirected graph with an estimated undirected graph in terms of True Positive Rate (TPR), False Positive Rate (FPR) and True Discovery Rate (TDR).

Usage

compareGraphs(gl, gt)

Arguments

Details

If the input graph is directed, the directions are omitted. Special cases:

• If the true graph contains no edges, the tpr is defined to be zero.

• Similarly, if the true graph contains no gaps, the fpr is defined to be one.

• If there are no edges in the true graph and there are none in the estimated graph, tdr is one. If there are none in the true graph but there are some in the estimated graph, tdr is zero.

Value

A named numeric vector with three numbers:

Author(s)

Markus Kalisch ([email protected]) and Martin Maechler

See Also

randomDAG for generating a random DAG.

Examples

## generate a graph with 4 nodes
V <- LETTERS[1:4]
edL2 <- vector("list", length=4)
names(edL2) <- V
edL2[[1]] <- list(edges= 2)
edL2[[2]] <- list(edges= c(1,3,4))
edL2[[3]] <- list(edges= c(2,4))
edL2[[4]] <- list(edges= c(2,3))
gt <- new("graphNEL", nodes=V, edgeL=edL2, edgemode="undirected")

## change graph
gl <- graph::addEdge("A","C", gt,1)

## compare the two graphs
if (require(Rgraphviz)) {
par(mfrow=c(2,1))
plot(gt) ; title("True graph")
plot(gl) ; title("Estimated graph")
(cg <- compareGraphs(gl,gt))

}

Description

Using Fisher's z-transformation of the partial correlation, test for zero partial correlation of sets of normally / Gaussian distributed random variables.

The gaussCItest() function, using zStat() to test for (conditional) independence between gaussian random variables, with an interface that can easily be used in skeleton, pc and fci.

Usage

condIndFisherZ(x, y, S, C, n, cutoff, verbose= )
zStat         (x, y, S, C, n)
gaussCItest   (x, y, S, suffStat)

Arguments

Details

For gaussian random variables and after performing Fisher's z-transformation of the partial correlation, the test statistic zStat() is (asymptotically for large enough n) standard normally distributed.

Partial correlation is tested in a two-sided hypothesis test, i.e., basically, condIndFisherZ(*) == abs(zStat(*)) > qnorm(1 - alpha/2). In a multivariate normal distribution, zero partial correlation is equivalent to conditional independence.

Value

zStat() gives a number

$Z = \sqrt{n - \left|S\right| - 3} \cdot \log((1+r)/(1-r))/2$

which is asymptotically normally distributed under the null hypothesis of correlation 0.

condIndFisherZ() returns a logical $L$ indicating whether the “partial correlation of x and y given S is zero” could not be rejected on the given significance level. More intuitively and for multivariate normal data, this means: If TRUE then it seems plausible, that x and y are conditionally independent given S. If FALSE then there was strong evidence found against this conditional independence statement.

gaussCItest() returns the p-value of the test.

Author(s)

Markus Kalisch ([email protected]) and Martin Maechler

References

M. Kalisch and P. Buehlmann (2007). Estimating high-dimensional directed acyclic graphs with the PC-algorithm. JMLR 8 613-636.

See Also

pcorOrder for computing a partial correlation given the correlation matrix in a recursive way.

dsepTest, disCItest and binCItest for similar functions for a d-separation oracle, a conditional independence test for discrete variables and a conditional independence test for binary variables, respectively.

Examples

set.seed(42)
## Generate four independent normal random variables
n <- 20
data <- matrix(rnorm(n*4),n,4)
## Compute corresponding correlation matrix
corMatrix <- cor(data)
## Test, whether variable 1 (col 1) and variable 2 (col 2) are
## independent given variable 3 (col 3) and variable 4 (col 4) on 0.05
## significance level
x <- 1
y <- 2
S <- c(3,4)
n <- 20
alpha <- 0.05
cutoff <- qnorm(1-alpha/2)
(b1 <- condIndFisherZ(x,y,S,corMatrix,n,cutoff))
   # -> 1 and 2 seem to be conditionally independent given 3,4

## Now an example with conditional dependence
data <- matrix(rnorm(n*3),n,3)
data[,3] <- 2*data[,1]
corMatrix <- cor(data)
(b2 <- condIndFisherZ(1,3,2,corMatrix,n,cutoff))
   # -> 1 and 3 seem to be conditionally dependent given 2

## simulate another dep.case: x -> y -> z
set.seed(29)
x <- rnorm(100)
y <- 3*x + rnorm(100)
z <- 2*y + rnorm(100)
dat <- cbind(x,y,z)

## analyze data
suffStat <- list(C = cor(dat), n = nrow(dat))
gaussCItest(1,3,NULL, suffStat) ## dependent [highly signif.]
gaussCItest(1,3,  2,  suffStat) ## independent | S

Description

Computes the correlation graph. This is the graph in which an edge is drawn between node i and node j, if the null hypothesis “Correlation between $X_i$ and $X_j$ is zero” can be rejected at the given significance level $\alpha (alpha)$.

Usage

corGraph(dm, alpha=0.05, Cmethod="pearson")

Arguments

Value

Undirected correlation graph, a graph-class object (package graph); getGraph for the “fitted” graph.

Author(s)

Markus Kalisch ([email protected]) and Martin Maechler

Examples

## create correlated samples
x1 <- rnorm(100)
x2 <- rnorm(100)
mat <- cbind(x1,x2, x3 = x1+x2)

if (require(Rgraphviz)) {
## ``analyze the data''
(g <- corGraph(mat)) # a 'graphNEL' graph, undirected
plot(g) # ==> (1) and (2) are each linked to (3)

## use different significance level and different method
(g2 <- corGraph(mat, alpha=0.01, Cmethod="kendall"))
plot(g2) ## same edges as 'g'
}

Description

Convert a DAG (Directed Acyclic Graph) to a Completed Partially Directed Acyclic Graph (CPDAG).

Usage

dag2cpdag(g)

Arguments

Details

This function converts a DAG into its corresponding (unique) CPDAG as follows. Because every DAG in the Markov equivalence class described by a CPDAG shares the same skeleton and the same v-structures, this function takes the skeleton and the v-structures of the given DAG g. Afterwards it simply uses the 3 orientation rules of the PC algorithm (see references) to orient as many of the remaining undirected edges as possible.

The function is a simple wrapper function for dag2essgraph which is more powerfull since it also allows the calculation of the Markov equivalence class in the presence of interventional data.

The output of this function is exactly the same as the one using

pc(suffStat, indepTest, alpha, labels)

using the true correlation matrix in the function gaussCItest with a large virtual sample size and a large alpha, but it is much faster.

Value

A graph object containing the CPDAG.

Author(s)

Markus Kalisch ([email protected]) and Alain Hauser([email protected])

References

C. Meek (1995). Causal inference and causal explanation with background knowledge. In Proceedings of the Eleventh Conference on Uncertainty in Artificial Intelligence (UAI-95), pp. 403-411. Morgan Kaufmann Publishers, Inc.

P. Spirtes, C. Glymour and R. Scheines (2000) Causation, Prediction, and Search, 2nd edition, The MIT Press.

See Also

dag2essgraph, randomDAG, pc

Examples

## A -> B <- C
am1 <- matrix(c(0,1,0, 0,0,0, 0,1,0), 3,3)
colnames(am1) <- rownames(am1) <- LETTERS[1:3]
g1 <- as(t(am1), "graphNEL") ## convert to graph
cpdag1 <- dag2cpdag(g1)

if(requireNamespace("Rgraphviz")) {
    par(mfrow = c(1,2))
    plot(g1)
    plot(cpdag1)
}

## A -> B -> C
am2 <- matrix(c(0,1,0, 0,0,1, 0,0,0), 3,3)
colnames(am2) <- rownames(am2) <- LETTERS[1:3]
g2 <- as(t(am2), "graphNEL") ## convert to graph
cpdag2 <- dag2cpdag(g2)

if(requireNamespace("Rgraphviz")) {
    par(mfrow = c(1,2))
    plot(g2)
    plot(cpdag2)
}

Description

Convert a DAG to an (interventional or observational) essential graph.

Usage

dag2essgraph(dag, targets = list(integer(0)))

Arguments

Details

This function converts a DAG to its corresponding (interventional or observational) essential graph, using the algorithm of Hauser and Bühlmann (2012).

The essential graph is a partially directed graph that represents the (interventional or observational) Markov equivalence class of a DAG. It has the same has the same skeleton as the DAG; a directed edge represents an arrow that has a common orientation in all representatives of the (interventional or observational) Markov equivalence class, whereas an undirected edge represents an arrow that has different orientations in different representatives of the equivalence class. In the observational case, the essential graph is also known as “CPDAG” (Spirtes et al., 2000).

In a purely observational setting (i.e., if targets = list(integer(0))), the function yields the same graph as dag2cpdag.

Value

Depending on the class of dag, the essential graph is returned as

• an instance of graphNEL, if dag is an instance of graphNEL,

• an instance of EssGraph, if dag is an instance of a class derived from ParDAG.

Author(s)

Alain Hauser ([email protected])

References

A. Hauser and P. Bühlmann (2012). Characterization and greedy learning of interventional Markov equivalence classes of directed acyclic graphs. Journal of Machine Learning Research 13, 2409–2464.

P. Spirtes, C.N. Glymour, and R. Scheines (2000). Causation, Prediction, and Search, MIT Press, Cambridge (MA).

See Also

dag2cpdag, Score, EssGraph

Examples

p <- 10     # Number of random variables
s <- 0.4    # Sparseness of the DAG

## Generate a random DAG
set.seed(42)
require(graph)
dag <- randomDAG(p, s)
nodes(dag) <- sprintf("V%d", 1:p)

## Calculate observational essential graph
res.obs <- dag2essgraph(dag)

## Different argument classes
res2 <- dag2essgraph(as(dag, "GaussParDAG"))
str(res2)

## Calculate interventional essential graph for intervention targets
## {1} and {3}
res.int <- dag2essgraph(dag, as.list(c(1, 3)))

Description

Convert a DAG with latent variables into its corresponding (unique) Partial Ancestral Graph (PAG).

Usage

dag2pag(suffStat, indepTest, graph, L, alpha, rules = rep(TRUE,10),
        verbose = FALSE)

Arguments

Details

This function converts a DAG (graph object) with latent variables into its corresponding (unique) PAG, an fciAlgo class object, using the ancestor information and conditional independence tests entailed in the true DAG. The output of this function is exactly the same as the one using

fci(suffStat, gaussCItest, p, alpha, rules = rep(TRUE, 10))

using the true correlation matrix in gaussCItest() with a large “virtual sample size” and a large alpha, but it is much faster, see the example.

Value

An object of class fciAlgo, containing the estimated graph (in the form of an adjacency matrix with various possible edge marks), the conditioning sets that lead to edge removals (sepset) and several other parameters.

Author(s)

Diego Colombo and Markus Kalisch [email protected].

References

Richardson, T. and Spirtes, P. (2002). Ancestral graph Markov models. Ann. Statist. 30, 962–1030; Theorem 4.2., page 983.

See Also

fci, pc

Examples

## create the graph
set.seed(78)
g <- randomDAG(10, prob = 0.25)
graph::nodes(g) # "1" "2" ... "10" % FIXME: should be kept in result!

## define nodes 2 and 6 to be latent variables
L <- c(2,6)

## compute the true covariance matrix of g
cov.mat <- trueCov(g)
## transform covariance matrix into a correlation matrix
true.corr <- cov2cor(cov.mat)

## Find PAG
## as dependence "oracle", we use the true correlation matrix in
## gaussCItest() with a large "virtual sample size" and a large alpha:
system.time(
true.pag <- dag2pag(suffStat = list(C = true.corr, n = 10^9),
                    indepTest = gaussCItest,
                    graph=g, L=L, alpha = 0.9999) )

### ---- Find PAG using fci-function --------------------------

## From trueCov(g), delete rows and columns belonging to latent variable L
true.cov1 <- cov.mat[-L,-L]
## transform covariance matrix into a correlation matrix
true.corr1 <- cov2cor(true.cov1)

## Find PAG with FCI algorithm
## as dependence "oracle", we use the true correlation matrix in
## gaussCItest() with a large "virtual sample size" and a large alpha:
system.time(
true.pag1 <- fci(suffStat = list(C = true.corr1, n = 10^9),
                 indepTest = gaussCItest,
                 p = ncol(true.corr1), alpha = 0.9999) )

## confirm that the outputs are equal
stopifnot(true.pag@amat == true.pag1@amat)

Description

$G^2$ test for (conditional) independence of discrete (each with a finite number of “levels”) variables $X$ and $Y$ given the (possibly empty) set of discrete variables $S$.

disCItest() is a wrapper of gSquareDis(), to be easily used in skeleton, pc and fci.

Usage

gSquareDis(x, y, S, dm, nlev, adaptDF = FALSE, n.min = 10*df, verbose = FALSE)
disCItest (x, y, S, suffStat)

Arguments

Details

The $G^2$ statistic is used to test for (conditional) independence of X and Y given a set S (can be NULL). If only binary variables are involved, gSquareBin is a specialized (a bit more efficient) alternative to gSquareDis().

Value

The p-value of the test.

Author(s)

Nicoletta Andri and Markus Kalisch ([email protected]).

References

R.E. Neapolitan (2004). Learning Bayesian Networks. Prentice Hall Series in Artificial Intelligence. Chapter 10.3.1

See Also

gSquareBin for a (conditional) independence test for binary variables.

dsepTest, gaussCItest and binCItest for similar functions for a d-separation oracle, a conditional independence test for gaussian variables and a conditional independence test for binary variables, respectively.

Examples

## Simulate data
n <- 100
set.seed(123)
x <- sample(0:2, n, TRUE) ## three levels
y <- sample(0:3, n, TRUE) ## four levels
z <- sample(0:1, n, TRUE) ## two levels
dat <- cbind(x,y,z)

## Analyze data
gSquareDis(1,3, S=2, dat, nlev = c(3,4,2)) # but nlev is optional:
gSquareDis(1,3, S=2, dat, verbose=TRUE, adaptDF=TRUE)
## with too little data, gives a warning (and p-value 1):
gSquareDis(1,3, S=2, dat[1:60,], nlev = c(3,4,2))

suffStat <- list(dm = dat, nlev = c(3,4,2), adaptDF = FALSE)
disCItest(1,3,2,suffStat)

Description

Let x and y be two distinct vertices in a mixed graph G. This function computes D-SEP(x,y,G), which is defined as follows:

A node v is in D-SEP(x,y,G) iff v is not equal to x and there is a collider path between x and v in G such that every vertex on this path is an ancestor of x or y in G.

See p.136 of Sprirtes et al (2000) or Definition 4.1 of Maathuis and Colombo (2015).

Usage

dreach(x, y, amat, verbose = FALSE)

Arguments

Value

Vector of column positions indicating the nodes in D-SEP(x,y,G).

Author(s)

Diego Colombo and Markus Kalisch ([email protected])

References

P. Spirtes, C. Glymour and R. Scheines (2000) Causation, Prediction, and Search, 2nd edition, The MIT Press.

M.H. Maathuis and D. Colombo (2015). A generalized back-door criterion. Annals of Statistics 43 1060-1088.

See Also

backdoor uses this function; pag2magAM.

Description

This function tests for d-separation of nodes in a DAG.

Usage

dsep(a, b, S=NULL, g, john.pairs = NULL)

Arguments

Details

This function checks separation in the moralized graph as explained in Lauritzen (2004).

Value

TRUE if a and b are d-separated by S in G, otherwise FALSE.

Author(s)

Markus Kalisch ([email protected])

References

S.L. Lauritzen (2004), Graphical Models, Oxford University Press, Chapter 3.2.2

See Also

dsepTest for a wrapper of this function that can easily be included into skeleton, pc, fci or fciPlus. dsepAM for a similar function for MAGs.

Examples

## generate random DAG
p <- 8
set.seed(45)
myDAG <- randomDAG(p, prob = 0.3)
if (require(Rgraphviz)) {
plot(myDAG)
}

## Examples for d-separation
dsep("1","7",NULL,myDAG)
dsep("4","5",NULL,myDAG)
dsep("4","5","2",myDAG)
dsep("4","5",c("2","3"),myDAG)

## Examples for d-connection
dsep("1","3",NULL,myDAG)
dsep("1","6","3",myDAG)
dsep("4","5","8",myDAG)

Description

This function tests for d-separation (also known as m-separation) of nodes X and nodes Y given nodes S in a MAG.

Usage

dsepAM(X, Y, S = NULL, amat, verbose=FALSE)

Arguments

Details

This function checks separation in the moralized graph as explained in Richardson and Spirtes (2002).

Value

TRUE if X and Y are d-separated by S in the MAG encoded by amat, otherwise FALSE.

Author(s)

Markus Kalisch ([email protected]), Joris Mooij

References

T.S. Richardson and P. Spirtes (2002). Ancestral graph Markov models. Annals of Statistics 30 962-1030.

See Also

dsepAMTest for a wrapper of this function that can easily be included into skeleton, fci or fciPlus. dsep for a similar function for DAGs.

Examples

# Y-structure MAG
# Encode as adjacency matrix
p <- 4 # total number of variables
V <- c("X1","X2","X3","X4") # variable labels
# amat[i,j] = 0 iff no edge btw i,j
# amat[i,j] = 1 iff i *-o j
# amat[i,j] = 2 iff i *-> j
# amat[i,j] = 3 iff i *-- j
amat <- rbind(c(0,0,2,0),
              c(0,0,2,0),
              c(3,3,0,2),
              c(0,0,3,0))
rownames(amat)<-V
colnames(amat)<-V

## d-separated
cat('X1 d-separated from X2? ', dsepAM(1,2,S=NULL,amat),'\n')
## not d-separated given node 3
cat('X1 d-separated from X2 given X4? ', dsepAM(1,2,S=4,amat),'\n')
## not d-separated by node 3 and 4
cat('X1 d-separated from X2 given X3 and X4? ', dsepAM(1,2,S=c(3,4),amat),'\n')

Description

This function tests for d-separation (also known as m-separation) of node x and node y given nodes S in a MAG. dsepAMTest() is written to be easily used in skeleton, fci, fciPlus.

Usage

dsepAMTest(x, y, S = NULL, suffStat)

Arguments

Details

The function is a wrapper for dsepAM, which checks separation in the moralized graph as explained in Richardson and Spirtes (2002).

Value

Returns 1 if x and y are d-separated by S in the MAG encoded by amat, otherwise 0.

This is analogous to the p-value of an ideal (without sampling error) conditional independence test on any distribution that is faithful to the MAG.

Author(s)

Markus Kalisch ([email protected]), Joris Mooij

References

T.S. Richardson and P. Spirtes (2002). Ancestral graph Markov models. Annals of Statistics 30 962-1030.

See Also

dsepTest for a similar function for DAGs. gaussCItest, disCItest and binCItest for similar functions for a conditional independence test for gaussian, discrete and binary variables, respectively.

Examples

# Y-structure MAG
# Encode as adjacency matrix
p <- 4 # total number of variables
V <- c("X1","X2","X3","X4") # variable labels
# amat[i,j] = 0 iff no edge btw i,j
# amat[i,j] = 1 iff i *-o j
# amat[i,j] = 2 iff i *-> j
# amat[i,j] = 3 iff i *-- j
amat <- rbind(c(0,0,2,0),
              c(0,0,2,0),
              c(3,3,0,2),
              c(0,0,3,0))
rownames(amat)<-V
colnames(amat)<-V

suffStat<-list(g=amat,verbose=FALSE)
## d-separated
cat('X1 d-separated from X2? ', dsepAMTest(1,2,S=NULL,suffStat),'\n')
## not d-separated given node 3
cat('X1 d-separated from X2 given X4? ',
dsepAMTest(1,2,S=4,suffStat),'\n')
## not d-separated by node 3 and 4
cat('X1 d-separated from X2 given X3 and X4? ', dsepAMTest(1,2,S=c(3,4),
suffStat),'\n')

# Derive PAG that represents the Markov equivalence class of the MAG with the FCI algorithm
# Make use of d-separation oracle as "independence test"
indepTest <- dsepAMTest
fci.pag <- fci(suffStat,indepTest,alpha = 0.5,labels = V,verbose=FALSE)
cat('True MAG:\n')
print(amat)
cat('PAG output by FCI:\n')
print(fci.pag@amat)

Description

Tests for d-separation of nodes in a DAG. dsepTest() is written to be easily used in skeleton, pc, fci.

Usage

dsepTest(x, y, S=NULL, suffStat)

Arguments

Details

The function is based on dsep. For details on d-separation see the reference Lauritzen (2004).

Value

If x and y are d-separated by S in DAG G the result is 1, otherwise it is 0. This is analogous to the p-value of an ideal (without sampling error) conditional independence test on any distribution that is faithful to the DAG G.

Author(s)

Markus Kalisch ([email protected])

References

S.L. Lauritzen (2004), Graphical Models, Oxford University Press.

See Also

dsepAMTest for a similar function for MAGs. gaussCItest, disCItest and binCItest for similar functions for a conditional independence test for gaussian, discrete and binary variables, respectively.

Examples

p <- 8
set.seed(45)
myDAG <- randomDAG(p, prob = 0.3)

if (require(Rgraphviz)) {
## plot the DAG
plot(myDAG, main = "randomDAG(10, prob = 0.2)")
}

## define sufficient statistics (d-separation oracle)
suffStat <- list(g = myDAG, jp = RBGL::johnson.all.pairs.sp(myDAG))

dsepTest(1,6, S= NULL,  suffStat) ## not d-separated
dsepTest(1,6, S= 3,     suffStat) ## not d-separated by node 3
dsepTest(1,6, S= c(3,4),suffStat) ## d-separated by node 3 and 4

Description

This class represents an (observentional or interventional) essential graph.

Details

An observational or interventional Markov equivalence class of DAGs can be uniquely represented by a partially directed graph, the essential graph. Its edges have the following interpretation:

1. a directed edge $a \longrightarrow b$ stands for an arrow that has the same orientation in all representatives of the Markov equivalence class;

2. an undirected edge $a - b$ stands for an arrow that is oriented in one way in some representatives of the equivalence class and in the other way in other representatives of the equivalence class.

Extends

All reference classes extend and inherit methods from "envRefClass".

Constructor

new("EssGraph", nodes, in.edges, ...)

nodes

Vector of node names; cf. also field .nodes.

in.edges

A list of length p consisting of index vectors indicating the edges pointing into the nodes of the DAG.

Fields

.nodes:

Vector of node names; defaults to as.character(1:p), where p denotes the number of nodes (variables) of the model.

.in.edges:

A list of length p consisting of index vectors indicating the edges pointing into the nodes of the DAG.

targets

List of mutually exclusive intervention targets with respect to which Markov equivalence is defined.

score:

Object of class Score; used internally for score-based causal inference.

Class-Based Methods

Most class-based methods are only for internal use. Methods of interest for the user are:

repr():

Yields a representative causal model of the equivalence class, an object of a class derived from Score. Since the representative is not only characterized by the DAG, but also by appropriate parameters, the field score must be assigned for this method to work. The DAG is drawn at random; note that all representatives are statistically indistinguishable under a given set of intervention targets.

node.count():

Yields the number of nodes of the essential graph.

edge.count():

Yields the number of edges of the essential graph. Note that unoriented edges count as 2, whereas oriented edges count as 1 due to the internal representation.

Methods

plot

signature(x = "EssGraph", y = "ANY"): plots the essential graph. In the plot, undirected and bidirected edges are equivalent.

Author(s)

Alain Hauser ([email protected])

See Also

ParDAG

Examples

showClass("EssGraph")

Description

Estimate a Partial Ancestral Graph (PAG) from observational data, using the FCI (Fast Causal Inference) algorithm, or from a combination of data from different (e.g., observational and interventional) contexts, using the FCI-JCI (Joint Causal Inference) extension.

Usage

fci(suffStat, indepTest, alpha, labels, p,
    skel.method = c("stable", "original", "stable.fast"),
    type = c("normal", "anytime", "adaptive"),
    fixedGaps = NULL, fixedEdges = NULL,
    NAdelete = TRUE, m.max = Inf, pdsep.max = Inf,
    rules = rep(TRUE, 10), doPdsep = TRUE, biCC = FALSE,
    conservative = FALSE, maj.rule = FALSE,
    numCores = 1, selectionBias = TRUE,
    jci = c("0","1","12","123"), contextVars = NULL, 
    verbose = FALSE)

Arguments

Details

This function is a generalization of the PC algorithm (see pc), in the sense that it allows arbitrarily many latent and selection variables. Under the assumption that the data are faithful to a DAG that includes all latent and selection variables, the FCI algorithm (Fast Causal Inference algorithm) (Spirtes, Glymour and Scheines, 2000) estimates the Markov equivalence class of MAGs that describe the conditional independence relationships between the observed variables. Under the assumption that the data are $\sigma$-faithful to a simple (possibly cyclic) SCM that allows for latent confounding (but selection bias is absent), the FCI algorithm estimates the $\sigma$-Markov equivalence class of the DMGs (directed mixed graphs) that describe the causal relations between the observed variables and the conditional independence relationships in the observed distribution through $\sigma$-separation (Mooij and Claassen, 2020). The FCI-JCI (Joint Causal Inference) extension allows the algorithm to combine data from different contexts, for example, observational and different types of interventional data (Mooij et al., 2020).

FCI estimates a partial ancestral graph (PAG). The PAG represents a Markov equivalence class of DAGs with latent and selection variables in the acyclic case (Zhang, 2008), and a Markov equivalence class of directed graphs with latent variables (but without selection variables) in the cyclic $\sigma$-separation case. A PAG contains the following types of edges: o-o, o–, o->, –>, <->, —. The bidirected edges come from latent confounders, and the undirected edges come from latent selection variables. The edges have the following interpretation: (i) there is an edge between x and y if and only if variables x and y are conditionally dependent given S for all sets S consisting of all selection variables and a subset of the observed variables (assuming the Markov property and faithfulness); (ii) a tail x --* y at x on an edge between x and y means that x is an ancestor of y or S; (iii) an arrowhead x <-* y at x on an edge between x and y means that x is not an ancestor of y, nor of S; (iv) a circle mark x o-* y at x on an edge between x and y means that there exists both a graph in the Markov equivalence class where x is ancestor of y or S, and one where x is not ancestor of y, nor of S. For further information on the interpretation of PAGs see e.g. (Zhang, 2008) and (Mooij and Claassen, 2020).

The first part of the FCI algorithm is analogous to the PC algorithm. It starts with a complete undirected graph and estimates an initial skeleton using skeleton(*, method="stable") which produces an initial order-independent skeleton, see skeleton for more details. All edges of this skeleton are of the form o-o. Due to the presence of hidden variables, it is no longer sufficient to consider only subsets of the neighborhoods of nodes x and y to decide whether the edge x o-o y should be removed. Therefore, the initial skeleton may contain some superfluous edges. These edges are removed in the next step of the algorithm which requires some orientations. Therefore, the v-structures are determined using the conservative method (see discussion on conservative below).

After the v-structures have been oriented, Possible-D-SEP sets for each node in the graph are computed at once. To decide whether edge x o-o y should be removed, one performs conditional indepedence tests of x and y given all possible subsets of Possible-D-SEP(x) and of Possible-D-SEP(y). The edge is removed if a conditional independence is found. This produces a fully order-independent final skeleton as explained in Colombo and Maathuis (2014). Subsequently, the v-structures are newly determined on the final skeleton (using information in sepset). Finally, as many as possible undetermined edge marks (o) are determined using (a subset of) the 10 orientation rules given by Zhang (2008).

The “Anytime FCI” algorithm was introduced by Spirtes (2001). It can be viewed as a modification of the FCI algorithm that only performs conditional independence tests up to and including order m.max when finding the initial skeleton, using the function skeleton, and the final skeleton, using the function pdsep. Thus, Anytime FCI performs fewer conditional independence tests than FCI. To use the Anytime algorithm, one sets type = "anytime" and needs to specify m.max, the maximum size of the conditioning sets.

The “Adaptive Anytime FCI” algorithm was introduced by Colombo et. al (2012). The first part of the algorithm is identical to the normal FCI described above. But in the second part when the final skeleton is estimated using the function pdsep, the Adaptive Anytime FCI algorithm only performs conditional independence tests up to and including order m.max, where m.max is the maximum size of the conditioning sets that were considered to determine the initial skeleton using the function skeleton. Thus, m.max is chosen adaptively and does not have to be specified by the user.

Conservative versions of FCI, Anytime FCI, and Adaptive Anytime FCI are computed if conservative = TRUE is specified. After the final skeleton is computed, all potential v-structures a-b-c are checked in the following way. We test whether a and c are independent conditioning on any subset of the neighbors of a or any subset of the neighbors of c. When a subset makes a and c conditionally independent, we call it a separating set. If b is in no such separating set or in all such separating sets, no further action is taken and the normal version of the FCI, Anytime FCI, or Adaptive Anytime FCI algorithm is continued. If, however, b is in only some separating sets, the triple a-b-c is marked ‘ambiguous’. If a is independent of c given some S in the skeleton (i.e., the edge a-c dropped out), but a and c remain dependent given all subsets of neighbors of either a or c, we will call all triples a-b-c ‘unambiguous’. This is because in the FCI algorithm, the true separating set might be outside the neighborhood of either a or c. An ambiguous triple is not oriented as a v-structure. Furthermore, no further orientation rule that needs to know whether a-b-c is a v-structure or not is applied. Instead of using the conservative version, which is quite strict towards the v-structures, Colombo and Maathuis (2014) introduced a less strict version for the v-structures called majority rule. This adaptation can be called using maj.rule = TRUE. In this case, the triple a-b-c is marked as ‘ambiguous’ if and only if b is in exactly 50 percent of such separating sets or no separating set was found. If b is in less than 50 percent of the separating sets it is set as a v-structure, and if in more than 50 percent it is set as a non v-structure (for more details see Colombo and Maathuis, 2014). Colombo and Maathuis (2014) showed that with both these modifications, the final skeleton and the decisions about the v-structures of the FCI algorithm are fully order-independent. Note that the order-dependence issues on the 10 orientation rules are still present, see Colombo and Maathuis (2014) for more details.

The FCI-JCI extension of FCI was introduced by Mooij et. al (2020). It is an implementation of the Joint Causal Inference (JCI) framework that reduces causal discovery from several data sets corresponding to different contexts (e.g., observational and different interventional settings, or data measured in different labs or in different countries) to causal discovery from the pooled data (treated as a single 'observational' data set). Two types of variables are distinguished in the JCI framework: system variables (describing aspects of the system in some environment) and context variables (describing aspects of the environment of the system). Different assumptions regarding the causal relations between context variables and system variables can be made. The most common assumption (JCI Assumption 1: 'exogeneity') is that no system variable can affect any context variable. The second, less common, assumption (JCI Assumption 2: 'complete randomized context') is that there is no latent confounding between system and context. The third assumption (JCI Assumption 3: 'generic context model') is in fact a faithfulness assumption on the context distribution. The FCI-JCI extension can be used by specifying the subset of the variables that are designated as context variables with the contextVars argument, and by specifying the combination of JCI assumptions that is used with the jci argument. For the default values of these arguments, the JCI extension is not used. The only difference between the FCI-JCI extension and the standard FCI algorithm is that the background knowledge about the PAG implied by the JCI assumptions is exploited at several points in the FCI-JCI algorithm to enforce adjacencies between context variables (in case of JCI Assumption 3) and to orient certain edges adjacent to context variables (in case of JCI Assumptions 1, 2 or 3). The current JCI framework assumes that no selection bias is present, and therefore the FCI-JCI extension should be called with selectionBias = FALSE. For more details on FCI-JCI, and the general JCI framework, see Mooij et al. (2020).

Value

An object of class fciAlgo (see fciAlgo) containing the estimated graph (in the form of an adjacency matrix with various possible edge marks), the conditioning sets that lead to edge removals (sepset) and several other parameters.

Author(s)

Diego Colombo, Markus Kalisch ([email protected]) and Joris Mooij.

References

D. Colombo and M.H. Maathuis (2014). Order-independent constraint-based causal structure learning. Journal of Machine Learning Research 15 3741-3782.

D. Colombo, M. H. Maathuis, M. Kalisch, T. S. Richardson (2012). Learning high-dimensional directed acyclic graphs with latent and selection variables. Ann. Statist. 40, 294-321.

M.Kalisch, M. Maechler, D. Colombo, M. H. Maathuis, P. Buehlmann (2012). Causal Inference Using Graphical Models with the R Package pcalg. Journal of Statistical Software 47(11) 1–26, doi:10.18637/jss.v047.i11.

J. M. Mooij, S. Magliacane, T. Claassen (2020). Joint Causal Inference from Multiple Contexts. Journal of Machine Learning Research 21(99), 1-108.

J. M. Mooij and T. Claassen (2020). Constraint-Based Causal Discovery using Partial Ancestral Graphs in the presence of Cycles. In Proc. of the 36th Conference on Uncertainty in Artificial Intelligence (UAI-20), 1159-1168.

P. Spirtes (2001). An anytime algorithm for causal inference. In Proc. of the Eighth International Workshop on Artificial Intelligence and Statistics 213-221. Morgan Kaufmann, San Francisco.

P. Spirtes, C. Glymour and R. Scheines (2000). Causation, Prediction, and Search, 2nd edition, MIT Press, Cambridge (MA).

P. Spirtes, C. Meek, T.S. Richardson (1999). In: Computation, Causation and Discovery. An algorithm for causal inference in the presence of latent variables and selection bias. Pages 211-252. MIT Press.

T.S. Richardson and P. Spirtes (2002). Ancestral graph Markov models. Annals of Statistics 30 962-1030.

J. Zhang (2008). On the completeness of orientation rules for causal discovery in the presence of latent confounders and selection bias. Artificial Intelligence 172 1873-1896.

See Also

fciPlus for a more efficient variation of FCI; skeleton for estimating a skeleton using the PC algorithm; pc for estimating a CPDAG using the PC algorithm; pdsep for computing Possible-D-SEP for each node and testing and adapting the graph accordingly; qreach for a fast way of finding Possible-D-SEP for a given node.

gaussCItest, disCItest, binCItest, dsepTest and dsepAMTest as examples for indepTest.

Examples

##################################################
## Example without latent variables
##################################################

set.seed(42)
p <- 7
## generate and draw random DAG :
myDAG <- randomDAG(p, prob = 0.4)

## find skeleton and PAG using the FCI algorithm
suffStat <- list(C = cov2cor(trueCov(myDAG)), n = 10^9)
res <- fci(suffStat, indepTest=gaussCItest,
           alpha = 0.9999, p=p, doPdsep = FALSE)

##################################################
## Example with hidden variables
## Zhang (2008), Fig. 6, p.1882
##################################################

## create the graph g
p <- 4
L <- 1 # '1' is latent
V <- c("Ghost", "Max","Urs","Anna","Eva")
edL <- setNames(vector("list", length=length(V)), V)
edL[[1]] <- list(edges=c(2,4),weights=c(1,1))
edL[[2]] <- list(edges=3,weights=c(1))
edL[[3]] <- list(edges=5,weights=c(1))
edL[[4]] <- list(edges=5,weights=c(1))
g <- new("graphNEL", nodes=V, edgeL=edL, edgemode="directed")

## compute the true covariance matrix of g
cov.mat <- trueCov(g)

## delete rows and columns belonging to latent variable L
true.cov <- cov.mat[-L,-L]
## transform covariance matrix into a correlation matrix
true.corr <- cov2cor(true.cov)

## The same, for the following three examples
indepTest <- gaussCItest
suffStat <- list(C = true.corr, n = 10^9)

## find PAG with FCI algorithm.
## As dependence "oracle", we use the true correlation matrix in
## gaussCItest() with a large "virtual sample size" and a large alpha:

normal.pag <- fci(suffStat, indepTest, alpha = 0.9999, labels = V[-L],
                  verbose=TRUE)

## find PAG with Anytime FCI algorithm with m.max = 1
## This means that only conditioning sets of size 0 and 1 are considered.
## As dependence "oracle", we use the true correlation matrix in the
## function gaussCItest with a large "virtual sample size" and a large
## alpha
anytime.pag <- fci(suffStat, indepTest, alpha = 0.9999, labels = V[-L],
                   type = "anytime", m.max = 1,
                   verbose=TRUE)

## find PAG with Adaptive Anytime FCI algorithm.
## This means that only conditining sets up to size K are considered
## in estimating the final skeleton, where K is the maximal size of a
## conditioning set found while estimating the initial skeleton.
## As dependence "oracle", we use the true correlation matrix in the
## function gaussCItest with a large "virtual sample size" and a large
## alpha
adaptive.pag <- fci(suffStat, indepTest, alpha = 0.9999, labels = V[-L],
                    type = "adaptive",
                    verbose=TRUE)

## define PAG given in Zhang (2008), Fig. 6, p.1882
corr.pag <- rbind(c(0,1,1,0),
                  c(1,0,0,2),
                  c(1,0,0,2),
                  c(0,3,3,0))

## check if estimated and correct PAG are in agreement
all(corr.pag == normal.pag  @ amat) # TRUE
all(corr.pag == anytime.pag @ amat) # FALSE
all(corr.pag == adaptive.pag@ amat) # TRUE
ij <- rbind(cbind(1:4,1:4), c(2,3), c(3,2))
all(corr.pag[ij] == anytime.pag @ amat[ij]) # TRUE



##################################################
## Joint Causal Inference Example
## Mooij et al. (2020), Fig. 43(a), p. 97
##################################################

# Encode MAG as adjacency matrix
p <- 8 # total number of variables
V <- c("Ca","Cb","Cc","X0","X1","X2","X3","X4") # 3 context variables, 5 system variables
# amat[i,j] = 0 iff no edge btw i,j
# amat[i,j] = 1 iff i *-o j
# amat[i,j] = 2 iff i *-> j
# amat[i,j] = 3 iff i *-- j
amat <- rbind(c(0,2,2,2,0,0,0,0),
              c(2,0,2,0,2,0,0,0),
              c(2,2,0,0,2,2,0,0),
              c(3,0,0,0,0,0,2,0),
              c(0,3,3,0,0,3,0,2),
              c(0,0,3,0,2,0,0,0),
              c(0,0,0,3,0,0,0,2),
              c(0,0,0,0,2,0,3,0))
rownames(amat)<-V
colnames(amat)<-V

# Make use of d-separation oracle as "independence test"
indepTest <- dsepAMTest
suffStat<-list(g=amat,verbose=FALSE)

# Derive PAG that represents the Markov equivalence class of the MAG with the FCI algorithm
# (assuming no selection bias)
fci.pag <- fci(suffStat, indepTest, alpha = 0.5, labels = V,
verbose=TRUE, selectionBias=FALSE)

# Derive PAG with FCI-JCI, the Joint Causal Inference extension of FCI
# (assuming no selection bias, and all three JCI assumptions)
fcijci.pag <- fci(suffStat, indepTest, alpha = 0.5, labels = V,
verbose=TRUE, contextVars=c(1,2,3), jci="123", selectionBias=FALSE)

# Report results
cat('True MAG:\n')
print(amat)
cat('PAG output by FCI:\n')
print(fci.pag@amat)
cat('PAG output by FCI-JCI:\n')
print(fcijci.pag@amat)

# Read off causal features from the FCI PAG
cat('Identified absence (-1) and presence (+1) of ancestral causal relations from FCI PAG:\n')
print(pag2anc(fci.pag@amat))
cat('Identified absence (-1) and presence (+1) of direct causal relations from FCI PAG:\n')
print(pag2edge(fci.pag@amat))
cat('Identified absence (-1) and presence (+1) of pairwise latent confounding from FCI PAG:\n')
print(pag2conf(fci.pag@amat))

# Read off causal features from the FCI-JCI PAG
cat('Identified absence (-1) and presence (+1) of ancestral causal relations from FCI-JCI PAG:\n')
print(pag2anc(fcijci.pag@amat))
cat('Identified absence (-1) and presence (+1) of direct causal relations from FCI-JCI PAG:\n')
print(pag2edge(fcijci.pag@amat))
cat('Identified absence (-1) and presence (+1) of pairwise latent confounding from FCI-JCI PAG:\n')
print(pag2conf(fcijci.pag@amat))

Description

This class of objects is returned by functions fci(), rfci(), fciPlus, and dag2pag and represent the estimated PAG (and sometimes properties of the algorithm). Objects of this class have methods for the functions plot, show and summary.

Usage

## S4 method for signature 'fciAlgo'
show(object)
## S3 method for class 'fciAlgo'
print(x, amat = FALSE, zero.print = ".", ...)

## S4 method for signature 'fciAlgo'
summary(object, amat = TRUE, zero.print = ".", ...)
## S4 method for signature 'fciAlgo,ANY'
plot(x, y, main = NULL, ...)

Arguments

Slots

The slots call, n, max.ord, n.edgetests, sepset, and pMax are inherited from class "gAlgo", see there.

In addition, "fciAlgo" has slots

amat:

adjacency matrix; for the coding of the adjacency matrix see amatType

allPdsep

a list: the ith entry of this list contains Possible D-SEP of node number i.

n.edgetestsPDSEP

the number of new conditional independence tests (i.e., tests that were not done in the first part of the algorithm) that were performed while checking subsets of Possible D-SEP.

max.ordPDSEP

an integer: the maximum size of the conditioning sets used in the new conditional independence that were performed when checking subsets of Possible D-SEP.

Extends

Class "gAlgo".

Methods

plot

signature(x = "fciAlgo"): Plot the resulting graph

show

signature(object = "fciAlgo"): Show basic properties of the fitted object

summary

signature(object = "fciAlgo"): Show details of the fitted object

Author(s)

Markus Kalisch and Martin Maechler

See Also

fci, fciPlus, etc (see above); pcAlgo

Examples

## look at slots of the class
showClass("fciAlgo")

## Also look at the extensive examples in ?fci , ?fciPlus, etc !

## Not run: 
## Suppose, fciObj is an object of class fciAlgo
## access slots by using the @ symbol
fciObj@amat   ## adjacency matrix
fciObj@sepset ## separation sets

## use show, summary and plot method
fciObj ## same as  show(fciObj)
show(fciObj)
summary(fciObj)
plot(fciObj)

## End(Not run)

Description

Estimate a Partial Ancestral Graph (PAG) from observational data, using the FCI+ (Fast Causal Inference) algorithm, or from a combination of data from different (e.g., observational and interventional) contexts, using the FCI+-JCI (Joint Causal Inference) extension.

Usage

fciPlus(suffStat, indepTest, alpha, labels, p, verbose=TRUE,
	selectionBias = TRUE, jci = c("0","1","12","123"), contextVars = NULL)

Arguments

Details

A (possibly much faster) variation of FCI (Fast Causal Inference). For details, please see the references, and also fci.

Value

An object of class fciAlgo (see fciAlgo) containing the estimated graph (in the form of an adjacency matrix with various possible edge marks), the conditioning sets that lead to edge removals (sepset) and several other parameters.

Author(s)

Emilija Perkovic, Markus Kalisch ([email protected]) and Joris Mooij.