MRPC: An R package for accurate inference of causal graphs
Abstract
We present MRPC, an R package that learns causal graphs with improved accuracy over existing packages, such as pcalg and bnlearn. Our algorithm builds on the powerful PC algorithm, the canonical algorithm in computer science for learning directed acyclic graphs. The improvement in accuracy results from online control of the false discovery rate (FDR) that reduces false positive edges, a more accurate approach to identifying vstructures (i.e., ), and robust estimation of the correlation matrix among nodes. For genomic data that contain genotypes and gene expression for each sample, MRPC incorporates the principle of Mendelian randomization to orient the edges. Our package can be applied to continuous and discrete data.
1 Introduction: Causal graph inference in R
Graphical models provide a powerful mathematical framework to represent dependence among variables. Directed edges in a graphical model further represent different types of dependencies that may be interpreted as causality [Lauritzen, 1996]. Directed Acyclic Graphs (DAGs), also known as Bayesian networks, are a class of graphical models with only directed edges and no cycles. However, with realworld problems, the true graph contains both directed and undirected edges, or the data may not provide enough information for direction and an edge can be inferred only as undirected. Graphs with both directed and undirected edges are mixed graphs. We treat undirected edges the same as bidirected ones, as not knowing the direction is equivalent to assuming that both directions are equally likely. We term these mixed graphs causal graphs (or equivalently, causal networks) to emphasize the causal interpretation of the directed edges in the graph.
Multiple DAGs may be equivalent in their likelihoods, and therefore belong to the same Markov equivalent class [Richardson, 1997]. Without additional information, inference methods can infer only these Markov equivalent classes. For example, for a simple graph of three nodes, namely , and , if and are conditionally independent given (i.e., ), three Markov equivalent graphs exist:
(1) 
Without additional information, it is not possible to determine which graph is the truth.
Due to Markov equivalence, existing methods typically can infer only a Markov equivalence class of DAGs that can be uniquely described by completed partially directed acyclic graph (CPDAG), which is a mixed graph. Many methods have been developed, which can be broadly classified as i) constraintbased methods, which perform statistical tests of marginal and conditional independence, ii) scoredbased methods, which optimizes the search according to a score function; and iii) hybrid methods that combine the former two approaches [Scutari, 2010]. The PC algorithm (named after its developers Peter Spirtes and Clark Glymour) is one of the first constraintbased algorithms [Spirtes et al., 2000]. Through a series of statistical tests of marginal and conditional independence, this algorithm makes it computationally feasible to infer graphs of high dimensions, and has been implemented in opensource software, such as the R package pcalg [Kalisch et al., 2012]. The PC algorithm consists of two main steps: skeleton identification, which infers an undirected graph from the data;, and edge orientation, which determines the direction of edges in the skeleton. The R package bnlearn implements a collection of graph learning methods from the three classes [Scutari, 2010]. In particular, the function mmhc() implements a hybrid algorithm that also consists of two steps: learning the neighbors (parent and child nodes) of a node, and finding the graph that is consistent with the data and the neighbors identified from the first step [Tsamardinos et al., 2006].
The methods described above are designed for generic scenarios. In genomics there is growing interest in learning causal graphs among genes or other biological entities. Biological constraints, such as the principle of Mendelian randomization (PMR), may be explored for efficient inference of causal graphs in biology. Loosely speaking, the PMR assumes that genotype causes phenotype and not the other way around. Specifically, people in a natural population may carry different genotypes at the same locus on the genome. Such locus is a genetic variant and can be used as an instrumental variable in causal inference. Different genotypes that exist in the population can be thought of as natural perturbations randomly assigned to individuals. Therefore, an association between these genotypes and a phenotype (e.g., expression of a nearby gene) can be interpreted as the genetic variant being causal to the phenotype. This principle has been explored in genetic epidemiology in recent years [Davey Smith and Hemani, 2014]. Under the PMR, five causal graphs involving a genetic variant node and two phenotype nodes are possible (Fig 1; also see Fig 1 in Badsha and Fu [2018]). These five basic causal graphs add constraints to the graph inference process, and allow us to develop the MRPC algorithm that learns a causal graph accurately for hundreds of nodes from genomic data [Badsha and Fu, 2018].
Our algorithm, namely MRPC, is essentially a variant of the PC algorithm that incorporates the PMR. Similar to other PC algorithms, our method also consists of two steps: skeleton identification and edge orientation. We incorporate the PMR in the second step. Meanwhile, our MRPC algorithm is not limited to genomic data. We have developed several improvements over existing methods, and these improvements enable us to obtain more accurate, stable and efficient inference for generic data sets compared to mmhc in the bnlearn package and pc in the pcalg package. In this paper, we explain the details of our MRPC algorithm, describe the main functions in the R package MRPC, demonstrate its functionality with examples, and highlight the key improvements over existing methods. See Badsha and Fu [2018] for comparison with other methods based on the PMR and application to genomic data.
2 Models and software
As PCbased algorithms have demonstrated efficiency in learning causal graphs, we built our algorithm on the pc function implemented in the R package pcalg. Similar to other PCbased algorithm, our MRPC algorithm consists of two steps: learning the graph skeleton, and orienting edges in the skeleton.
The key input of the MRPC function includes the following:

A data matrix with each row being one of individuals, and each column one of features (i.e., a node). Each entry may be continuous or discrete.

The correlation matrix calculated for the features in the data matrix. If the correlation matrix is provided, then the data matrix is ignored. Otherwise, the correlation matrix is calculated from the data matrix. MRPC() allows for calculation of the robust correlation matrix for continuous variables.

A desirable overall false discovery rate (FDR).

The number of genetic variants. This number is 0 if the data do not contain such information.
2.1 Step I: skeleton identification
In order to learn the graph skeleton, which is an undirected graph, we use the procedure implemented in R function pc() in pcalg, but apply the LOND algorithm to control the false discovery rate [Javanmard and Montanari, 2018]. The procedure in pc() is standard in all the PCbased algorithms: it starts with a fully connected graph, and then conducts a series of statistical tests for pairs of nodes, removing the edge if the test for independence is not rejected. The statistical tests consist of those testing marginal independence between two nodes, followed by those of conditional independence tests between two nodes, given one other node, two other nodes, and so on. Each of the tests produces a value, and if the value is larger than the threshold, then the two nodes are independent and the edge between them is removed. Following removal of this edge, the two nodes involved will not be tested again in this step.
The first step involves a sequence of many tests, with the number of tests unknown beforehand. We have implemented the LOND method in our algorithm to control the overall false discovery rate in an online manner [Javanmard and Montanari, 2018]. Specifically, each time we need to conduct a test, we use the LOND method to calculate a value threshold for this test. Depending how many tests have been performed and how many rejections have been made, the value thresholds tend to be large at the beginning and decrease as more tests are performed.
2.2 Step II: edge orientation
We design the following steps for edge orientation:

If the data contain genotype information and there are edges connecting a genotype node to a nongenotype node, then the edge should always point to the nongenotype node, reflecting the assumption in the PMR that genotypes influence phenotypes, but not the other way round.

Next, we search for possible triplets that may form a vstructure (e.g., ). We check conditional test results from step I to see whether and are independent given . If they are, then this is not a vstructure; alternative models for the triplet may be any one of the following three: , , and . If this test is not performed in the first step, we conduct it now.

If there are undirected edges after steps (1) and (2), we look for triplets of nodes with at least one directed edge and no more than one undirected edge. We check the marginal and conditional test results from step I to determine which of the basic models is consistent with the test results (Fig. 1). If we can identify such a basic model, we determine the direction of the undirected edge. We then move on to another candidate triplet. We go through all undirected edges in this step until all of them have been examined. It is plausible that some undirected edges cannot be oriented (see examples below), and we leave them as undirected. Thus, the resulting graph may have both directed and undirected edges, but no directed cycles.
2.3 Hypothesis testing for marginal and conditional independence
A variety of tests may be performed in this step. pcalg includes a commonlyused test for continuous values, which is the (partial) correlation test using Fisher’s z transformation [Kalisch and Bühlmann, 2007]. Briefly, the test statistic between nodes and given , a set of other nodes in the graph, is
(2) 
where is the sample size, the number of nodes in the set , the estimated (partial) correlation between and given . For marginal correlation, is an empty set, and therefore and is the sample correlation between and . This test statistic has a standard normal distribution under the null.
For discrete values, pcalg applies the G test. We use the same function in our MRPC package as well. Briefly, in a contingency table for two discrete variables of levels and of levels, we have the observed count and expected count , where and . The test statistic is
(3) 
which follows a distribution with degrees of freedom. When testing for conditional independence between and given with levels, the test statistic is
(4) 
where and are the observed and expected count for the th level in and the th level in , given the th level in . This statistic has a distribution with degrees of freedom under the null.
2.4 Pearson and robust correlation
The main input to the MRPC() function is the correlation matrix calculated from the data matrix. Pearson correlation is typically used. When outliers are suspected in the data with most nodes having continuous measurements, one may provide the data matrix as the input, and the MRPC() function calculates a robust correlation matrix internally. We implemented the robust correlation calculation method in Badsha et al. [2013]. Specifically, for data that are normal or approximately normal, we calculate iteratively the robust mean vector of length and the robust covariance matrix of until convergence. At the st iteration,
(5) 
and
(6) 
where
(7) 
In the equations above, is the vector of continuous measurements (e.g., gene expression) in the th sample, the sample size, and the tuning parameter. Equation 7 downweighs the outliers through , which takes values in . Larger leads to smaller weights on the outliers. When , Equation 6 becomes a biased estimator of the variance, with the scalar . When the data matrix contains missing values, we perform imputation using the R package mice [Buuren and GroothuisOudshoorn, 2010]. Alternatively, one may impute the data using other appropriate methods, and calculate Pearson or robust correlation matrix as the input for MRPC().
2.5 Metrics for inference accuracy
We provide three metrics for assessing the accuracy of the graph inference.

Recall and precision: Recall (i.e., power, or sensitivity) measures how many edges from the true graph a method can recover, whereas precisions (i.e., 1FDR) measures how many correct edges are recovered in the inferred graph., as follows:
Recall (8) Precision (9) where is the count of edges correctly recovered by the method, the number of edges in the true graph, and the number of edges in the inferred graph.However, we consider it more important to be able to identify the presence of an edge than to also get the direct correct. Therefore, we assign 1.0 to an edge with the correct direction and 0.5 to an edge with the wrong direction or no direction. For example, when the true graph is with 2 true edges, and the inferred graphs are i) , and ; ii) ; and iii) , the number of correctly identified edges is then 2, 1.5 and 1.5, respectively. Recall is calculated to be 2/2=100%, 1.5/2=75%, and 1.5/2=75%, respectively, whereas precision is 2/3=67%, 1.5/2=75%, and 1.5/2=75%, respectively.

Adjusted Structural Hamming Distance (aSHD): The SHD, as is implemented in the R package pcalg and bnlearn, counts how many differences exist between two directed graphs. This distance is 1 if an edge exists in one graph but missing in the other, or if the direction of an edge is different in the two graphs. The larger this distance, the more different the two graphs are. Similar to our approach to recall and precision, we adjusted the SHD to reduce the penalty on the wrong direction of an edge to 0.5.

Number of unique graphs: in our simulation we typically generate multiple data sets under the same true graph. We are also interested in how many unique graphs a method can infer from these data sets. For example, when the true graph is , two inferred graphs may be and . However, the number of wrongly inferred edges is one in both cases, leading to identical recall, precision and aSHD. To distinguish these two graphs, we start with the adjacency matrix for the inferred graph, denoted by , where takes on value 1 if there is a directed edge from node to node , and 0 otherwise. Next we convert this matrix to a binary string, concatenating all the values by row, and treat this binary string as a binary number. We then convert the binary number to decimal, and subtract the decimal number corresponding to the true graph. This decimal difference from the truth serves as a unique identification number for the graph. Additionally, this value being 0 means that the inferred graph is identical to the truth. We count the uniquely inferred graphs in simulation to examine the variation in the graph inference procedure. For example, the true and inferred graphs mentioned above are converted to binary strings and then decimal values as follows:
(10) (11) (12) The difference between the two inferred graphs and the truth is and , respectively, thus representing two distinct graphs with edges wrongly inferred.
2.6 Simulating continuous and discrete data
When simulating data under a true graph, we generated data first for the nodes without parents from a marginal distribution, and then for other nodes from a conditional distribution. When genetic variants are present in the true graph, they are nodes without parents. We assume that each genetic variant node is a biallelic single nucleotide polymorphism (SNP); that is the variant has two alleles (denoted by 0 and 1) in the population. Since an individual receives an allele from their mother and an allele from their father, the genotype at this variant may be 0 for both alleles, 0 for one allele and 1 for the other, or 1 for both alleles. If we count the number of allele 1 in this individual, then the genotype may take on one of three values: 0, 1, and 2. Let be the probability of allele 1 in the population. Assume that the probability of one allele is not affected by that of the other allele in the same individual (i.e., the genotypes are in HardyWeinberg equilibrium). Then the genotype of the node follows a multinomial distribution:
(13) 
Other types of nodes in the graph that are not genetic variants are denoted by . Denote the th nongenotype node by and the set of its parent nodes by , which may be empty, or may include genotype nodes or other nongenotype nodes. We assume that the values of follows a normal distribution
(14) 
The variance may be different for different nodes. For simplicity, we use the same value for all the nodes.
We treat undirected edges as bidirected edges and interpret such an edge as an average of the two directions with equal weights. For example, for graph in Fig. 1, we consider that the undirected edge is a mixture of the two possible directions. Therefore, we generate data for :
(15) 
and separately for :
(16) 
We then randomly choose a pair of values with 50:50 probability for each sample. For simplicity in simulation, we set and all the other s to take the same value, which reflects the strength of the association signal.
The procedure above describes how continuous data are generated. For discrete data, we use the same procedure and then discretize the generated continuous values into categories.
2.7 The LOND method for online FDR control
The LOND algorithm controls FDR in an online manner [Javanmard and Montanari, 2018]. This is particularly useful when we do not know the number of tests beforehand in learning the causal graph. Specifically, consider a sequence of null hypotheses (marginal or conditional independence between two nodes), denoted as , with corresponding pvalue, , . The LOND algorithm aims to determine a sequence of significance level , such that the decision for the th test is
(17) 
The number of rejections over the first tests is then
(18) 
For the overall FDR to be , we use a series of nonnegative values , and set the significance level as
(19) 
where the FDR for the th test is
(20) 
with
(21) 
for integer and constant . The default value for is set to 2 in MRPC. At an FDR of 0.05 and , we have
(22) 
Then
(23) 
The larger is, the more conservative the LOND method, which means that fewer rejections will be made. We therefore set throughout simulation and real data analyses.
3 Main functions of the MRPC package
3.1 Causal graph inference

MRPC(): this function is central to our package and implements the causal graph inference algorithm described in the previous section. The required arguments are data,
suffStat
,GV
,FDR
,indepTest
, andlabels
.data
is the data matrix that may contain discrete or continuous values, and may or may not be genomic data. suffStat is a list of two elements: the (robust) correlation matrix calculated from the data matrix, and the sample size. GV is the number of genetic variants, and should be set to 0 if the data is not genomic or does not have genotype information. FDR is the desirable false discovery rate, a value between 0 and 1. indepTest specifies the statistical test for marginal and conditional independence; typically, this argument is set to gaussCItest for continuous values, or todisCItest
for discrete data. Lastly, labels is the names of the nodes used for visualization. 
ModiSkeleton(): this is a modified version of the skeleton function in the pcalg package, and is called within the MRPC function. Similar to skeleton(), this function also learns the skeleton of a causal graph. However, we have incorporated the LOND algorithm into this function such that the overall FDR is controlled.

EdgeOrientation(): this function is called by
MRPC()
and performs Step II of the MRPC algorithm. 
RobustCor(): this function calculates the robust correlation matrix for a data matrix that does not contain genotype information and have continuous values in all columns (nodes), or that contains genotype information and have continuous values in the phenotype columns. This function may be used when the above data matrix may contain outliers. It returns a correlation matrix that may be used as an argument for MRPC().
3.2 Simulation

SimulatedData(): this function generates genomic data that contains genetic variants.
3.3 Visualization

plot(): this function inherits the functionality of the plot function in the pcalg package. Therefore, if the graph does not display correctly, the user may need to run library (pcalg) to ensure that the pcalg package is loaded into R.

DendroModuleGraph(): this function is modified from the function with the name in the R package WGCNA [Langfelder and Horvath, 2008]. Similar to the WGCNA function, this function also generates a dendrogram of all the nodes in the inferred causal graph based on the distance (i.e., the number of edges) between two nodes, cluster the nodes into modules according to the dendrogram and the minimum module size (i.e., the number of nodes in a module), and draws the inferred causal graph, where nodes of different modules are colored differently. For genomic data, this function draws the genotype nodes in filled triangles and phenotype nodes in filled circles.
3.4 Assessment of inferred graphs

Recall_Precision(): this function compares an inferred graph to a true one, counts the number of true positives and false positives in the inferred graph, and calculates recall and precision.

aSHD(): this function calculates the (adjusted) SHD between two graphs.
4 Illustrations
4.1 Continuous data with and without genetic information
The first example is a genomic data set with two genotype nodes, and , and 4 phenotype nodes, through (Fig. 2a). The genotype nodes are discrete, whereas the phenotype nodes are continuous.
R> data ("ExampleMRPC") R> n < nrow (ExampleMRPC$simple$cont$withGV$data) R> V < colnames(ExampleMRPC$simple$cont$withGV$data) R> Rcor_R < RobustCor(ExampleMRPC$simple$cont$withGV$data, Beta=0.005) R> suffStat_R < list(C = Rcor_R$RR, n = n) R> data.mrpc.cont.withGV < MRPC(data=ExampleMRPC$simple$cont$withGV$data, +ΨsuffStat=suffStat_R, GV=2, FDR=0.05, indepTest=’gaussCItest’, +Ψlabels=V, verbose = TRUE) R> par (mfrow=c(1,2)) R> plot (ExampleMRPC$simple$cont$withGV$graph, main="truth") R> plot (data.mrpc.cont.withGV$graph, main="inferred")
The second example is the data set from the pcalg package, in which none of the nodes are genotype nodes, and all the nodes take on continuous values (Fig. 2b). The command lines for generating the output are similar to those above, with ExampleMRPC$simple$cont$withGV$data being replaced by ExampleMRPC$simple$cont$withoutGV$data, and GV is set to 0 when calling MRPC().
4.2 Discrete data with and without genetic information
The first example here is a genomic data set with one genotype node and five phenotype nodes (Fig. 3a). All the nodes have discrete values.
R> data ("ExampleMRPC") R> n < nrow (ExampleMRPC$simple$disc$withGV$data) R> V < colnames(ExampleMRPC$simple$disc$withGV$data) R> Rcor_R < RobustCor(ExampleMRPC$simple$disc$withGV$data, Beta=0.005) R> suffStat_R < list(C = Rcor_R$RR, n = n) R> data.mrpc.disc.withGV < MRPC(data=ExampleMRPC$simple$disc$withGV$data, +ΨsuffStat=suffStat_R, GV=1, FDR=0.05, indepTest=’gaussCItest’, +Ψlabels=V, verbose = TRUE) R> par (mfrow=c(1,2)) R> plot (ExampleMRPC$simple$disc$withGV$graph, main="truth") R> plot (data.mrpc.disc.withGV$graph, main="inferred")
The second example is a generic data set with five nodes of discrete values, where no nodes are genetic variants (Fig. 3b). The command lines for
generating the output are similar to those above, with ExampleMRPC$simple$disc$withGV$data being replaced by
ExampleMRPC$simple$disc$withoutGV$data, and with GV set to 0 when calling MRPC().
4.3 Simulation and performance assessment
We include a small example to demonstrate how simulation is performed and how we use functions implemented in MRPC to compare the inferred graph to the truth. The user may run the following command line to see the demo:
R> simu < SimulationDemo(N=1000, model="truth1", signal=1.0, n_data=2, +Ψn_nodeordering=6) R> apply (simu, 2, unique)
In this demo, “truth1” represents the graph (also see Figure 6a), which gives rise to six permutations of the three nodes (we do not permute all the nodes, as the nodes, which correspond to the genetic variants, need to be at the beginning of the node list): , , , , , . For each node list, the demo simulates two independent data sets with 1000 observations, applies MRPC, mmhc, and pc to each data set, and calculates the difference between the inferred graph and the truth, after converting each graph into a binary string and then a decimal number. The output is a matrix of 2 rows and columns. The second line above lists the unique numbers (representing uniquely inferred graphs) for each permutation inferred by each method.
4.4 Visualizing a complex graph
For complex graphs with many nodes, it may help interpretation to cluster the nodes into modules. Below, we demonstrate this functionality on a complex graph of 22 nodes. This graph may be a true graph, or an inferred one from running the MRPC function. Figure 4 shows the graph without clustering, the dendrogram of the nodes, and the graph with clustering.
R> plot(ExampleMRPC$complex$cont$withGV$graph) R> Adj_directed < as(ExampleMRPC$complex$cont$withGV$graph,"matrix") R> DendroModuleGraph(Adj_directed, minModuleSize = 5, GV=14)
5 Key improvements over existing methods
5.1 Online FDR control
Existing methods for graph inference, such as methods implemented in packages pcalg and bnlearn, use only the type I error rate for individual tests and do not control for multiple testing. A key improvement in MRPC over these methods is the incorporation of the LOND method for controlling the false discovery rate in an online manner. As the marginal and conditional independence tests are performed sequentially, the LOND method calculates a desirable type I error rate for each test, using the number of tests that have been performed and the number of rejections so far, with the aim of controlling the overall FDR at a predetermined level.
5.2 Vstructure identification
For any triplet with two edges (i.e., ), only two relationships may exist: conditional independence , and conditional dependence (we add the parentheses for clarity; these parentheses are typically not used in standard notation). Three equivalent graphs exist for conditional independence:
(24) 
As mentioned in the Introduction, these three graphs are Markov equivalent and cannot be distinguished. Only one graph is possible for conditional dependence:
(25) 
which is the vstructure. As the only unambiguous graph one may infer for the same skeleton, it is therefore critical to correctly identify all the vstructures in the data.
However, both pc and mmhc may wrongly identify vstructures when there is not one (Figure 5). With pc, the false positive is due to incorrect interpretation of the lack of the edge . Specifically, when testing for marginal independence in the first step of skeleton identification, the null hypothesis of and being independent is likely not rejected, leading to the removal of the edge . Once an edge is removed in the inference, the two nodes are not going to be tested again. This means that there is no evidence for or against . However, pc interprets this lack of evidence as support for and claims a vstructure. As a result, pc usually identifies vstructures when they are present, but also falsely calls vstructures when there are none. It is unclear why mmhc also falsely identifies vstructures. With MRPC, once we identify a candidate triplet of this skeleton, we examine whether the conditional independence test between and given has been performed before, and if not, conduct this test to determine whether the data supports a vstructure.
5.3 Nodeordering independence
It is desirable to learn the same graph even after the nodes in the data set are permuted. Our MRPC algorithm demonstrates stability across different simulations, better than pc, which implements a nodeordering independent algorithm, and much better than mmhc. For example, for each of the two graphs in Figure 6, we simulated a single data set and permuted the nodes (columns). Our method processes the genotype nodes () and the phenotype nodes () separately, and therefore we permuted only the s. Three s have three different permutations. We then ran MRPC, mmhc and pc on each data set to infer the causal graph. Whereas MRPC and pc produced the same graph on three permutations, mmhc generated a different graph for each permutation in Figure 6a and two graphs for three permutations in Figure 6b.
We next generated 200 data sets from each of the graphs in Figure 6. For each data set, we permuted the nodes to derive six data sets. We then applied MRPC, mmhc and pc to all the data sets, counted the number of unique graphs among the six inferred graphs for each data set, and summarized the interquartile range, median and max across all 200 data sets for each true graph (Table 1). mmhc appears to be more unstable compared to the other two methods. Note that we focus on the variation in inferred graphs, not the accuracy. Although pc is stable, it may not infer the graph correctly (Figure 6 (a)).
Graph (a) in Fig 6  

1st Quartile  Median  3rd Quartile  Max  
MRPC  1.00  1.00  1.00  2.00 
mmhc  2.00  3.00  3.00  4.00 
pc  1.00  1.00  1.00  2.00 
Graph (b) in Fig 6  
1st Quartile  Median  3rd Quartile  Max  
MRPC  1.00  1.00  1.00  2.00 
mmhc  2.00  2.00  2.00  3.00 
pc  1.00  1.00  1.00  1.00 
5.4 Robustness in the presence of outliers
When the data may contain outliers, our MRPC package allows for calculation of a robust correlation matrix from the data matrix. Each sample in the data matrix is assigned a weight in the correlation calculation. Outliers tend to receive a weight near 0, thus their contribution to correlation is reduced. Neither pc or mmhc deals with outliers. As our simulation results show, MRPC infers the same graph with or without the presence of outliers, and the inferred graph is close to the truth. By contrast, pc and mmhc wrongly place edges where they do not belong (Figure 7).
Our current implementation of the robust correlation calculation is limited to continuous data for all the columns if there is no genotype information, and for the phenotype columns if there is genotype information.
6 Summary and discussion
Here, we introduce our R package MRPC that implements our novel algorithm for causal graph inference. Our algorithm builds on existing PC algorithms and incorporates the principle of Mendelian randomization for genomic data when genotype and molecular phenotype data are both available at the individual level. Our algorithm also controls the overall FDR for the inferred graph, improves the vstructure identification, reduces dependence on the node ordering, and deals with outliers in the data. These improvements are not limited to genomic data. Therefore, our MRPC algorithm is a much improved algorithm for causal graph learning for generic data. Our MRPC package contains the main function MRPC() for causal graph learning, as well as functions for simulating genomic and nongenomic data from a wide range of graphs, for visualizing the graphs, and for calculating several metrics for performance assessment.
Through simulation, we demonstrated that our method is stable, accurate and efficient on relatively small graphs. However, due to the online FDR control implemented in our package, the inference can be slow when the number of nodes reaches thousands. We have yet to develop more efficient search algorithms while retaining accuracy and stability.
Computational details
The results in this paper were obtained using R 3.4.1 with the MASS 7.3.47 package. The MRPC package is available at https://github.com/audreyqyfu/mrpc. R itself and all packages used are available from the Comprehensive R Archive Network (CRAN) at https://CRAN.Rproject.org/.
Acknowledgments
This research is supported by NIH/NHGRI R00HG007368 (to A.Q.F.) and partially by the NIH/NIGMS grant P20GM104420 to the Center for Modeling Complex Interactions at the University of Idaho.
References
 Lauritzen [1996] Steffen L Lauritzen. Graphical Models, volume 17 of Oxford Statistical Science Series. Clarendon Press, 1996.
 Richardson [1997] Thomas Richardson. A characterization of Markov equivalence for directed cyclic graphs. International Journal of Approximate Reasoning, 17(23):107–162, 1997.
 Scutari [2010] Marco Scutari. Learning bayesian networks with the bnlearn R package. Journal of Statistical Software, 35, 2010. doi: 10.18637/jss.v035.i03.
 Spirtes et al. [2000] P. Spirtes, C. N. Glymour, and R. Scheines. Causation, Prediction, and Search. MIT Press, 2000.
 Kalisch et al. [2012] Markus Kalisch, Martin Mächler, Diego Colombo, Marloes H Maathuis, and Peter Bühlmann. Causal inference using graphical models with the R package pcalg. Journal of Statistical Software, 47(11):1–26, 2012. doi: 10.18637/jss.v047.i11.
 Tsamardinos et al. [2006] I. Tsamardinos, L.E. Brown, and C.F. Aliferis. The maxmin hillclimbing Bayesian network structure learning algorithm. Machine Learning, 65:31–78, 2006. doi: 10.1007/s1099400668897.
 Davey Smith and Hemani [2014] George Davey Smith and Gibran Hemani. Mendelian randomization: genetic anchors for causal inference in epidemiological studies. Human Molecular Genetics, 23(R1):R89–R98, 2014.
 Badsha and Fu [2018] Md Bahadur Badsha and Audrey Qiuyan Fu. Learning causal biological networks with the principle of mendelian randomization. bioRxiv, 2018. doi: 10.1101/171348.
 Javanmard and Montanari [2018] Adel Javanmard and Andrea Montanari. Online rules for control of false discovery rate and false discovery exceedance. The Annals of statistics, 46(2):526–554, 2018. doi: 10.1214/17AOS1559.
 Kalisch and Bühlmann [2007] Markus Kalisch and Peter Bühlmann. Estimating highdimensional directed acyclic graphs with the PCalgorithm. Journal of Machine Learning Research, 8:613–636, 2007.
 Badsha et al. [2013] M. B. Badsha, M. N. Mollah, N. Jahan, and H. Kurata. Robust complementary hierarchical clustering for gene expression data analysis by betadivergence. Journal of Bioscience and Bioengineering, 2013. doi: 10.1016/j.jbiosc.2013.03.010.
 Buuren and GroothuisOudshoorn [2010] S van Buuren and Karin GroothuisOudshoorn. mice: Multivariate imputation by chained equations in R. Journal of statistical software, pages 1–68, 2010. doi: 10.18637/jss.v045.i03.
 Langfelder and Horvath [2008] Peter Langfelder and Steve Horvath. WGCNA: an R package for weighted correlation network analysis. BMC bioinformatics, 9(1):559, 2008. doi: 10.1186/147121059559.