Scoring Bayesian Networks of Mixed Variables

2.1 The Bayesian Information Criterion

The Bayesian Information Criterion (BIC) is a well studied and widely applied marginal log likelihood approximation that is generally used for model selection. Let M be a model we wish to score given a dataset $\mathcal{D}$. We can write the probability of model M given $\mathcal{D}$ using Bayes’ rule as,

However, since the data are fixed, $p(\mathcal{D})$ will remain constant across different model choices. Thus, for model selection, we use,

BIC aims to approximate $p(M|\mathcal{D})$ in (1). For now, we assume p(M) is distributed uniformly and thus drop it. Later in section 5.1, we introduce an alternative distribution for p(M), which we find performs well in practice. Raftery [17] shows that when assuming a flat prior over the parameters, the logarithm of $p(\mathcal{D}|M)$ can be approximated as:

where ℓ(θ) is the maximum log likelihood of the data, df are the degrees of freedom, and n is the sample size. The approximation on the right hand side of (2) characterizes the BIC, introduced by Schwarz [23]. BIC is decomposable and can be readily applied to score Bayesian networks. In sections 3 and 4, we detail how to calculate the log likelihood and degrees of freedom terms for BIC using our proposed scoring methods. We score DAGs using BIC given such log likelihood and degrees of freedom calculations.

In general, the Conditional Gaussian (CG) score calculates conditional Gaussian mixtures using the ratios of joint distributions. Since CG uses BIC as a framework to evaluate its approximations, the score is decomposable into a sum of parent-child relationships. In order to outline such a relationship, we introduce continuous variables C₁, C₂ and discrete variables D₁, D₂ Below we detail the CG score using these four variables, however this procedure straightforwardly generalizes to any number of variables. We compute the conditional distribution where C₁ is a child with parents C₂, D₁, and D₂ as,

and the conditional distribution where D₁ is a child with parents C₁, C₂, and D₂ as,

In (3) and (4), we can readily calculate p(C₁, C₂|D₁, D₂) and p(C₁, C₂|D₂) using Gaussian distributions partitioned on the discrete variables and p(D₁, D₂), p(D₂) using multinomial distributions. This raises the first of CG’s assumptions.

3.1 Modeling a Set of Variables

When using CG, we have three different kinds of sets to model: Y_i ⋃ Pa_i where Y_i is continuous, Y_i ⋃ Pa_i where Y_i is discrete, and Pa_i. They all follow the same generic format so we will describe the process in general while pointing out any subtle differences where they apply.

First we partition the data with respect to a partitioning set ∏_i generated according to the discrete variables Pd_i. Note that if our set includes a discrete child Y_i, then the discrete variables are comprised of Y_i ⋃ Pd_i and we partition according to these variables. ∏_i contains a partition for every combination of values in the discrete variables. We define the partitioning set ∏_i using a Cartesian product of the discrete variables. Let |Pd_i| = d, then partitioning set ∏_i = (Y_i)×Pd_i(1)×Pd_i(2)×⋯ ×Pd_i(d) where Y_i is the set of values for the child (included only if Y_i is discrete), Pd_i(1) is the set of values for the first discrete parent, Pd_i(2) is the set of values for the second discrete parent, and so forth.

Let |Pc_i| = c, then for each partition p ∈ ∏_i we define a design matrix X_p with n_p observations and c variables. Here, if our set includes a continuous child Y_i, then we instead define X_p with c + 1 variables corresponding to the variables in Y_i ⋃ Pc_i. That is,

where x_jk is the j^th value with respect to partition p of the k^th variable in Pc_i and y_j is the j^th value with respect to p of the child Y_i (included only if Y_i is continuous) for j ∈ {1, 2,…, n_p} and k ∈ {1, 2,…, c}.

3.2 Calculating the Log Likelihood and Degrees of Freedom

The calculations for the three aforementioned sets are identical in formulation, so without loss of generality, we demonstrate the log likelihood and degrees of freedom calculations for the set Y_i ⋃ Pa_i. The log likelihood for a set is calculated component-wise over each partition and summed together as follows,

The degrees of freedom are calculated in a similar manner,

where the minus 1 term accounts for the redundant mixing component. Let n be the number of observations in the unpartitioned dataset. For each partition p ∈ ∏_i, let d be the number of variables in X_p and x_p,j be the j^th observation from X_p. From [2], we calculate the Gaussian log likelihood for partition p as,

where μ_p, Σ_p are the mean and variance of the Gaussian distribution respectively. The maximum likelihood estimate ${\widehat{\sum }}_p$ is computed as,

Let ${\mathit{\mu }}_p={\overline{x}}_p$. Note that ${\overline{\mathit{x}}}_p$ will converge quickly to μ_p. Using the estimate in (8), the log likelihood in (7) simplifies to,

We use (9) to compute the log likelihood of a Gaussian conditioned on discrete variables. However, we still must add the log probability of an instance being from partition p to calculate the desired joint log likelihood. These probabilities are computed using the maximum likelihood estimate of variables distributed according to a multinomial. This estimate is the count of instances in partition p denoted n_p over the total count n of all instances: $\frac{n_p}{n}$. Thus, we calculate the log likelihood for partition p as,

We use (10) to calculate ℓ_p(θ_p|X_p) in (5). To find the number of parameters in partition p, we count the number of unique terms in ${\widehat{\sum }}_p$ plus one for the mixing component. Therefore,

We use (11) to calculate df_p(θ_p) in (6).

Using the form of (3) and (4), we calculate the log likelihood and degrees of freedom terms as,

BIC uses (12) and (13) to compute the score for the parent-child relationship of Y_i given Pa_i.

The Mixed Variable Polynomial (MVP) score uses higher order polynomial functions to approximate relationships between any number of continuous and discrete variables. Since MVP uses BIC as a framework to evalute its approximations, the score is decomposable into a sum of parent-child relationships. The MVP method scores the decomposed local components of a DAG $\mathcal{G}$ using approximating polynomial functions. To motivate the ideas underlying this approach, we note the implications of Weierstrass’s approximation theorem.

4.1 Partitioned Regression

Let Y_i be the i^th variable in a DAG $\mathcal{G}$ and Pa_i be the set containing the parents of Y_i in $\mathcal{G}$. Furthermore, let Pa_i consist of two mutually exclusive subsets Pc_i and Pd_i such that Pc_i and Pd_i hold the continuous and discrete parents of Y_i respectively. In general, to evaluate the local score component between Y_i and its parents Pa_i, MVP first partitions the data with respect to the discrete parents Pd_i and performs least squares regression using the continuous parents Pc_i. The log likelihood and degrees of freedom for the model are calculated depending on the variable type of Y_i. BIC uses the log likelihood and degrees of freedom terms to compute the score.

A partitioning set ∏_i partitions $\mathcal{D}$ with respect to the discrete parents Pd_i and contains a partition for every combination of values in the discrete parents. We define ∏_i using a Cartesian product of the discrete parents Pd_i. Let |Pd_i| = d, then partitioning set ∏_i = Pd_i (1) × Pd_i (2) × ⋯ × Pd_i (d) where Pd_i (1) is the set of values for the first discrete parent, Pd_i (2) is the set of values for the second discrete parent, and so forth.

Let |Pc_i| = c, then for each partition p ∈ ∏_i we define a design matrix X_p with n_p observations and c variables. Additionally, we add a bias term and higher order polynomial terms for each variable in Pc_i, stopping at a maximum polynomial order specified by g(n_p),

where x_jk is the j^th value with respect to partition p of the k^th variable in Pc_i for j ∈ {1, 2, …,n_p} and k ∈ {1, 2,…, c}. In this paper, we report two choices for g(n_p): g(n_p) = 1, and g(n_p) = ⌊log n_p⌋. We have tried other choices, such as g(n_p) = 3, but found the above options provide the best solutions. Define x_j and y_p as,

where x_p,j is the j^th observation in X_p and y_j is the j^th value with respect to partition p of Y_i for j ∈ {1, 2,…, n_p}.

We calculate the log likelihood of a variable Y_i given a set of parents Pa_i as a sum over the log likelihoods from each partition p,

where ℓ_p(θ_p |X_p, y_p) is defined depending on whether Y_i is discrete or continuous. Similarly, the degrees of freedom for Y_i are calculated as a sum over the parameter counts in each partition p,

BIC computes the local score component using the log likelihood ℓ_i(θ|X, y) and degrees of freedom df_i (θ).

4.2 Modeling a Continuous Child

In the case where Y_i is continuous, for partition p with design matrix X_p and target vector y_p, we use the maxmum likelihood estimate to determine the log likelihood ℓ_p(θ_p|X_p, y_p) and degrees of freedom df_p(θ_p). By Assumption 2, for partition p, each y_j = f_p(x_p,j) + ε_p where f_p is a continuous function defined on a bounded interval and ${\epsilon }_p~N(0,{\sigma }_p^2)$ is additive Gaussian noise with variance ${\sigma }_p^2$. By the Weierstrass Approximation Theorem, there exists a polynomial function ${\widehat{f}}_p$ which approximates f_p such that y_j ≈ f_p(x_p,j) + ε_p. Estimating the parameters of $y_j\approx {\widehat{f}}_p({\mathit{x}}_{p,j})+{\epsilon }_p$ using least squares regression, we have ${\widehat{y}}_j~N({\mathit{X}}_p{\mathit{\beta }}_p,{\sigma }_p^2)$. Therefore the log likelihood partition p becomes,

where the maximum likelihood estimates are computed as,

Using the estimates in (17) and (18), the log likelihood in (16) simplifies to,

which we use to calculate ℓ_p(θ_p|X_p, y_p) in (14). To find the number of parameters in partition p to calculate df_p(θ_p) for (15), we count the number of terms in ${\widehat{\mathit{\beta }}}_p$,

The BIC uses (14) and (15) to compute the parent-child relationship for Y_i given Pa_i.

4.3 Modeling a Discrete Child

In the case where Y_i is discrete, for partition p with design matrix X_p and target vector y_p, we calculate the log likelihood ℓ_p(θ_p|X_p, y_p) and degrees of freedom df_p(θ_p|X_p, y_p) using least squares regression. Suppose Y_i consists of d categories. Let f_p,_h calculate the probability of the h^th category in Y_i given the values of the continuous parents Pc_i where h ∈ {1,…,d}. By Assumption 2 and the Weierstrass Approximation Theorem, there exists a polynomial function ${\widehat{f}}_{p,h}$ which approximates f_p,_h arbitrarily well. With this in mind, we aim to approximate each f_p,_h with the polynomial function ${\mathit{X}}_p{\widehat{\mathit{\beta }}}_{p,h}$ where ${\widehat{\mathit{\beta }}}_{p,h}$ are the polynomial coefficients calculated from least squares regression. Our end goal is to use the approximations of each f_p,_h as components of a conditional probability mass function in order to calculate the log likelihood and degrees of freedom terms.

Define the categories of y_p such that each yj ∈ {1,…, d}. Expand y_p into d binary vectors where the h^th vector represents the h^th category. That is, the j^th element from the h^th vector asserts whether or not the j^th observation of y_p is category h. To represent these binary vectors in matrix notation, we define 1{_condition} as an indicator variable which is 1 if condition is true and 0 otherwise. The h^th binary vector is then defined as,

We further define 1_p to represent an n_p × 1 vector of ones. Using the binary vectors as our targets, we calculate the least squares regression estimate, which yields,

If we want to interpret the results of least squares regression ${\mathit{X}}_p{\widehat{\mathit{\beta }}}_{p,h}$ as probabilities, we first must ensure that,

1. ${\displaystyle {\sum }_{h=1}^d{\mathit{X}}_p{\widehat{\mathit{\beta }}}_{p,h}={\mathbf{1}}_p}$

2. ${\mathit{x}}_{p,j}{\widehat{\mathit{\beta }}}_{p,h}\ge 0$, ∀ j ∈ {1, …,n_p}, h ∈ {1, …,d}.

We prove condition 1 is necessarily true and that condition 2 holds for an adequate maximum polynomial degree in the sample limit; see Appendix A for the proofs.

Unfortunately, there is no guarantee the values given by least squares regression will be strictly non-negative for finite samples. Instead we define a procedure that maps the values calculated for each y_j to a positive value. In the procedure, we avoid setting values directly to zero in order to prevent assigning zero probability to any observed instances. Instead, we use a value which tends towards zero; we insure the minimum value of each y_j·is at least $\frac{1}{n_p}$ and maintain the condition that each set of estimates sums to one. Therefore, we can treat the mapped estimates as probabilities. Algorithm 1 outlines how the procedure accomplishes this mapping. Figure 1 shows the effect of Algorithm 1 when applied to the least squares estimates of the conditional probability mass function for a particular category of a discrete variable.

In words, our procedure is as follows:

1. Shift the estimates such that they are centered about a non-informative center by subtracting $\frac{1}{d}$ (line 6).

2. Scale the estimates such that the smallest final values will be at least $\frac{1}{n_p}$ (line 7).

3. Shift the scaled estimates back to the original center by adding $\frac{1}{d}$ (line 8).

Algorithm 1.

Maps least squares estimates to valid probability distributions.

Open in a new tab

Since we only want to perform this procedure if one of the least squares estimates is negative, we define $m_{p,j}=min\{\frac{1}{n_p},{\mathit{x}}_{p,j}{\widehat{\mathit{\beta }}}_{p,h}\forall h\}$ so that m_p,j is either the minimum estimate for y_j or $\frac{1}{n_p}$ (line 4). We calculate the scaling factor α_p,j (line 5) by noting that we want,

Solving for α_p,j we find,

Note, that if $m_{p,j}=\frac{1}{n_p}$, then α_p,j = 1 and we do not transform the estimates. We compute the log likelihood in the discrete case as

We use (21) to calculate ℓ(θ_p|X_p, y_p) in (14). To find the number of parameters in p, we count the number of terms across all ${\widehat{\mathit{\beta }}}_{p,j}$. Each ${\widehat{\mathit{\beta }}}_{p,j}$ has c · g(n_p) + 1 parameters and j ranges over d categories. However, since Proposition 1 shows the estimated probabilities sum to one, the number of free parameters is

As before, BIC uses (14) and (15) to compute the parent-child relationship for Y_i given Pa_i.

5.1 Binomial Structure Prior

We introduce a structure prior inspired by the binomial distribution. The idea is to give a prior distribution over the number of parents of each variable. We view the addition of each edge as an independent event that occurs with probability q. Therefore, we expect to see q • m′ parents for any given variable where m′ is the total number of possible parents. Then we have,

where π(k) is the prior probability that any given variable Y_i in DAG $\mathcal{G}$ has k parents. Often it is more convenient to work in log space. Thus we calculate the log prior probability as,

Note that since DAGs encode an ordering over the variables, the total number of possible parents is not necessarily all the variables in the data excluding the variable currently acting as a child. It is usually the case that m′ ≠ m − 1 where m is the total number of variables in the data.

In section 6, we let m′ = m in order to calculate the binomial structure prior more efficiency.

We calculate q as $q=\frac{r}{(m-1)}$, where r represents a user-specified upper bound on the expected number of parents of any given node.

Usually, BIC assumes the prior probability of models in equation (1) is distributed uniformly. By using the binomial structure prior instead, we adapt BIC to further penalize networks with complex structure. There are other approaches that use a non-uniform prior for BIC, notability, the extended BIC (EBIC) [4], an similar modification to BIC which aims to address the small-n-large-P situation. In section 6.2 we compare both the binomial structure prior and EBIC against the use of a uniform prior.

5.2 Multinomial Scoring with Continuous Parents

Both scores presented in this paper reduce to multinomial scoring in the case of a discrete child with exclusively discrete parents. As an alternative, we explore the use of multinomial scoring when there are discrete children and any combination of parents. Before starting a search, we create discretized versions of each continuous variable using equal frequency binning with a predefined number of bins b. Whenever scoring a discrete child, we replace any continuous parents with the precomputed discretized versions of those variables. This allows us to quickly and efficiently perform multinomial scoring for all discrete children. We will henceforth refer to this adaptation as the discretization heuristic and report our finding when choosing b = 3 as a modification to CG in section 6.4.

5.3 As a Conditional Independence Test

We can readily adapt CG and MVP to produce conitional independence tests; to do so, we calculate the log likelihood and degrees of freedom as usual, but perform a likelihood ratio test instead of scoring with BIC. Suppose we wish to test Y₀ ⫫ Y₁|Z where Y₀ and Y₁ are variables (nodes) and Z is a conditioning set of variables. Define ℓ₀ and df₀ respectively as the log likelihood and degrees of freedom for Y₀ given Pa₀ where Pa₀ = Y₁ ⋃ Z. Further, define ${\mathit{\ell }}_0^{\prime }$ and $d{f}_0^{\prime }$ respectively as the log likelihood and degrees of freedom for Y₀ given $P{a}_0^{\prime }$ where $P{a}_0^{\prime }=Z$. Perform a likelihood ratio test with test statistic $2({\mathit{\ell }}_0-{\mathit{\ell }}_0^{\prime })$ and $d{f}_0-d{f}_0^{\prime }$ degrees of freedom. This tests whether the model encoding Y₀ ⫫ Y₁|Z or the model encoding Y₀ ⫫ Y₁|Z fits the data better. If the scoring method used is not score equivalent, then we must also perform a likelihood ratio test with test statistic $2({\mathit{\ell }}_1-{\mathit{\ell }}_1^{\prime })$ and $d{f}_1-d{f}_1^{\prime }$ degrees of freedom where Y₀ and Y₁ are swapped. In this case we decide the variables are dependent if there is enough evidence in either test to support that hypothesis.

6.1 The Conditional Gaussian Approximation

We empirically evaluated the choice of approximating a mixture of Gaussians with a single Gaussian for CG (Assumption 3) in Table 1. We denote the use of a single Gaussian as Approx and the use of the correct mixture calculation as Exact. Originally the results did not appear comparable as the approximate method output a much denser graph than the exact method. In the results shown, we use the binomial structure prior proposed in section 5.1 and achieve comparable results. We see that the approximation performs better in term of precision and comparably in term of recall when compared to the exact method. In the comparisons, we simulate graphs of average degree 2 and 4 with 200 and 1,000 samples and 100 measured variables using fGES. Results are given with the binomial structure prior adustment set to 1.

6.2 Binomial Structure Prior

We tested the usefulness of the binomial structure prior by simulating 200 and 1,000 samples from graphs of average degree 2 and 4 with 100 measured variables using fGES. We compare our scoring functions with and without the binomial structure prior. Additionally we compare against extended BIC (EBIC). In these experiments the binomial structure prior is set to 1 and EBIC’s gamma parameter is set to 0.5 upon suggestion of the authors [4]. In Tables 2 and 3 report findings when the average degrees of the graphs are 2 and 4 respectively.

While we set the binomial structure prior’s parameter to 1 for the experiments presented in this paper, it is important to note that this parameter can be chosen to be any value greater than 0. By varying the expected number parents, we can influence how sparse or dense the output graph will be. The choice of a low value results in a relatively sparse graph and a high value in a denser one.

From Table 2 and 3, for both the binomial structure prior and EBIC, we see boosts in precision with a reduction in recall. Additionally, we see vast reductions in the computation times. In general, EBIC seems to work better with small sample sizes. This makes sense, since EBIC is aimed at the small-n-large-P situation. However, for 1,000 samples, we find the binomial structure prior relatively well. We use the binomial structure prior for the remainder of our score based experiments.

6.3 Conditional Independence Tests

We tested the usefulness of the CG and MVP scores as conditional independence tests by simulating, 200 and 1,000 samples from graphs of average degree 2 and 4 with 100 measured variables. As a search algorithm, we used CPC Stable [19], which is a modified version of PC [26] that treats ambiguous triples as non-colliders. For independence testing, we set the significance level α = 0.001. Here we also use the discretization heuristic with b = 3 for CG, denoted CGd, however we do not use a structure prior since we are no longer scoring a full Bayesian network in this paradigm. We did not include results for a version of MVP which uses the discretization heuristic because it had little effect. The results are shown in Table 4.

In general, we find that our methods perform better as scores, but still perform reasonably well as conditional independence tests. This is promising for use in algorithms, such as FCI, that model the possibility of latent confounding [26].

6.4 Tests Against Baseline Scores

We used two simple baseline scores as a point of comparison for our methods. The first, which we denote MN, uses multinomial scoring for all cases. In order to do so, we essentially extend the discretization heuristic to the continuous child case so that we are always scoring with a multinomial. The second, which we denote as LR, uses partitioned linear regression in the continuous child case and partitioned logistic regression in the discrete child case. In our experiments, we applied Lib Linear [7], a widely used and efficient toolkit for logistic regression which uses truncated Newton optimization [9]. In a recent paper, Zaidi et al. [27] note that among the many optimization methods that have been evaluated, the truncated Newton method has been shown to converge the fastest, which provides support that Lib Linear is a competitive, state-of-the-art method to apply in our evaluation, as a baseline point of comparison. As with MVP, the appended term on LR denotes the maximum polynomial degree of the regressors.

We compared CG, CGd, MVP 1, LR 1, MVP log n, LR log n, and MN by simulating 200 and 1,000 samples from graphs of average degree 2 and 4 with 100 measured variables. As a search algorithm, we again used fGES. Here we also use the discretization heuristic with b = 3 for CGd and the binomial structure prior set to 1 for all scores. Additionally, boldface text highlights the best performing score for each statistic in each column of the table. The results are shown in Tables 5 and 6. For the results in Table 6, we extended our method for simulating data. Since MVP is designed to handle non-linearity, while CG is not, we modified the continuous child phase of data generation to allow for non-linearities. To do so, we additionally generate second, and third order polynomial terms. However, because of the nature of these non-linear functions, the values of the data often become unmanageably large. To correct for this issue, we resample a variable with square-root and cube-root relationships if the values are too large. Appendix B contains details about how the data were simulated and the parameters used.

Table 5 shows the results when using linearly generated data and 100 variables. As a general pattern, MN had better precision than the CG methods which had better precision than the MVP and LR methods. For recall, just the opposite pattern tended to occur. In terms of timing, in general, MN was faster than the CG methods, which were faster than the MVP methods, which were considerably faster than LR.

Table 6 shows the results when using non-linearly generated data and 100 variables. MN tended to have a higher precision than the MVP and LR methods, which often had higher precision than the CG methods. The relatively good performance of MN is surprising; although multinomial distributions can represent non-linear relationships, the process of discretizing continuous variables loses information; the manner in which we generated the data (see the beginning of section 6) when there is a discrete child and continuous parents may play a role in producing this result. The relatively better precision performance of the MVP methods compared to CG methods is not surprising, given that MVP can model non-linear relationships and CG cannot. In terms of recall, MVP and CG performed comparably, while both performed better than MN. The relative timing results in Table 6 are similar to those in Table 5.

In Tables 5 and 6, there is almost no difference in precision and recall performance between MVP and LR. This result is understandable, since MVP is using an approximation to logistic regression in the case of a discrete child with continuous parents and performing all other cases identically. However, MVP is often 10-fold or more faster than LR.

Table 7 shows the results of assessing the scalability of the methods. We simulated linear data on 500 variables. For average degree 4, no MVP results are shown because our machine ran out of memory while searching. Also, LR is not included at all in Table 7, because LR (as implemented) cannot scale to networks of this size due to time complexity. Table 7 shows that the CG methods had similar precision to MN, which generally had better precision than MVP. For the results shown, the recall of the MVP and CG methods were similar, which were generally better than the recall for MN. MN and the CG methods had similar timing results, which were faster than those of MVP.

In Table 7, we see that the CG and MVP methods are capable of scaling to graphs containing 500 measured variables, albeit sparse ones. CG was able to scale to a slightly denser graph of 500 variables. In general, we see the same performance on these larger networks as before on the networks of 100 measured variables. Additionally, for the smaller sample size of 200, MN performed comparably to CDd, but with a slightly higher precision and lower recall.