{'fraction_non_alphanumeric': 0.023446969696969695, 'fraction_numerical': 0.014545454545454545, 'mean_word_length': 4.88126531521497, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Impact of Variable Ordering on Bayesian Network Structure Learning

Causal Bayesian Networks provide an important tool for reasoning under uncertainty with potential application to many complex causal systems. Structure learning algorithms that can tell us something about the causal structure of these systems are becoming increasingly important. In the literature, the validity of these algorithms is often tested for sensitivity over varying sample sizes, hyper-parameters, and occasionally objective functions. In this paper, we show that the order in which the variables are read from data can have much greater impact on the accuracy of the algorithm than these factors. Because the variable ordering is arbitrary, any significant effect it has on learnt graph accuracy is concerning, and this raises questions about the validity of the results produced by algorithms that are sensitive to, but have not been assessed against, different variable orderings.


Introduction
Understanding and modelling causal relationships in complex systems is fundamental to understanding how these systems behave, and in particular how we can intervene to change that behaviour. Such systems can be modelled by Causal Bayesian Networks (CBNs) which are probabilistical graphical models where variables of interest are represented as nodes in a Directed Acyclic Graph (DAG). A directed arc between two nodes is taken to represent a direct causal influence of one variable on the other (Spirtes et al., 2000;Pearl, 2009). The strength and form of the causal relationship is expressed using model parameters that specify conditional distributions.
The key challenge in defining a CBN is to specify the DAG accurately, and here we focus on algorithms that are used to learn the DAG from observational data. There are two main classes of algorithm which can be used to learn graphs from observational data. Constraint-based algorithms, such as the PC algorithm (Spirtes and Glymour, 1991), use statistical tests of marginal and conditional independence to identify dependencies and interdependencies in the data and deduce the possible causal relations from these. This class is often considered the most principled approach to learning causal relationships from data.
The second major class of learning consists of score-based algorithms that search over a set of possible graphs and have an objective score which indicates how well each graph visited fits the data. The objective function is usually either a Bayesian score such as BDeu (Heckerman et al., 1995) which measures how likely the data would be to be generated by that DAG, or an entropy-based score such as BIC (Suzuki, 1999) which balances modelfitting with the dimensionality of the model. These scores are decomposable so that the score for the whole DAG is the sum of the individual scores for each node and its parents, which facilitates computing the scores of neighbouring DAGs. They are also score equivalent which means that all DAGs in an equivalence class have the same score.
The score-based algorithms may be further divided into two subgroups. The first subgroup is referred to as exact algorithms and guarantees to return the highest scoring graph over a set of graphs, an example being GOBNILP (Cussens, 2011). The second subgroup is known as approximate algorithms and often involves heuristic search over the search space that returns a local maximum solution. Examples of approximate algorithms include the hill-climbing (HC) algorithm (Bouckaert, 1994) and variants of greedy search, as well as genetic algorithms (Larranaga et al., 1996). The number of possible DAGs rises superexponentially with the number of variables (Robinson, 1977), and so exact algorithms are typically limited to problems with a few tens of variables. There are also hybrid algorithms which combine constraint-based with score-based solutions. Kitson et al. (2021) provide a recent review of over sixty algorithms covering score-based, constraint-based and hybrid algorithms.
The learnt CBN is only useful if it accurately reflects the true causal relationships in the system under study. There are a number of reasons why the learnt graphical structure might not be an accurate or complete representation of the real world. One group of reasons relates to weaknesses in the algorithms themselves. For example, many approximate scorebased algorithms are greedy and end up returning a graph with only a locally maximum score. Constraint-based algorithms rely on statistical tests of independence which can only say how unlikely the variables are to be independent, and so mistakes are made in the learning process and these mistakes tend to propagate in the algorithms (Spirtes et al., 2000). The accuracy of algorithms is also affected by the objective score or independence tests used (Scutari, 2016), and hyper-parameters which the user must specify such as the independence test threshold p-value, maximum node in-degree, or Equivalent Sample Size (ESS) associated with Bayesian scores (Heckerman et al., 1995).
A second group of causes of inaccuracy relates to the assumptions that many algorithms rely on, and which frequently do not apply in the real world. These include assuming that there are no missing data values, latent confounders or measurement noise, as well as assumptions about the underlying statistical distributions that the variables belong to. Constraint-based algorithms often require another assumption known as faithfulness, which requires that there are no independencies in the data that are not implied by the DAG. Lastly, it is important to note that it is usually not possible to determine the direction of every causal relationship from observational data alone, since the data is generally consistent with a set of DAGs known as the Markov equivalence class (Verma and Pearl, 1990). The equivalence class can be represented by a Completed Partially Directed Acyclic Graph (CPDAG) where undirected edges represent reversible edges that are consistent with the conditional independencies found in the data irrespective of the orientation.
The primary aim of this study is to quantify the effect of variable ordering on the accuracy of the learnt graphs, compared to other factors on which the algorithms are known to be sensitive, such as sample size, objective score and hyper-parameters. We define variable ordering to mean the order that the variables appear as columns in the tabular data structure presented to the structure learning algorithm. The variable ordering is arbitrary and so any significant effect it has on learnt graph accuracy is concerning. Some score-based algorithms, for example K2 (Cooper and Herskovits, 1992), only search across DAGs which are consistent with a predefined ordering, and Hruschka Jr and Ebecken (2007) investigate the use of feature selection to choose that ordering. However, we are not aware of any assessment of the effect of variable ordering across different classes of structure learning algorithms which do not require a variable ordering to be specified.

Methodology
We evaluate the effect of variable ordering on 16 discrete networks ranging from 8 to 109 variables. Many of these networks are widely used in the literature as case studies to evaluate structure learning algorithms. Table 1 lists them and their key characteristics. The Sports, Property, Formed and Diaorrhea networks are available from the Bayesys repository (Constantinou et al., 2020), and the rest obtained from the bnlearn repository (Scutari, 2021b). The network definitions specify the causal structure and the Conditional Probability Tables (CPTs) 1 for the discrete variables. We use these definitions to generate synthetic datasets with sample sizes of {10, 20, 40, 50, 80, 100, 200, . . . ., 8x10 6 , 10 7 }. The default variable ordering of these datasets is alphabetic. Table 1 presents the number of variables and arcs in each network, as well as mean and maximum degree and in-degree. The final column shows what fraction of the arcs are reversible, that is, those edges which are undirected in the CPDAG.
We investigate the effect of variable ordering on two score-based algorithms. The first of these, HC, starts from an empty DAG and compares the scores of each neighbouring DAG, and moves to the neighbouring DAG with the highest score at each iteration. The process continues until there is no neighbouring DAG which improves on the current DAG's score. In the simple variant described here, neighbouring DAGs are defined as ones created by making a single arc addition, deletion or reversal to the current DAG. The final DAG is returned as the learnt DAG. Since HC generally has just a locally maximum score, it represents an approximate algorithm. Variants of the basic HC algorithm have been proposed to try to escape local minima, notably restarts with randomised initial DAGs, and the Tabu list which allows for a limited number of changes that minimally decrease the objective score (Bouckaert, 1995). The Tabu approach is the second score-based algorithm investigated here.
Three constraint-based algorithms are also investigated. The first of these is the PC-Stable algorithm (Colombo et al., 2014), which starts with a complete undirected graph and performs marginal and Conditional Independence (CI) tests to remove edges from the graph. For efficiency, the algorithm starts with marginal independence tests and moves onto conditional independence tests with conditioning sets of increasing size. Conditional The PC-Stable algorithm considers independencies between all variables and constructs the graph "globally", whereas local constraint-based algorithms construct the structure relating to each node separately first, and then merge these to construct the final graph. The local constraint-based GS (Margaritis and Thrun, 1999) and Inter-IAMB (Tsamardinos et al., 2003) algorithms are included in this study.
Two hybrid algorithms are investigated as well. These have an initial local constraintbased algorithm which determines the skeleton of the graph, and then a second HC phase which only considers adding arcs consistent with that skeleton. The first hybrid algorithm MMHC (Tsamardinos et al., 2006) uses a local constraint-based algorithm MMPC (Tsamardinos et al., 2003) in its first constraint-based phase, and the second one, H2PC, uses HPC (Hybrid Parents and Children) (Gasse et al., 2014) in that phase.
The results in this study are obtained using algorithm implementations provided in version 4.7 of the bnlearn package (Scutari, 2009(Scutari, , 2021a. We use the default objective scores, conditional independence test functions and hyper-parameters 2 except for experiments where we explicitly report changing one of those factors. The constraint-based algo-2. For score-based and hybrid algorithms, the default is to use the BIC score with a complexity scaling of 1. For constraint-based algorithms and hybrid algorithms, the Mutual Information Conditional Independence test is used with a p-value threshold of 0.05. (a) Smaller networks (b) Larger networks Figure 1: Arbitrary DAG changes in HC algorithm using 10,000 rows of data rithms investigated produce Partially Directed Acyclic Graphs (PDAGs) 3 , and the scorebased and hybrid algorithms produce a DAG. In order to get consistent evaluations across all algorithm classes, we extend the DAG or PDAG output to a CPDAG where possible, and compare this with the CPDAG to which the reference DAG belongs. The learnt and reference CPDAGs are compared using the F1 metric. The F1 score is widely used in the literature, and is intuitive in that it ranges from 0 to 1 where a higher score represents closer agreement between the learnt and reference graph. Moreover, the F1 score is comparable across networks with different numbers of variables. The F1 metric and CPDAG extension methods provided in bnlearn (Scutari, 2009(Scutari, , 2021a) are used.

Arbitrary DAG Changes in the HC Algorithm
We first examine the individual changes -arc addition, reversal or removal -which the HC algorithm makes at each iteration as it learns the DAG structure. In particular, we note where changes are arbitrary; that is, where two neighbouring DAGs are Markov equivalent. Figure 1 shows the proportion of graphical modifications which are arbitrary for each network with a sample size of 10,000. The first arc added must always be in an arbitrary direction, and hence all lines on the chart start with a proportion of 1.0 arbitrary changes at iteration 1. For the Pathfinder dataset in Figure 1b, the greatest score improvement at iteration 2 happens to be provided by adding an arc onto the first arc to create a chain. The alternative orientation of that edge, which would create a collider, has a smaller score improvement and so this change is not arbitrary. Thus, the proportion of arbitrary changes at iteration 2 for Pathfinder drops to 0.5. The next highest scoring change is an arc addition that does not connect onto this initial chain. Hence that change is arbitrary too and so the proportion of arbitrary changes rises to 2/3 at iteration 3.
We can contrast this behaviour with that of Hailfinder in Figure 1b. In that case, the first eight highest scoring arc additions are all between completely separate pairs of nodes, so that the DAG at iteration eight consists of eight unconnected arcs. Thus, at that iteration, the proportion of arbitrary arcs remains at 1.0. The highest scoring change at iteration 9 happens to be an arc addition which creates a collider with one of the existing arcs and so this non-arbitrary change reduces the proportion to 8/9. The following seven changes are also arbitrary. Figure 1 thus provides an overview of the HC learning process and the proportion of edge modifications whose orientation is determined arbitrarily by the variable order. Based on these results, it is reasonable to conclude that variable ordering has a significant effect in the initial iterations, and that part of this effect might propagate to the final graph. We compare the sensitivity of the F1 score relative to the default alphabetic ordering and two other orderings which we term "optimal" and "worst". In optimal ordering, the variables are ordered so that they are consistent with the topological ordering of the nodes in the reference graph. This optimal ordering ensures that the algorithm always chooses the correct orientation when faced with adding an arc where the two orientations give the same score improvement. Worst ordering is the reverse of optimal, and ensures that the algorithm always makes the wrong decision when adding an arc in the same circumstances. Figure 2 illustrates how the F1 score for each variable ordering changes, as we vary the sample sizes for each of the 16 networks. For many networks, we see that variable ordering makes a large difference in learnt accuracy. In the most extreme cases, such as Asia, Formed, Property, and Hailfinder, there is a difference of more than 0.5 in F1 at large sample sizes. These networks tend to be those which have a high proportion of arbitrary changes as shown in Figure 1. In contrast, there are some networks, particularly Pathfinder, Diarrhoea and Hepar2, where the variable ordering has much less influence on accuracy. These tend to be networks with a smaller proportion of arbitrary changes as the DAG is learnt. The performance of the alphabetic ordering can be viewed as a random ordering, which is why the scores produced for this ordering are mostly between that of the worst and optimal ordering. Figure 2 also provides a comprehensive view into how HC's accuracy varies over a wide range of sample sizes and networks. The broad trend for most networks is that accuracy rises with sample size but then tends to reach a plateau. The height of this plateau varies with the variable ordering and with the network, as does the sample size at which accuracy begins to level off. For some smaller networks such as Asia, Child and Property, maximum accuracy is achieved when sample sizes are of the order of 10 3 to 10 4 rows, whereas larger networks such as Win95pts obtain maximum accuracy at higher sample sizes. Interestingly, some networks such as Water and Pathfinder have not yet plateaued at 10 7 rows. The picture is also complicated because there are some cases where accuracy drops as sample size increases, for example, the Formed and Alarm networks using alphabetic ordering.

The Impact of Variable Ordering on HC Algorithm Graphs
Another issue is the erratic behaviour of some of the plots where there are some large increases or, less often, decreases which occur between two sample sizes. Examples include the Sachs network at around 2,000 rows and Mildew at 10 6 rows. It is suggested that possibly some differences in the learning process have a 'cascading' effect which causes the sequence of changes to be dramatically altered. Further research in this area would be useful. Figure 3 summarises the impact of variable ordering on the accuracy of HC, and compares this with the impact of other factors which affect learning accuracy. Note that for this figure and all subsequent results, we only consider sample sizes in the range 10 3 to 10 6 . This allows us to use the same range of sample sizes across all networks and algorithms 4 . Figure 3 presents a series of 'cat's whiskers' plots showing the distribution of absolute changes to F1 accuracy, where one factor is changed at a time. For example, the leftmost, dark blue plot shows the change in accuracy when the sample size is increased ten fold. To compute this, the F1 from sample size 1,000 is compared with that of sample size 10,000, sample size 2,000 compared with 20,000 and so on, over sample sizes and networks. Similarly, we compare accuracy from sample sizes which are a multiple of 100 times the other in the orange cat's whiskers plot.
The remaining plots show the impact where the sample size is kept the same but where another factor is varied. The green plot shows the change in accuracy resulting from using Figure 3: Impact on F1 scores (CPDAG) of changing sample size, variable ordering, score or hyper-parameters across all networks using the HC algorithm. Each plot shows the mean change as a number, the median change as a horizontal black line, the interquartile range as the coloured rectangle, and the minimum and maximum values as whiskers.
the optimal ordering rather than the alphabetic one, and the red plot shows the change in accuracy moving from the worst to optimal ordering. The last set of plots in Figure 3 shows the impact of some hyper-parameter or objective score changes. These examine the effect of changing the objective score from BIC to BDeu, and changing the objective score from BIC to a more recent Bayesian score, BDS (Scutari, 2016). The last two plots show the impact of changing the complexity scaling factor in the BIC score from 1 to 5 and the Equivalent Sample Size (ESS) in the BDeu score from 1 to 5. Figure 3 illustrates that changing the variable ordering from worst to optimal has the largest overall impact on F1, with a mean improvement of 0.412 and an interquartile range from 0.233 to 0.584. The mean improvement going from alphabetic to optimal ordering is 0.215, with an interquartile range from 0.062 to 0.348. Comparing this with improvements due to other factors, Figure 3 illustrates that changing from alphabetic to optimal ordering results in a higher mean improvement than increasing the sample size by one hundred times (mean improvement 0.161), and considerably higher improvement than increasing sample size by ten times. This observation includes cases where increasing sample size worsens F1 corresponding to the negative gradients seen in Figure 2. Changes to F1 due to changing the objective score from BIC to BDeu and from BIC to BDS have small mean changes of 0.028 and 0.038 respectively. Changing the BIC complexity penalty multiplier hyperparameter and BDeu ESS from 1 to 5 produce a mean change in F1 of -0.080 and -0.004 respectively. Thus, the mean changes and interquartile ranges associated with changing the objective score or hyper-parameters are smaller than the changes associated with changing the sample size, and considerably smaller than the changes associated with changing from optimal to worst or alphabetical variable ordering. Figure 4: The sensitivity of different algorithms to variable ordering compared to other factors in terms of impact on F1 score (CPDAG). The comparisons are: sample size increased by 100 times; variable ordering changed from worst to optimal; score function changed from BIC to BDeu for score-based and hybrid algorithms or the CI test from Mutual Information to Chi-squared for constraint-based algorithms; and, the complexity scaling hyper-parameter is changed from 1 to 5 for score-based and hybrid algorithms or the p-value significance hyper-parameter changed from 0.05 to 0.01 for constraint-based algorithms. Figure 4 shows the sensitivity to variable ordering for the algorithms described in section 2, and compares it with their sensitivity to other selected factors (details are given in the figure caption). TABU is a variant of HC and is, as one might expect, sensitive to variable ordering, but the mean F1 change of 0.278 is considerably less than that of 0.412 for HC. One would expect hybrid algorithms such as H2PC and MMHC, which utilise the HC algorithm in their score-based learning phase to be affected by variable ordering also, and indeed this is the case. However, the constraint-based phase appears to moderate the effect of variable ordering, so that H2PC and MMHC have mean F1 changes of 0.256 and 0.261 respectively.

The Effect of Variable Ordering on other Structure Learning Algorithms
Overall, variable ordering is found to be the most impactful factor on F1 score in three out of the seven algorithms investigated, and the second most impactful factor for the remaining four algorithms. For score-based and hybrid algorithms, variable order has a larger impact than whether BIC or BDeu score is used or setting the BIC complexity penalty to 5. Figure 4 also shows that the constraint-based algorithms are sensitive to variable ordering, with a mean change in F1 of 0.036, 0.104 and 0.021 for PC-Stable, GS and Inter-IAMB respectively. Although these sensitivities are small relative to those observed in other algorithms, they are still larger than the sensitivity of the constraint-based algorithms to the change in CI test and p-value threshold investigated. Note that because constraint-based learning tends to involve dependency tests across sets of variables in triples, it might be possible that other variable orderings have a larger effect on constraint-based algorithms. The sensitivity of PC-Stable is surprising given that it is designed not be sensitive to variable ordering. On the other hand, we note that we find large sensitivities in algorithms outside the PC family, which does not seem to have been considered much in the literature. Table 2 shows the mean F1 scores each of the other algorithms generate relative to the F1 scores of HC, averaged across all networks and sample sizes, for each of the three variable orderings. It thus provides a ranking of the algorithms for each variable ordering. We see that HC is the highest ranked algorithm when using an optimal ordering, since the other six algorithms produce lower F1 scores than HC for that ordering, but that the rankings change for the other two variable orderings. Overall, variable ordering is found to also affect algorithm ranking considerably.

Conclusions
By examining the way a DAG develops iteration by iteration in the simple HC algorithm, we find that arbitrary decisions about edge modifications play an important role in determining the accuracy of the learnt graph and thus, in judging the structure learning capability of an algorithm. This is particularly so in algorithms which 'grow' from several initially unconnected arcs, as this means a higher proportion of the initial arc orientations are arbitrary.
We start by investigating the bnlearn implementation of HC which is widely used in the literature (Scutari, 2009(Scutari, , 2021a and find that these arbitrary edge orientations are made on the basis of variable ordering, which therefore has an impact on the accuracy of the learnt CPDAG. For HC, variable ordering is found to typically have a larger effect than sample size, objective score or hyper-parameters used. We also examine the effect of variable ordering on score-based TABU, on hybrid MMHC and H2PC which use HC in their score-based phase, as well as on constraint-based PC-Stable, Inter-IAMB, and GS. We find that whilst variable ordering has a smaller effect on these algorithms compared to the effect it has on HC, it is still considerable and often stronger than the effect of other factors such as sample size, objective score, and hyper-parameter changes.
Whilst it might be argued that these results are somewhat specific to this particular bnlearn implementation of score-based or constraint-based learning, simple hill-climbing search has to make these choices on some arbitrary basis and so these concerns are likely to be relevant to other implementations, as well as in other algorithms. We find smaller sensitivity for the constraint-based algorithms, though still some, which is surprising in the case of PC-Stable which is designed not to be sensitive to variable ordering.
The effect of variable ordering on structure learning might be well-known for some in the research community, but we argue that its importance is largely underestimated. This is because whilst it is typical for structure learning algorithms to be assessed across different objective functions, varied sample sizes, and different hyper-parameters to convince readers about their validity, almost none of the algorithms published in the literature is tested for sensitivity to variable order. This paper suggests that variable ordering has a much stronger effect on structure learning performance than previously assumed, especially when compared to the other factors on which algorithms are tested for sensitivity, and raises questions about the validity of algorithms which are sensitive to variable ordering and at the same time are oblivious to this effect.