Estimate the final skeleton in the FCI algorithm (Spirtes et al, 2000), as described in Steps 2 and 3 of Algorithm 3.1 in Colombo et al. (2012). The input of this function consists of an initial skeleton that was estimated by the PC algorithm (Step 1 of Algorithm 3.1 in Colombo et al. (2012)).

Given the initial skeleton, all unshielded triples are considered and
oriented as colliders when appropriate. Then, for all nodes x in the
resulting partially directed graph G, Possible-D-SEP(x,G) is computed,
using the function `qreach`

. Finally, for any edge y-z that is
present in G, conditional independence between Y and Z is tested given
all subsets of Possible-D-SEP(y,G) and all subsets of
Possible-D-SEP(z,G). These tests are done at level alpha, using
`indepTest`

. If the pair of nodes is judged to be independent
given some set S, then S is recorded in sepset(y,z) and sepset(z,y)
and the edge y-z is deleted. Otherwise, the edge remains and there is
no change to sepset.

1 2 3 |

`skel` |
Graph object returned by |

`suffStat` |
Sufficient statistic: A list containing all necessary
elements for making conditional independence decisions using
function |

`indepTest` |
Predefined function for testing conditional independence. The
function is internally called as |

`p` |
Number of variables. |

`sepset` |
List of length |

`alpha` |
Significance level for the individual conditional independence tests. |

`pMax` |
Matrix with the maximal p-values of conditional
independence tests in a previous call of |

`m.max` |
Maximum size of the conditioning sets that are considered in the conditional independence tests. |

`pdsep.max` |
Maximum size of Possible-D-SEP for which subsets are
considered as conditioning sets in the conditional independence
tests. If the nodes |

`NAdelete` |
If indepTest returns |

`unfVect` |
Vector containing numbers that encode the unfaithful
triple (as returned by |

`biCC` |
Logical; if |

`verbose` |
Logical indicating that detailed output is to be provided. |

To make the code more efficient, we only perform tests that were not performed in the estimation of the initial skeleton.

Note that the Possible-D-SEP sets are computed once in the beginning. They are not updated after edge deletions, in order to make sure that the output of the algorithm does not depend on the ordering of the variables (see also Colombo and Maathuis (2014)).

A list with the following elements:

`G` |
Updated adjacency matrix representing the final skeleton |

`sepset` |
Updated sepsets |

`pMax` |
Updated matrix containing maximal p-values |

`allPdsep` |
Possible-D-Sep for each node |

`max.ord` |
Maximal order of conditioning sets during independence tests |

`n.edgetests` |
Number of conditional edgetests performed, grouped by the size of the conditioning set. |

Markus Kalisch (kalisch@stat.math.ethz.ch) and Diego Colombo.

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

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

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

`qreach`

to find Possible-D-SEP(x,G);
`fci`

.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 | ```
p <- 10
## generate and draw random DAG:
set.seed(44)
myDAG <- randomDAG(p, prob = 0.2)
## generate 10000 samples of DAG using gaussian distribution
library(RBGL)
n <- 10000
d.mat <- rmvDAG(n, myDAG, errDist = "normal")
## estimate skeleton
indepTest <- gaussCItest
suffStat <- list(C = cor(d.mat), n = n)
alpha <- 0.01
skel <- skeleton(suffStat, indepTest, alpha=alpha, p=p)
## prepare input for pdsep
sepset <- skel@sepset
pMax <- skel@pMax
## call pdsep to find Possible-D-Sep and enhance the skeleton
pdsepRes <- pdsep(skel@graph, suffStat, indepTest, p, sepset, alpha,
pMax, verbose = TRUE)
## call pdsep with biconnected components to find Possible-D-Sep and enhance the skeleton
pdsepResBicc <- pdsep(skel@graph, suffStat, indepTest, p, sepset, alpha,
pMax, biCC= TRUE, verbose = TRUE)
``` |

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

All documentation is copyright its authors; we didn't write any of that.