# Chain Graph Models to Elicit the Structure of a Bayesian Network

Authors: Federico M. Stefanini
Journal: The Scientific World Journal
 (2014)

---

## Abstract

Bayesian networks are possibly the most successful graphical models to build decision support systems. Building the structure of large networks is still a challenging task, but Bayesian methods are particularly suited to exploit experts’ degree of belief in a quantitative way
while learning the network structure from data. In this paper details are provided about how to build a prior distribution on the space of network structures
by eliciting a chain graph model on structural reference features. Several structural features expected to be often useful during the elicitation are described.
The statistical background needed to effectively use this approach is summarized, and some potential pitfalls are illustrated. Finally, a few seminal contributions from the literature are reformulated in terms of structural features.


---

## Body

## 1. Introduction

Bayesian networks (BNs) are possibly the most successful graphical models to represent probabilistic and causal relationships [1, 2]. BNs are used in very different fields including medical domains [3], engineering [4], ecology [5], bioinformatics [6], and many others.The core of this class of models is made by a directed acyclic graph (DAG)𝒢, where nodes in the graph are labels of modeled variables (elements of vector X), and oriented edges (arrows) capture probabilistic and/or causal relationships. The joint distribution of X is represented by the product of conditional distribution functions following from the structure of 𝒢. If substantial prior information is available on a given problem domain, it is possible that an expert defines the structure of 𝒢 and even the parameters inside the conditional distribution functions at a reasonable extent. Otherwise, structure and parameters have to be estimated from data (structural and parameter learning), for example, from a collection of exchangeable observations 𝒟=(x1,…,xn).It is often the case that an expert knows some features of DAG𝒢, but knowledge is not enough to define a DAG because several of its aspects are affected by relevant uncertainty. Following the Bayesian paradigm [7, 8], the expert is invited to state his/her degree of belief about 𝒢 by eliciting a prior distribution on the set of DAGs which can be considered given a fixed set of variables (nodes). All other unknown quantities, like missing values and parameters, are considered in the prior distribution [9].Learning the structure of a BN remains nowadays still challenging for the combinatorial explosion of candidate structures with the increase in the number of considered nodes. Following Robinson [10], the number of possible DAGs on 6 nodes is 3781503; thus the enumeration of all structures while eliciting expert beliefs is unfeasible. Computational difficulties in the full enumeration of DAGs follow from about a dozen of nodes on. For these reasons, several restrictions and simplifications in stating prior beliefs were considered in the past, with the aim of making structural learning tasks affordable in large networks. Widely adopted elicitation techniques are based on restrictions like a total ordering of nodes, the presence of sharp order constraints on nodes, the marginal independence of unknowns, or the existence of a prior network which is a good summary of expert beliefs.In a recent work [11], graphical models were proposed to elicit beliefs about the structure of a BN. The approach is characterized by the possibility of expressing beliefs about limited (but relevant) aspects of the problem domain, called structural features [12], that by their own nature are often only indirectly related to oriented edges in the unknown DAG.Here the most general approach based on chain graph (CG) models is reconsidered with the aim of establishing connections with seminal contributions from the literature on structural learning, of detailing a general parameterization, and of describing one way to perform the refinement of elicited beliefs. Common useful structural features are defined and some issues related to their implementation and revision are examined.


## 2. Methods

Notation and some background information on graphs and Markov properties are provided below. A comprehensive account may be found in [13–15]. An approach to the elicitation of structural features by CG models is described, together with methods for the revision.


### 2.1. Graphs

A graph𝒢 is a pair (V,E) where V={v1,v2,…,vK} is a finite set of nodes and E⊂V×V is the set of edges. The set E represents the structure of the graph because it defines which nodes are linked by an edge and if such edge is directed (arrow) or not. If (vi,vj)∈E and (vj,vi)∈E then an undirected edge joins vi and vj, indicated as vi—vj, and vj is neighbor of vi. If (vi,vj)∈E but (vj,vi)∉E the ordered pair corresponds to the directed edge vi→vj; vj is said to be the child of vi and vi is a parent of vj. The set pa(vj) includes all parents of node vi, that is, all nodes originating arrows with end in vi, while the set ch(vi) collects all children of vi, namely, all nodes in which arrows originated from vi end.A path is a sequence of vertices such that there is an edge for each pair of subsequent nodes in the sequence, that is,vi—vi+1 or vi←vi+1 or vi→vi+1. A directed path is a path in which all edges maintain the head-to-tail orientation, for example, (vi,vj,vk) with vi→vj→vk.In a directed graph all edges are directed. The ancestorsan(vi) of node vi are nodes on a directed path reaching vi, while descent nodes de(vi) are nodes on a directed path starting from node vi. Note that vi∈de(vi) and that vi∈an(vi). The extension of the above definition to an(A) with A⊂V is obtained by union of sets an(vi) for each vi∈A. A similar extension holds for de(A).In a directed graph without cycles it is not possible to visit the same node more than one time by following a directed path, and in this case the graph is called directed acyclic graph. A moralized DAG is an undirected graph obtained by joining pairs of nodes sharing children (if not yet connected) with an undirected edge and by removing the direction of all edges. A subgraph𝒢A on A⊂V is obtained by removing all nodes in V∖A and all edges in which at least one node is in V∖A from the graph 𝒢.A graph without directed edges is called undirected graph (UG). An UG withE=V×V is said to be complete. A subgraph on a subset S⊂V of nodes is obtained by removing nodes not in S and all edges reaching-leaving nodes not in S. Note that a subset of nodes is indicated by capital letters. A clique C is a maximal complete subgraph of an UG.A chain graph, also called partially directed acyclic graph (PDAG), is made by an ordered collection of chain componentsτ=(τ1,τ2,…) which are undirected graphs and by directed edges between nodes located in different chain components, so that the arrow vi→vj is allowed only if vi belongs to a chain component preceding the chain component in which vj is located. Therefore directed edges are forbidden within a chain component and in the direction from τj to τj-k, with k>0.The moralization of a chain graph mimics the moralization of a DAG. A moralized CG, indicated as𝒢m, is obtained by the following steps:(1)
letk=2;(2)
join with undirected edges all pairs of nodes inpa(τk), with pa(τk) being the union of parents sets for each node in chain component τk;(3)
iterate step (2) fork=3,4,…;(4)
remove directions to all edges.


### 2.2. Some Markov Properties

Conditional independence [16] is fundamental to reason out highly structured stochastic systems and to simplify the representation of high dimensional distributions.In this paper the random vectorX refers to random variables included in the BN. The notation used hereafter is based on nodes in V; thus X is also indicated as XV=Xv1,…,vK=(Xv1,…,XvK) and the sample space ΩXV=⊗vi∈VΩXvi is the Cartesian product of sample spaces of considered variables.The joint probability distribution of a random vectorXV is Markov with respect to an UG 𝒢 if
(1)p(xv1,xv2,…,xvK)=σ-1∏C∈𝒞‍ϕC(xC)
with 𝒞 being a collection of graph cliques in 𝒢 and with ϕC nonnegative functions called clique potentials; note that xC is the coordinate projection of vector xV on the subset of coordinates defined by C; σ=∑ΩXV‍∏C∈𝒞‍ϕC(xC) is the partition function that normalizes the product of potentials in (1).Markov properties for positive distributions with respect to an undirected graph may be read using the separation theorem [14]. Let A⊂V, B⊂V, S⊂V be disjoint subsets of nodes. The separation theorem states that subvectors XA and XB are conditionally independent given the subvector XS if and only if all paths from a node in A to a node in B include nodes located in S; thus nodes in S separate nodes in A from nodes in B.The joint probability distribution of random variables indexed inV is Markov with respect to a DAG 𝒢 if the following factorization holds:
(2)p(xv1,xv2,…,xvK)=∏vi∈V‍p(xvi∣xpa(vi)),
where xpa(vi) is the random vector made by variables whose labels belong to the parents set of vi.Markov properties on a DAG may be read using the separation theorem for directed graphs [14]. Given three disjoint sets of nodes A,B,S, consider the subgraph 𝒢 defined on an(A∪B∪S) after moralization, say 𝒢an(A,B,S)m. Random subvectors XA and XB are conditionally independent given the subvector XS if and only if nodes in S separate nodes in A from nodes in B in 𝒢an(A,B,S)m.A joint probability distribution is Markov with respect to a CG𝒢 if
(3)p(xv1,xv2,…,xvK)=∏τi∈τ‍p(xτi∣xpa(τi))
with τi∈τ, the chain components of 𝒢. Furthermore factors on r.h.s. of (3) may be factorized by considering the subgraph on nodes defined by τi∪pa(τi):
(4)p(xτi∣xpa(τi))=σt-1∏C∈𝒞‍ϕC(xC),
where σt-1, t∈Ωpa(τi) are normalization constants, one for each conditioning value of the random subvector Xpa(τi). The working UG in (4) is obtained by removing the orientation of edges from parents of nodes in τi and by joining them into a complete undirected subgraph.Conditional independence relationships in CGs are also obtained from an extension of the separation theorem in UGs and DAGs for positive distributions [17]. Let 𝒢 be a chain graph and A⊂V, B⊂V, S⊂V be three disjoint subsets of V. The separation theorem states that subvectors XA and XB are conditionally independent given the subvector XS if and only if all paths from a node in A to a node in B in 𝒢an(A∪B∪S)m include nodes of S; thus nodes in S separate nodes in A from nodes in B. Note that 𝒢an(A∪B∪S)m is the moral graph of the smallest ancestral set for A∪B∪S, that is, a subgraph of 𝒢m described in Section 2.1.


### 2.3. Causal DAGs

A DAG may represent causal relations among variables. According to the causal semantic, an arrowvi→vj indicates that vi is a direct cause of vj with respect to nodes included in V, that is, at the considered model granularity. In principle the intervention on variable xvi may bear an effect on xvj. The intervention on a subset of variables D⊂V indicates the external setting of variables in XD to prescribed values; thus the system or process is perturbed, not merely observed.Pearl [2, pp. 27–32] starts with the definition of functional causal models, which are deterministic in nature, and he demonstrates [2, theorem 1.4.1] that such formulation induces the Markov factorization in (2), the so-called Markov causal assumption. An equivalent representation embeds exogenous variables into the node of interest and transforms the deterministic relationships into probabilistic conditioning, thus leading to Bayesian networks.A key property of a casual DAG𝒢 is the stability under external intervention: if a variable xvi is manipulated all the other variables maintain their relationships as represented by 𝒢. In other terms the intervention is local on manipulated variables and it does not break all the other relationships represented in a causal DAG. The intervention regime is in contrast with the plain observation of values taken by the random vector XV.The granularity of a causal DAG depends on the variables included in the model. A variableXL not included in a causal DAG may eventually affect just one variable Xv, with v∈V; otherwise if several variables are affected then a more general class of models is needed, called semi-Markovian networks (not considered further in this work).For an updated presentation of the approach see [18] while [19] warns against the blind definition of DAGs within the causal semantic in observational studies. He also reconsiders the foundations of causal inference by anchoring them to the extended conditional independence (C.I.), both at algebraic level and with a graphical counterpart based on influence diagrams.


### 2.4. Inference about the Structure of Bayesian Networks

In all cases in which strong prior information is absent, structural learning is performed by means of a database𝒟=(x1,…,xn) of n exchangeable observations of the random vector XV.Several algorithms have been proposed to infer causal and probabilistic relations, but in Bayesian inference key quantities enter into the joint probability distribution of𝒟 and network’s unknowns given the context ξ:
(5)p(𝒟,θ,z∣ξ)=p(𝒟∣θ,z,ξ)·p(θ∣z,ξ)·p(z∣ξ),
where θ=(θv1,pa(v1),…,θvK,pa(vK)) are vectors of parameters characterizing the conditional probability distributions of Xvi given Xpa(vi); variable Z indicates the unknown DAG, and it is built as a bijection from the set of DAGs (fixed V) to a subset ΩZ of natural numbers.The likelihood functionp(𝒟∣θ,z,ξ) is typically expressed as a product of multinomials by using sufficient statistics. Whenever expert beliefs are reasonably captured by Dirichlet prior distributions for elements of θ, and they are elicited as marginally independent, closed-form integration marginalizes out thetas. The resulting marginal distribution p(𝒟∣z,ξ) has a reduced dimensionality and may be optimized with respect to z while looking for optimal structures characterized by the highest posterior probability values [9]:
(6)p(𝒟,z∣ξ)=∫‍p(𝒟∣θ,z,ξ)·p(θ∣z,ξ)·p(z∣ξ)·dθ=p(z∣ξ)·∏i=1K‍∏j=1qi‍Γ(αi,j)Γ(αi,j+Ni,j)·∏s=1ri‍Γ(αi,j,s+Ni,j,s)Γ(αi,j,s)
with ri being the number of states taken by Xi and qi of Xpa(vi); sufficient statistics are Ni,j,s for the ith variable taking the sth state while its parent configuration is in the jth state; αi,j=∑s=1ri‍αi,j,s and Ni,j=∑s=1ri‍Ni,j,s. A clever choice of hyperparameters guarantees the likelihood equivalence; that is, all BNs equivalent as regards the set of encoded conditional independence relationships are equally scored by (6).Equation (6) shows that the elicitation of the prior distribution p(z∣ξ) is a key step to perform Bayesian structural learning using 𝒟.


### 2.5. Plausible Network Features

A candidate structurez∈ΩZ on a fixed set of nodes V is plausible if it has got structural features (SFs) believed to be relevant by the expert. Following Stefanini [11, Definition  1], we recall the formal definition.Definition 1 (SF).
A structural feature (SF)ℛj(z,w) in a reference set ℛ for the set of DAGs on V is a predicate describing a plausible probabilistic or causal characteristic of the unknown directed acyclic graph z∈ΩZ. Argument w is in the partition 𝒲 of a given numeric domain ΩW of variable W. An atomic structural feature (ASF) ℛj(z) does not depend on any auxiliary variable W.A reference setℛ is a collection ℛ={ℛj:j∈J} of SFs indexed in a set J, with nf being the number of considered SFs. A proposition might be defined to carry disbelief to a candidate structure, but the feature-rises-belief direction is conveniently adopted here.An example makes the previous definition operational. Let us defineℛ1(z,w)= “The number of immediate causes of Xv3 is in w” to capture the expert belief about the number of parents of node v3 in the unknown DAG. An expert may consider ΩW={0,1,…,12} and 𝒲=({0},{1,2,3},{4,5},{6,…,12}). Given a candidate structure z, its plausibility will be determined by the element w∈𝒲 that makes the predicate true. A simple representation of configurations taken by reference features is obtained by descriptors [11, Definition  2].Definition 2 (descriptors).
A descriptorRi for the SF ℛi is a map:
(7)Ri:𝒲i×{false,true}⟶{0,1,2,…,|𝒲i|}
so that (w,false)↦0 for all w∈𝒲i and (w,true)↦hw for all w∈𝒲i; that is, a different integer is associated with each w if true. The vector R=(R1,R2,…,Rnf) defined on the Cartesian product ΩR=⊗i=1nfΩRi=⊗i=1nf{0,1,…,|𝒲i|} is called vector of descriptors. The descriptor of an ASF is defined by false↦0 and true↦1.Thejth configuration of descriptors in R is indicated as rj=(r1,j,r2,j,…,rnf,j)∈ΩR while a generic configuration is indicated as r∈ΩR. The jth configuration of a subvector of R defined by indexes in A is rA,j; for example, r{1,3},j=(r1,j,r3,j).VectorR induces equivalence classes on ΩZ; that is, 𝒵={𝒵r:r∈ΩR} contains sets of structures, with 𝒵r made by all those DAGs sharing the same configuration r.Note that members of the same equivalence class must be associated with the same degree of belief by the principle of insufficient reason: they differ in irrelevant and unconsidered ways by construction.Despite the generality of Definition1 some features are expected to occur more often than others. Below, some of them are described without pretending to be exhaustive.


#### 2.5.1. Indegree and Outdegree

In applications characterized by a small sample size, it is useful to impose a sharp constraint during the greedy search of top-scored candidate structures. The maximum number of arrows entering into a node, the indegree, is set to a small integer, say2 to 5, to exclude structures with a large number of parents from the consideration. A large number of parents implies a CPT with a huge number of parameters which are affected by large uncertainty after conditioning to observed data because of sampling zeros. A similar constraint may be set on the number of arrows leaving a node.Two reference features naturally embed this kind of information.(i)
“The maximum number of arrows reaching a node is innid for all nodes in V,” where nid is a small set of integers close to 1 elicited from the expert.(ii)
“The maximum number of arrows leaving a node is innod for all nodes in V,” where nod is a small set of integers close to 1 elicited from the expert.The above two features may be exploited to increase the plausibility of candidate structures which are sparsely connected. Higher control on connectivity is obtained by considering the fraction of nodes showing a given degree of connectivity, as described below.


#### 2.5.2. Partitioned Connectivity

A different way of characterizing the connectivity is obtained by eliciting the minimum fraction of nodes in DAG showing a given number of parents (children) and by iterating the elicitation from0 parents (children) up to a small integer s.Lets be a small integer representing the maximum number of parents (children) to be considered. Let W=(W0,…,Ws) be a vector of nondecreasing numbers in [0,1], with
(8)ΩW={(ω0,ω1,ω2,…,ωb,…,ωs):0hhh≤ωb≤ωb+1<ωs=1}
being the sample space. The elicitation of this feature is based on a vector a=(a0,a1,a2,…,ab,…,as), with 0≤ab≤ab+1<as=1, and on the induced partition 𝒲={w0,w1} where
(9)w1={(ω0,ω1,…,ωs):ωb≥ab,b=0,1,…,s}
and w0=ΩW∖w1.Two reference features naturally embed the extended evaluation of connectivity:(i)
“Thea-partitioned inconnectivity of degree s is in w1,” with a, s, and w1 defined above; a candidate structure z shows this feature if the cumulative relative frequency Fb of nodes in z with number of parents equal or less than b is greater than ab; that is, Fb≥ab, b=1,2,…,s.(ii)
“Thea-partitioned outconnectivity of degree s is in w1,” with a, s, and w1 defined above; a candidate structure z shows this feature if the cumulative relative frequency Fb of nodes in z with number of children equal or less than b is greater than ab; that is, Fb≥ab, b=1,2,…,s.Trained experts could prefer a conventional total number of nodes equal to100 to elicit cumulative percentages instead of cumulative fractions of nodes.In Table1, an example is shown where s=5 and the elicited vector a is defined on fractions. A candidate structure has the partitioned inconnectivity feature if the fraction of root nodes is equal or above 0.01, while the proportion of nodes with at least 1 parent is equal or above 0.2 and so on.Table 1
Elicited vectora, with s=5. In the top row, the number b of parents is reported, while in the second row the correspondent minimum fraction of nodes ab is shown.
Number of parents
0
1
2
3
4
5
Cumulative fraction of nodes
0.01
0.2
0.4
0.6
0.8
1.0


#### 2.5.3. Direct Cause and Direct Effect

The reference feature “The variablexvi is an immediate cause of variable xvj” refers to xvi as a parent of xvj so that by setting (intervening on) the value of the variable xvi to a given value, the distribution of Xvj is modified. This relation holds at the level of selected granularity; thus it may change if the collection of variables (nodes in V) is modified.


#### 2.5.4. Causal Ancestors

The “direct cause” feature may be extended by considering a variablexvi which is on a causal (directed) path reaching node xvj. In this case the reference feature is “The variable xvi is an indirect cause of variable xvj”; thus the expert believes that one or more variables mediate the effect of xvi on xvj.


#### 2.5.5. Causal Hubs

A hub node in a network is characterized by a high number of arrows leaving it. A reference feature which captures the local connectivity of nodevi is “Node vi is a hub node of at least outdegree w,” with the outdegree indicating the number of arrows originated in vi and w a set of integers. Note that the defined feature is a localized version of the outdegree feature.The expert might believe that a hub node should be present, but without indicating a specific node. In this case thea-partitioned outconnectivity feature (with a large s) can be exploited for this purpose.


#### 2.5.6. Conditional Independence Relationships

A statement about C.I. among three disjoint subsets of random variables may take the following form: “The random vectorXA is conditionally independent from the random vector XB given vector XS,” with A,B,S being disjoint subsets of nodes in V.


### 2.6. The Degree of Belief

The prior distribution on the space of structures onV is obtained by “extending the argument;” that is,
(10)P[Z=z∣ξ]=∑r∈ΩR‍P[Z=z∣R=r,ξ]·P[R=r∣ξ].
By recognizing that R induces the partition 𝒵, it follows that
(11)p(z∣ξ)=1nr[z]·P[R=r[z]∣ξ]=p(r1[z],r2[z],…,rnf[z]∣ξ)nr[z],
where nr[z] is the cardinality of the equivalence class in which z is located and where the structural configuration of DAG z is r[z]. In [12] the size of each equivalence class 𝒵r was estimated by Monte Carlo simulation to face the combinatorial explosion in the number of DAGs to be enumerated with the increase of the number of nodes:
(12)|𝒵^r|=NV(Nr+1)NT+1,
where NT is the total number of DAGs uniformly sampled from the space of all DAGs on V, NV is the size of such space (see [10]), and Nr≤NT is number of sampled DAGs showing configuration r.The numerator on the right of (11) represents the joint belief on each configuration of descriptors. While the elicitation in full generality becomes cognitively and numerically overwhelming around 7 descriptors on, some parsimony is achieved if a small number of descriptors may be considered at one time, so that conditional independence relationships among descriptors may be exploited. This is a choice available to the expert through the definition of an order relation on descriptors.Definition 3 (ordered partition).
The ordered partition(13)𝒪=(𝒪1,𝒪2,…)
of descriptors is defined by the expert to indicate disjoints subsets of SFs to be jointly considered during the elicitation, from the first subset 𝒪1 to the last.Otherwise stated, the expert decomposes the whole elicitation problem following an order taken from the substantive content of the specific problem domain: features are grouped according to the priority in the elicitation.The elicitation of the ordered partition is performed by formulating questions in the language typical of a given problem domain. It is indeed difficult to define those questions in general terms, because they result to be quite abstract and they are likely to be obscure for the domain expert (see [11]). We assume here that questions are properly phrased and that an ordered partition is defined.If a strict order relation is elicited, each element of𝒪 contains one descriptor: this is a special case addressed in [20]. Another special case is represented by a trivial partition of just one subset that contains all descriptors [21]. In the two special cases above it is possible to define, respectively, a Bayesian network and a Markov network on descriptors. The general case made by several subsets, each one of cardinality two or more, was addressed by using CG models in [11] for ASFs. Below details are provided about a reference set of nonatomic features.The joint probability distribution of descriptorsp(r1,r2,…,rnf∣𝒪,ξ) is assumed to be Markov with respect to the elicited partition where descriptors within group 𝒪j define the chain component τj of the CG model (see (3)). The elicitation keeps on by asking in the language of the domain expert which descriptors in subset 𝒪j should be jointly considered while defining the degree of belief. Nodes corresponding to related descriptors are joined by undirected edges and the collection of cliques defining the chain component are found. The elicitation is iterated for all elements in the ordered partition 𝒪.The resulting CG may support the elicitation if model parameters are cognitively suited for the quantitative step, that is, if they are interpretable and easy to assess for the expert, at least after some training.The first chain component is elicited as a marginal distribution which is not conditioned on other descriptors. An undirected graphical model on descriptors in𝒪1 is defined through the multiplicative model in (4), under empty conditioning. Nevertheless some care is needed because potential functions on cliques are not uniquely defined. Here a log-linear parameterization is suggested, following [13]. For example, if the first chain component has just one clique made by three descriptors, say R1,R2,R3, we define the multiplicative model as below, after exponentiation and rearrangement of log-linear terms:
(14)p(r1,r2,r3∣ξ)p(0,0,0∣ξ)=ϕ1(r1)ϕ2(r2)ϕ3(r3)ϕ1,2(r1r2)×ϕ1,3(r1r3)ϕ2,3(r2r3)ϕ1,2,3(r1r2r3).
Thus the odds value with respect to the no-feature configuration r1,2,3=(0,0,0) is explained by a multiplicative model where each factor, for example, ϕ2,3(r2r3), is equal to one if one or more descriptors are null, otherwise they are positive (the so-called treatment parameterization).The elicitation in (14) is performed on the odds scale by asking the domain expert how many times the configuration (r1,r2,r3) is more plausible than (0,0,0) for descriptors R1,R2,R3. This question is iterated for all configurations in the Cartesian product ΩR1×ΩR2×ΩR3 after exclusion of (0,0,0). Questions are posed from single main effects towards higher order interactions terms, so that one factor at a time is considered (see algorithm below).The procedure is iterated for all cliques in the first CG componentτ1, and indeed factors already elicited in previously considered cliques are not reconsidered anymore. For example, the chain component R1—R2—R3 is made by two cliques, R1,R2 and R2,R3; thus the shared factor ϕ2(r2) is elicited just one time.The general algorithm for the first chain component is summarized below.(1)
Consider the undirected graph on descriptors elicited as the first chain component. Control questions include the following: “Are all the relevant features included in the elicitation?”, “Are all pairs of features jointly affecting the probability of a structure linked by an undirected edge?”. If needed, revise the order relation and (or) links within chain component.(2)
Check out that each descriptorRvi given its neighbors Rne(vi) is independent of all other descriptors in the first CG component.(3)
Find the cliques of such UG (model generator) [13].(4)
For each clique, elicit parameter values by using odds.
(i)
Elicit main effects{ϕi(ri,j)}, one at a time for all configurations, by assigning an odds value with respect to the baseline with no features; that is,
(15)ϕi(ri,j)=p(0,…,0,ri,j,0,…,0∣ξ)p(0,…,0∣ξ).A control question for this step is: “How much above one is the odds value for the sole presence of featureRi=ri,j>0?”(ii)
Elicit first order interactions{ϕi,l(ri,j′,rl,j′′)}, one at a time, by assigning a multiplicative term for each configuration under the question: “Which is the value of the multiplicative term {ϕi,l(ri,j′,rl,j′′)} needed to account for the interaction of feature Ri=ri,j′>0 with feature Rl=rl,j′′>0?” The expression helping in this step is
(16)p(0,…,ri,j′,0,…,rl,j′′,0,…,0∣ξ)p(0,…,0∣ξ)=ϕi(ri,j′)ϕl(rl,j′′)ϕi,l(ri,j′,rl,j′′),whereϕi(ri,j′), ϕl(rl,j′′) are already elicited; this is the cross product ratio of features Ri,Rl, after dividing by ϕi(ri,j′), ϕl(rl,j′′). It is clear that if the interaction is absent ϕi,l(ri,j′,rl,j′′)=1, while ϕi,l(ri,j′,rl,j′′)∈(0,1) means that the interaction reduces the plausibility and ϕi,l(ri,j′,rl,j′′)>1 raises the plausibility of the considered configuration.(iii)
Iterate the step above with higher order interaction terms (for all configurations) with two constraints: before moving to higher interaction terms all the terms of the same degree must have been already elicited; moreover the maximum degree of interaction among a subset of features is defined by the size of the clique they belong to (model generator). For example, with the interaction of order two{ϕi,l,k(ri,j′,rl,j′′,rk,j′′′)}, the control question might be: “Having already elicited values of ϕi(ri,j′), ϕl(rl,j′′), ϕk(rk,j′′′), ϕi,l(ri,j′,rl,j′′), ϕi,k(ri,j′,rk,j′′′), and ϕl,k(rl,j′′,rk,j′′′), which is the value of the multiplicative term {ϕi,l,k(ri,j′,rl,j′′,rk,j′′′)} needed to adjust the odds value after considering the interaction among features whose configurations are ri,j′>0, rl,j′′>0, rk,j′′′>0?”; the helper expression is
(17)p(0,…,ri,j′,…,rl,j′′,…,rk,j′′′,…,0∣ξ)p(0,…,0∣ξ)=Qϕi,l,k(ri,j′,rl,j′′,rk,j′′′),where factor(18)Q=ϕi(ri,j′)ϕl(rl,j′′)ϕk(rk,j′′′)ϕi,l(ri,j′,rl,j′′)×ϕi,k(ri,j′,rk,j′′′)ϕl,k(rl,j′′,rk,j′′′)is already elicited. Similar expressions follow straightforwardly for higher order interactions.(5)
Calculate the baseline probabilityp(0,…,0∣ξ) of a structure without features by exploiting ∑r∈ΩR‍p(r∣ξ)=1.(6)
Revise the elicited values following the guidelines of Section2.7.(7)
Move to the next chain component.The algorithm presented above may be also applied to chain componentsτ2,τ3,…, that is, with a nonempty set of parents. The elicitation of parameters in a conditional Markov network is performed by iterating the algorithm for the first chain component over each configuration of conditioning parents. It follows that the elicitation burden depends on the number of parents, more precisely, on the cardinality of the Cartesian product where factors are samples spaces of parents.The algorithms defined above lead to a prior distribution on configurations, but the elicitation does not end before the revision of elicited values takes place.


### 2.7. Implementation and Revision Issues

The flexibility achieved by defining predicates entails potential pitfalls that should be considered during an elicitation run with structural features.The revision of elicited beliefs is always needed because the expert is not expected to provide unbiased elicited values, especially with limited training or in complex problem domains. Revision and elaboration of the prior distribution [22, and references therein] are applied to control elicitation bias and other causes of poor elicitation.Overelicitation is an important practice to check for the presence of elicited quantities that do not reflect the expert degree of belief. Nevertheless, in large networks elaboration of elicited probability values is the main resource for checking the quality of elicitation. In the proposed framework, the elicitation through reference features produces proper probability distributions; therefore one elaboration consists of inspecting one or more margins of the probability distributionp(r1,…,rnf∣ξ), for example, the bivariate margin p(ri,rj∣ξ), to look for configurations whose plausibility causes surprise or disbelief in the expert. This operation is particularly meaningful if inspected margins do not straightforwardly relate to elicited odds, for example, taking descriptors belonging to different cliques. Surprise and disbelief ask for the revision of elicited values. If an association between selected margins and bias is suspected, random selection of margins is an option.The stability of elicited values against different order relations on descriptors (reference features) should not be assumed. Nevertheless, in complex problem domains overelicitation made by the repetition of the interview with other ordered partitions not only seems unacceptable as regards the work load, but it could even be cognitively unfeasible, for example, if the original ordered partition selected by the expert is induced by a scientific hypothesis.The core of the proposed approach is based on predicates representing structural features. The expert might believe that some configurations of features are plausible although they are incompatible with DAGs, for example, because they imply the presence of cycles. A positive probability value would be assigned to such a configurationr*, but this is not an instance of elicitation bias if the elicited distribution properly matches expert beliefs. If this is the case, p(z∣ξ)=0 because P[Z=z∣R=r*]=0.The way a predicate is specified determines the granularity of the elicitation and the cardinality of equivalence classes in𝒵. For example, let us consider two nodes vi∈V and vj∈V and the reference feature ℛ1= “Nodes vi, vj are not descendent of other nodes in VR”.Reworking the original proposition we have two simpler predicates:(i)
ℛ1′= “Node vi is not a descendent of other nodes in VR”;(ii)
ℛ1′′= “Node vj is not a descendent of other nodes in VR.”and they can be considered by conjunction. In Table 2, the relation between the original reference feature (right) and the conjoint components (left) is shown: ¬ℛ1 collects three configurations generated by the conjunction of simpler predicates. It follows that simpler predicates are needed if the expert degree of belief changes over collapsed configurations in ℛ1, that is, ¬ℛ1′∩ℛ1′′, ℛ1′∩¬ℛ1′′, ¬ℛ1′∩¬ℛ1′′. While nothing prevents the expert from defining rich predicates, care should be taken to select a granularity suited to properly represent expert beliefs.Table 2
Relation between the original feature (right) and the conjunction of two simpler structural features (left).
Simpler feature1
Simpler feature2
Original feature
ℛ
1
′
ℛ
1
′
′
ℛ
1
¬
ℛ
1
′
ℛ
1
′
′
¬
ℛ
1
ℛ
1
′
¬
ℛ
1
′
′
¬
ℛ
1
¬
ℛ
1
′
¬
ℛ
1
′
′
¬
ℛ
1There are indeed several different ways of formulating a predicate. If two reference sets of features induce the same set of equivalence classes𝒵 then they are operationally equivalent. Nevertheless, from a cognitive standpoint, substantial differences in the ease of elicitation might depend on the way propositions are formulated [22]. Let us consider the reference feature ℛ1= “Nodes vi,vj precede all other nodes.” Does the expert use “precede” as “come before” all other nodes in the order relationship of nodes in V? Given a DAG with vj disconnected from other nodes how to answer? Does the expert use “precede” meaning “are ancestors of all other nodes in V” or to indicate that those nodes do not receive arrows coming from other nodes? This example makes clear that some training is mandatory in order to make an expert effective in the elicitation: a trained expert is expected to choose the right granularity in the elicitation and to define meaningful predicates, that is, statements straightforwardly true/false when applied to any DAG defined on V.Several reference features are jointly considered in actual applications, a number typically far beyond what the expert may simultaneously consider with success. Thus there is the possibility that exclusion and implication relations are not recognized. For example, let us consider two features:ℛ3= “The indegree is three or less” and ℛ25= “vi is a sink node.” After reworking the last feature, the expert reformulates the statement as “Node vi has indegree ten or more.” Clearly the plausibility of ℛ3∧ℛ25 should be null; otherwise the actual interpretation of ℛ3 is “The indegree of all but vi is three or less.”Implication must be recognized to maintain the probabilistic coherence; for example, letℛ3= “Variable Xvi is a direct cause of Xvj” and ℛ5= “Variable Xvi is an ancestral cause of Xvj” be two reference features. Clearly ℛ3 implies ℛ5 and the joint plausibility of both features is bound to the plausibility of ℛ3. In a normalized reference set logical relationships among features are properly handled.


## 2.1. Graphs

A graph𝒢 is a pair (V,E) where V={v1,v2,…,vK} is a finite set of nodes and E⊂V×V is the set of edges. The set E represents the structure of the graph because it defines which nodes are linked by an edge and if such edge is directed (arrow) or not. If (vi,vj)∈E and (vj,vi)∈E then an undirected edge joins vi and vj, indicated as vi—vj, and vj is neighbor of vi. If (vi,vj)∈E but (vj,vi)∉E the ordered pair corresponds to the directed edge vi→vj; vj is said to be the child of vi and vi is a parent of vj. The set pa(vj) includes all parents of node vi, that is, all nodes originating arrows with end in vi, while the set ch(vi) collects all children of vi, namely, all nodes in which arrows originated from vi end.A path is a sequence of vertices such that there is an edge for each pair of subsequent nodes in the sequence, that is,vi—vi+1 or vi←vi+1 or vi→vi+1. A directed path is a path in which all edges maintain the head-to-tail orientation, for example, (vi,vj,vk) with vi→vj→vk.In a directed graph all edges are directed. The ancestorsan(vi) of node vi are nodes on a directed path reaching vi, while descent nodes de(vi) are nodes on a directed path starting from node vi. Note that vi∈de(vi) and that vi∈an(vi). The extension of the above definition to an(A) with A⊂V is obtained by union of sets an(vi) for each vi∈A. A similar extension holds for de(A).In a directed graph without cycles it is not possible to visit the same node more than one time by following a directed path, and in this case the graph is called directed acyclic graph. A moralized DAG is an undirected graph obtained by joining pairs of nodes sharing children (if not yet connected) with an undirected edge and by removing the direction of all edges. A subgraph𝒢A on A⊂V is obtained by removing all nodes in V∖A and all edges in which at least one node is in V∖A from the graph 𝒢.A graph without directed edges is called undirected graph (UG). An UG withE=V×V is said to be complete. A subgraph on a subset S⊂V of nodes is obtained by removing nodes not in S and all edges reaching-leaving nodes not in S. Note that a subset of nodes is indicated by capital letters. A clique C is a maximal complete subgraph of an UG.A chain graph, also called partially directed acyclic graph (PDAG), is made by an ordered collection of chain componentsτ=(τ1,τ2,…) which are undirected graphs and by directed edges between nodes located in different chain components, so that the arrow vi→vj is allowed only if vi belongs to a chain component preceding the chain component in which vj is located. Therefore directed edges are forbidden within a chain component and in the direction from τj to τj-k, with k>0.The moralization of a chain graph mimics the moralization of a DAG. A moralized CG, indicated as𝒢m, is obtained by the following steps:(1)
letk=2;(2)
join with undirected edges all pairs of nodes inpa(τk), with pa(τk) being the union of parents sets for each node in chain component τk;(3)
iterate step (2) fork=3,4,…;(4)
remove directions to all edges.


## 2.2. Some Markov Properties

Conditional independence [16] is fundamental to reason out highly structured stochastic systems and to simplify the representation of high dimensional distributions.In this paper the random vectorX refers to random variables included in the BN. The notation used hereafter is based on nodes in V; thus X is also indicated as XV=Xv1,…,vK=(Xv1,…,XvK) and the sample space ΩXV=⊗vi∈VΩXvi is the Cartesian product of sample spaces of considered variables.The joint probability distribution of a random vectorXV is Markov with respect to an UG 𝒢 if
(1)p(xv1,xv2,…,xvK)=σ-1∏C∈𝒞‍ϕC(xC)
with 𝒞 being a collection of graph cliques in 𝒢 and with ϕC nonnegative functions called clique potentials; note that xC is the coordinate projection of vector xV on the subset of coordinates defined by C; σ=∑ΩXV‍∏C∈𝒞‍ϕC(xC) is the partition function that normalizes the product of potentials in (1).Markov properties for positive distributions with respect to an undirected graph may be read using the separation theorem [14]. Let A⊂V, B⊂V, S⊂V be disjoint subsets of nodes. The separation theorem states that subvectors XA and XB are conditionally independent given the subvector XS if and only if all paths from a node in A to a node in B include nodes located in S; thus nodes in S separate nodes in A from nodes in B.The joint probability distribution of random variables indexed inV is Markov with respect to a DAG 𝒢 if the following factorization holds:
(2)p(xv1,xv2,…,xvK)=∏vi∈V‍p(xvi∣xpa(vi)),
where xpa(vi) is the random vector made by variables whose labels belong to the parents set of vi.Markov properties on a DAG may be read using the separation theorem for directed graphs [14]. Given three disjoint sets of nodes A,B,S, consider the subgraph 𝒢 defined on an(A∪B∪S) after moralization, say 𝒢an(A,B,S)m. Random subvectors XA and XB are conditionally independent given the subvector XS if and only if nodes in S separate nodes in A from nodes in B in 𝒢an(A,B,S)m.A joint probability distribution is Markov with respect to a CG𝒢 if
(3)p(xv1,xv2,…,xvK)=∏τi∈τ‍p(xτi∣xpa(τi))
with τi∈τ, the chain components of 𝒢. Furthermore factors on r.h.s. of (3) may be factorized by considering the subgraph on nodes defined by τi∪pa(τi):
(4)p(xτi∣xpa(τi))=σt-1∏C∈𝒞‍ϕC(xC),
where σt-1, t∈Ωpa(τi) are normalization constants, one for each conditioning value of the random subvector Xpa(τi). The working UG in (4) is obtained by removing the orientation of edges from parents of nodes in τi and by joining them into a complete undirected subgraph.Conditional independence relationships in CGs are also obtained from an extension of the separation theorem in UGs and DAGs for positive distributions [17]. Let 𝒢 be a chain graph and A⊂V, B⊂V, S⊂V be three disjoint subsets of V. The separation theorem states that subvectors XA and XB are conditionally independent given the subvector XS if and only if all paths from a node in A to a node in B in 𝒢an(A∪B∪S)m include nodes of S; thus nodes in S separate nodes in A from nodes in B. Note that 𝒢an(A∪B∪S)m is the moral graph of the smallest ancestral set for A∪B∪S, that is, a subgraph of 𝒢m described in Section 2.1.


## 2.3. Causal DAGs

A DAG may represent causal relations among variables. According to the causal semantic, an arrowvi→vj indicates that vi is a direct cause of vj with respect to nodes included in V, that is, at the considered model granularity. In principle the intervention on variable xvi may bear an effect on xvj. The intervention on a subset of variables D⊂V indicates the external setting of variables in XD to prescribed values; thus the system or process is perturbed, not merely observed.Pearl [2, pp. 27–32] starts with the definition of functional causal models, which are deterministic in nature, and he demonstrates [2, theorem 1.4.1] that such formulation induces the Markov factorization in (2), the so-called Markov causal assumption. An equivalent representation embeds exogenous variables into the node of interest and transforms the deterministic relationships into probabilistic conditioning, thus leading to Bayesian networks.A key property of a casual DAG𝒢 is the stability under external intervention: if a variable xvi is manipulated all the other variables maintain their relationships as represented by 𝒢. In other terms the intervention is local on manipulated variables and it does not break all the other relationships represented in a causal DAG. The intervention regime is in contrast with the plain observation of values taken by the random vector XV.The granularity of a causal DAG depends on the variables included in the model. A variableXL not included in a causal DAG may eventually affect just one variable Xv, with v∈V; otherwise if several variables are affected then a more general class of models is needed, called semi-Markovian networks (not considered further in this work).For an updated presentation of the approach see [18] while [19] warns against the blind definition of DAGs within the causal semantic in observational studies. He also reconsiders the foundations of causal inference by anchoring them to the extended conditional independence (C.I.), both at algebraic level and with a graphical counterpart based on influence diagrams.


## 2.4. Inference about the Structure of Bayesian Networks

In all cases in which strong prior information is absent, structural learning is performed by means of a database𝒟=(x1,…,xn) of n exchangeable observations of the random vector XV.Several algorithms have been proposed to infer causal and probabilistic relations, but in Bayesian inference key quantities enter into the joint probability distribution of𝒟 and network’s unknowns given the context ξ:
(5)p(𝒟,θ,z∣ξ)=p(𝒟∣θ,z,ξ)·p(θ∣z,ξ)·p(z∣ξ),
where θ=(θv1,pa(v1),…,θvK,pa(vK)) are vectors of parameters characterizing the conditional probability distributions of Xvi given Xpa(vi); variable Z indicates the unknown DAG, and it is built as a bijection from the set of DAGs (fixed V) to a subset ΩZ of natural numbers.The likelihood functionp(𝒟∣θ,z,ξ) is typically expressed as a product of multinomials by using sufficient statistics. Whenever expert beliefs are reasonably captured by Dirichlet prior distributions for elements of θ, and they are elicited as marginally independent, closed-form integration marginalizes out thetas. The resulting marginal distribution p(𝒟∣z,ξ) has a reduced dimensionality and may be optimized with respect to z while looking for optimal structures characterized by the highest posterior probability values [9]:
(6)p(𝒟,z∣ξ)=∫‍p(𝒟∣θ,z,ξ)·p(θ∣z,ξ)·p(z∣ξ)·dθ=p(z∣ξ)·∏i=1K‍∏j=1qi‍Γ(αi,j)Γ(αi,j+Ni,j)·∏s=1ri‍Γ(αi,j,s+Ni,j,s)Γ(αi,j,s)
with ri being the number of states taken by Xi and qi of Xpa(vi); sufficient statistics are Ni,j,s for the ith variable taking the sth state while its parent configuration is in the jth state; αi,j=∑s=1ri‍αi,j,s and Ni,j=∑s=1ri‍Ni,j,s. A clever choice of hyperparameters guarantees the likelihood equivalence; that is, all BNs equivalent as regards the set of encoded conditional independence relationships are equally scored by (6).Equation (6) shows that the elicitation of the prior distribution p(z∣ξ) is a key step to perform Bayesian structural learning using 𝒟.


## 2.5. Plausible Network Features

A candidate structurez∈ΩZ on a fixed set of nodes V is plausible if it has got structural features (SFs) believed to be relevant by the expert. Following Stefanini [11, Definition  1], we recall the formal definition.Definition 1 (SF).
A structural feature (SF)ℛj(z,w) in a reference set ℛ for the set of DAGs on V is a predicate describing a plausible probabilistic or causal characteristic of the unknown directed acyclic graph z∈ΩZ. Argument w is in the partition 𝒲 of a given numeric domain ΩW of variable W. An atomic structural feature (ASF) ℛj(z) does not depend on any auxiliary variable W.A reference setℛ is a collection ℛ={ℛj:j∈J} of SFs indexed in a set J, with nf being the number of considered SFs. A proposition might be defined to carry disbelief to a candidate structure, but the feature-rises-belief direction is conveniently adopted here.An example makes the previous definition operational. Let us defineℛ1(z,w)= “The number of immediate causes of Xv3 is in w” to capture the expert belief about the number of parents of node v3 in the unknown DAG. An expert may consider ΩW={0,1,…,12} and 𝒲=({0},{1,2,3},{4,5},{6,…,12}). Given a candidate structure z, its plausibility will be determined by the element w∈𝒲 that makes the predicate true. A simple representation of configurations taken by reference features is obtained by descriptors [11, Definition  2].Definition 2 (descriptors).
A descriptorRi for the SF ℛi is a map:
(7)Ri:𝒲i×{false,true}⟶{0,1,2,…,|𝒲i|}
so that (w,false)↦0 for all w∈𝒲i and (w,true)↦hw for all w∈𝒲i; that is, a different integer is associated with each w if true. The vector R=(R1,R2,…,Rnf) defined on the Cartesian product ΩR=⊗i=1nfΩRi=⊗i=1nf{0,1,…,|𝒲i|} is called vector of descriptors. The descriptor of an ASF is defined by false↦0 and true↦1.Thejth configuration of descriptors in R is indicated as rj=(r1,j,r2,j,…,rnf,j)∈ΩR while a generic configuration is indicated as r∈ΩR. The jth configuration of a subvector of R defined by indexes in A is rA,j; for example, r{1,3},j=(r1,j,r3,j).VectorR induces equivalence classes on ΩZ; that is, 𝒵={𝒵r:r∈ΩR} contains sets of structures, with 𝒵r made by all those DAGs sharing the same configuration r.Note that members of the same equivalence class must be associated with the same degree of belief by the principle of insufficient reason: they differ in irrelevant and unconsidered ways by construction.Despite the generality of Definition1 some features are expected to occur more often than others. Below, some of them are described without pretending to be exhaustive.


### 2.5.1. Indegree and Outdegree

In applications characterized by a small sample size, it is useful to impose a sharp constraint during the greedy search of top-scored candidate structures. The maximum number of arrows entering into a node, the indegree, is set to a small integer, say2 to 5, to exclude structures with a large number of parents from the consideration. A large number of parents implies a CPT with a huge number of parameters which are affected by large uncertainty after conditioning to observed data because of sampling zeros. A similar constraint may be set on the number of arrows leaving a node.Two reference features naturally embed this kind of information.(i)
“The maximum number of arrows reaching a node is innid for all nodes in V,” where nid is a small set of integers close to 1 elicited from the expert.(ii)
“The maximum number of arrows leaving a node is innod for all nodes in V,” where nod is a small set of integers close to 1 elicited from the expert.The above two features may be exploited to increase the plausibility of candidate structures which are sparsely connected. Higher control on connectivity is obtained by considering the fraction of nodes showing a given degree of connectivity, as described below.


### 2.5.2. Partitioned Connectivity

A different way of characterizing the connectivity is obtained by eliciting the minimum fraction of nodes in DAG showing a given number of parents (children) and by iterating the elicitation from0 parents (children) up to a small integer s.Lets be a small integer representing the maximum number of parents (children) to be considered. Let W=(W0,…,Ws) be a vector of nondecreasing numbers in [0,1], with
(8)ΩW={(ω0,ω1,ω2,…,ωb,…,ωs):0hhh≤ωb≤ωb+1<ωs=1}
being the sample space. The elicitation of this feature is based on a vector a=(a0,a1,a2,…,ab,…,as), with 0≤ab≤ab+1<as=1, and on the induced partition 𝒲={w0,w1} where
(9)w1={(ω0,ω1,…,ωs):ωb≥ab,b=0,1,…,s}
and w0=ΩW∖w1.Two reference features naturally embed the extended evaluation of connectivity:(i)
“Thea-partitioned inconnectivity of degree s is in w1,” with a, s, and w1 defined above; a candidate structure z shows this feature if the cumulative relative frequency Fb of nodes in z with number of parents equal or less than b is greater than ab; that is, Fb≥ab, b=1,2,…,s.(ii)
“Thea-partitioned outconnectivity of degree s is in w1,” with a, s, and w1 defined above; a candidate structure z shows this feature if the cumulative relative frequency Fb of nodes in z with number of children equal or less than b is greater than ab; that is, Fb≥ab, b=1,2,…,s.Trained experts could prefer a conventional total number of nodes equal to100 to elicit cumulative percentages instead of cumulative fractions of nodes.In Table1, an example is shown where s=5 and the elicited vector a is defined on fractions. A candidate structure has the partitioned inconnectivity feature if the fraction of root nodes is equal or above 0.01, while the proportion of nodes with at least 1 parent is equal or above 0.2 and so on.Table 1
Elicited vectora, with s=5. In the top row, the number b of parents is reported, while in the second row the correspondent minimum fraction of nodes ab is shown.
Number of parents
0
1
2
3
4
5
Cumulative fraction of nodes
0.01
0.2
0.4
0.6
0.8
1.0


### 2.5.3. Direct Cause and Direct Effect

The reference feature “The variablexvi is an immediate cause of variable xvj” refers to xvi as a parent of xvj so that by setting (intervening on) the value of the variable xvi to a given value, the distribution of Xvj is modified. This relation holds at the level of selected granularity; thus it may change if the collection of variables (nodes in V) is modified.


### 2.5.4. Causal Ancestors

The “direct cause” feature may be extended by considering a variablexvi which is on a causal (directed) path reaching node xvj. In this case the reference feature is “The variable xvi is an indirect cause of variable xvj”; thus the expert believes that one or more variables mediate the effect of xvi on xvj.


### 2.5.5. Causal Hubs

A hub node in a network is characterized by a high number of arrows leaving it. A reference feature which captures the local connectivity of nodevi is “Node vi is a hub node of at least outdegree w,” with the outdegree indicating the number of arrows originated in vi and w a set of integers. Note that the defined feature is a localized version of the outdegree feature.The expert might believe that a hub node should be present, but without indicating a specific node. In this case thea-partitioned outconnectivity feature (with a large s) can be exploited for this purpose.


### 2.5.6. Conditional Independence Relationships

A statement about C.I. among three disjoint subsets of random variables may take the following form: “The random vectorXA is conditionally independent from the random vector XB given vector XS,” with A,B,S being disjoint subsets of nodes in V.


## 2.5.1. Indegree and Outdegree

In applications characterized by a small sample size, it is useful to impose a sharp constraint during the greedy search of top-scored candidate structures. The maximum number of arrows entering into a node, the indegree, is set to a small integer, say2 to 5, to exclude structures with a large number of parents from the consideration. A large number of parents implies a CPT with a huge number of parameters which are affected by large uncertainty after conditioning to observed data because of sampling zeros. A similar constraint may be set on the number of arrows leaving a node.Two reference features naturally embed this kind of information.(i)
“The maximum number of arrows reaching a node is innid for all nodes in V,” where nid is a small set of integers close to 1 elicited from the expert.(ii)
“The maximum number of arrows leaving a node is innod for all nodes in V,” where nod is a small set of integers close to 1 elicited from the expert.The above two features may be exploited to increase the plausibility of candidate structures which are sparsely connected. Higher control on connectivity is obtained by considering the fraction of nodes showing a given degree of connectivity, as described below.


## 2.5.2. Partitioned Connectivity

A different way of characterizing the connectivity is obtained by eliciting the minimum fraction of nodes in DAG showing a given number of parents (children) and by iterating the elicitation from0 parents (children) up to a small integer s.Lets be a small integer representing the maximum number of parents (children) to be considered. Let W=(W0,…,Ws) be a vector of nondecreasing numbers in [0,1], with
(8)ΩW={(ω0,ω1,ω2,…,ωb,…,ωs):0hhh≤ωb≤ωb+1<ωs=1}
being the sample space. The elicitation of this feature is based on a vector a=(a0,a1,a2,…,ab,…,as), with 0≤ab≤ab+1<as=1, and on the induced partition 𝒲={w0,w1} where
(9)w1={(ω0,ω1,…,ωs):ωb≥ab,b=0,1,…,s}
and w0=ΩW∖w1.Two reference features naturally embed the extended evaluation of connectivity:(i)
“Thea-partitioned inconnectivity of degree s is in w1,” with a, s, and w1 defined above; a candidate structure z shows this feature if the cumulative relative frequency Fb of nodes in z with number of parents equal or less than b is greater than ab; that is, Fb≥ab, b=1,2,…,s.(ii)
“Thea-partitioned outconnectivity of degree s is in w1,” with a, s, and w1 defined above; a candidate structure z shows this feature if the cumulative relative frequency Fb of nodes in z with number of children equal or less than b is greater than ab; that is, Fb≥ab, b=1,2,…,s.Trained experts could prefer a conventional total number of nodes equal to100 to elicit cumulative percentages instead of cumulative fractions of nodes.In Table1, an example is shown where s=5 and the elicited vector a is defined on fractions. A candidate structure has the partitioned inconnectivity feature if the fraction of root nodes is equal or above 0.01, while the proportion of nodes with at least 1 parent is equal or above 0.2 and so on.Table 1
Elicited vectora, with s=5. In the top row, the number b of parents is reported, while in the second row the correspondent minimum fraction of nodes ab is shown.
Number of parents
0
1
2
3
4
5
Cumulative fraction of nodes
0.01
0.2
0.4
0.6
0.8
1.0


## 2.5.3. Direct Cause and Direct Effect

The reference feature “The variablexvi is an immediate cause of variable xvj” refers to xvi as a parent of xvj so that by setting (intervening on) the value of the variable xvi to a given value, the distribution of Xvj is modified. This relation holds at the level of selected granularity; thus it may change if the collection of variables (nodes in V) is modified.


## 2.5.4. Causal Ancestors

The “direct cause” feature may be extended by considering a variablexvi which is on a causal (directed) path reaching node xvj. In this case the reference feature is “The variable xvi is an indirect cause of variable xvj”; thus the expert believes that one or more variables mediate the effect of xvi on xvj.


## 2.5.5. Causal Hubs

A hub node in a network is characterized by a high number of arrows leaving it. A reference feature which captures the local connectivity of nodevi is “Node vi is a hub node of at least outdegree w,” with the outdegree indicating the number of arrows originated in vi and w a set of integers. Note that the defined feature is a localized version of the outdegree feature.The expert might believe that a hub node should be present, but without indicating a specific node. In this case thea-partitioned outconnectivity feature (with a large s) can be exploited for this purpose.


## 2.5.6. Conditional Independence Relationships

A statement about C.I. among three disjoint subsets of random variables may take the following form: “The random vectorXA is conditionally independent from the random vector XB given vector XS,” with A,B,S being disjoint subsets of nodes in V.


## 2.6. The Degree of Belief

The prior distribution on the space of structures onV is obtained by “extending the argument;” that is,
(10)P[Z=z∣ξ]=∑r∈ΩR‍P[Z=z∣R=r,ξ]·P[R=r∣ξ].
By recognizing that R induces the partition 𝒵, it follows that
(11)p(z∣ξ)=1nr[z]·P[R=r[z]∣ξ]=p(r1[z],r2[z],…,rnf[z]∣ξ)nr[z],
where nr[z] is the cardinality of the equivalence class in which z is located and where the structural configuration of DAG z is r[z]. In [12] the size of each equivalence class 𝒵r was estimated by Monte Carlo simulation to face the combinatorial explosion in the number of DAGs to be enumerated with the increase of the number of nodes:
(12)|𝒵^r|=NV(Nr+1)NT+1,
where NT is the total number of DAGs uniformly sampled from the space of all DAGs on V, NV is the size of such space (see [10]), and Nr≤NT is number of sampled DAGs showing configuration r.The numerator on the right of (11) represents the joint belief on each configuration of descriptors. While the elicitation in full generality becomes cognitively and numerically overwhelming around 7 descriptors on, some parsimony is achieved if a small number of descriptors may be considered at one time, so that conditional independence relationships among descriptors may be exploited. This is a choice available to the expert through the definition of an order relation on descriptors.Definition 3 (ordered partition).
The ordered partition(13)𝒪=(𝒪1,𝒪2,…)
of descriptors is defined by the expert to indicate disjoints subsets of SFs to be jointly considered during the elicitation, from the first subset 𝒪1 to the last.Otherwise stated, the expert decomposes the whole elicitation problem following an order taken from the substantive content of the specific problem domain: features are grouped according to the priority in the elicitation.The elicitation of the ordered partition is performed by formulating questions in the language typical of a given problem domain. It is indeed difficult to define those questions in general terms, because they result to be quite abstract and they are likely to be obscure for the domain expert (see [11]). We assume here that questions are properly phrased and that an ordered partition is defined.If a strict order relation is elicited, each element of𝒪 contains one descriptor: this is a special case addressed in [20]. Another special case is represented by a trivial partition of just one subset that contains all descriptors [21]. In the two special cases above it is possible to define, respectively, a Bayesian network and a Markov network on descriptors. The general case made by several subsets, each one of cardinality two or more, was addressed by using CG models in [11] for ASFs. Below details are provided about a reference set of nonatomic features.The joint probability distribution of descriptorsp(r1,r2,…,rnf∣𝒪,ξ) is assumed to be Markov with respect to the elicited partition where descriptors within group 𝒪j define the chain component τj of the CG model (see (3)). The elicitation keeps on by asking in the language of the domain expert which descriptors in subset 𝒪j should be jointly considered while defining the degree of belief. Nodes corresponding to related descriptors are joined by undirected edges and the collection of cliques defining the chain component are found. The elicitation is iterated for all elements in the ordered partition 𝒪.The resulting CG may support the elicitation if model parameters are cognitively suited for the quantitative step, that is, if they are interpretable and easy to assess for the expert, at least after some training.The first chain component is elicited as a marginal distribution which is not conditioned on other descriptors. An undirected graphical model on descriptors in𝒪1 is defined through the multiplicative model in (4), under empty conditioning. Nevertheless some care is needed because potential functions on cliques are not uniquely defined. Here a log-linear parameterization is suggested, following [13]. For example, if the first chain component has just one clique made by three descriptors, say R1,R2,R3, we define the multiplicative model as below, after exponentiation and rearrangement of log-linear terms:
(14)p(r1,r2,r3∣ξ)p(0,0,0∣ξ)=ϕ1(r1)ϕ2(r2)ϕ3(r3)ϕ1,2(r1r2)×ϕ1,3(r1r3)ϕ2,3(r2r3)ϕ1,2,3(r1r2r3).
Thus the odds value with respect to the no-feature configuration r1,2,3=(0,0,0) is explained by a multiplicative model where each factor, for example, ϕ2,3(r2r3), is equal to one if one or more descriptors are null, otherwise they are positive (the so-called treatment parameterization).The elicitation in (14) is performed on the odds scale by asking the domain expert how many times the configuration (r1,r2,r3) is more plausible than (0,0,0) for descriptors R1,R2,R3. This question is iterated for all configurations in the Cartesian product ΩR1×ΩR2×ΩR3 after exclusion of (0,0,0). Questions are posed from single main effects towards higher order interactions terms, so that one factor at a time is considered (see algorithm below).The procedure is iterated for all cliques in the first CG componentτ1, and indeed factors already elicited in previously considered cliques are not reconsidered anymore. For example, the chain component R1—R2—R3 is made by two cliques, R1,R2 and R2,R3; thus the shared factor ϕ2(r2) is elicited just one time.The general algorithm for the first chain component is summarized below.(1)
Consider the undirected graph on descriptors elicited as the first chain component. Control questions include the following: “Are all the relevant features included in the elicitation?”, “Are all pairs of features jointly affecting the probability of a structure linked by an undirected edge?”. If needed, revise the order relation and (or) links within chain component.(2)
Check out that each descriptorRvi given its neighbors Rne(vi) is independent of all other descriptors in the first CG component.(3)
Find the cliques of such UG (model generator) [13].(4)
For each clique, elicit parameter values by using odds.
(i)
Elicit main effects{ϕi(ri,j)}, one at a time for all configurations, by assigning an odds value with respect to the baseline with no features; that is,
(15)ϕi(ri,j)=p(0,…,0,ri,j,0,…,0∣ξ)p(0,…,0∣ξ).A control question for this step is: “How much above one is the odds value for the sole presence of featureRi=ri,j>0?”(ii)
Elicit first order interactions{ϕi,l(ri,j′,rl,j′′)}, one at a time, by assigning a multiplicative term for each configuration under the question: “Which is the value of the multiplicative term {ϕi,l(ri,j′,rl,j′′)} needed to account for the interaction of feature Ri=ri,j′>0 with feature Rl=rl,j′′>0?” The expression helping in this step is
(16)p(0,…,ri,j′,0,…,rl,j′′,0,…,0∣ξ)p(0,…,0∣ξ)=ϕi(ri,j′)ϕl(rl,j′′)ϕi,l(ri,j′,rl,j′′),whereϕi(ri,j′), ϕl(rl,j′′) are already elicited; this is the cross product ratio of features Ri,Rl, after dividing by ϕi(ri,j′), ϕl(rl,j′′). It is clear that if the interaction is absent ϕi,l(ri,j′,rl,j′′)=1, while ϕi,l(ri,j′,rl,j′′)∈(0,1) means that the interaction reduces the plausibility and ϕi,l(ri,j′,rl,j′′)>1 raises the plausibility of the considered configuration.(iii)
Iterate the step above with higher order interaction terms (for all configurations) with two constraints: before moving to higher interaction terms all the terms of the same degree must have been already elicited; moreover the maximum degree of interaction among a subset of features is defined by the size of the clique they belong to (model generator). For example, with the interaction of order two{ϕi,l,k(ri,j′,rl,j′′,rk,j′′′)}, the control question might be: “Having already elicited values of ϕi(ri,j′), ϕl(rl,j′′), ϕk(rk,j′′′), ϕi,l(ri,j′,rl,j′′), ϕi,k(ri,j′,rk,j′′′), and ϕl,k(rl,j′′,rk,j′′′), which is the value of the multiplicative term {ϕi,l,k(ri,j′,rl,j′′,rk,j′′′)} needed to adjust the odds value after considering the interaction among features whose configurations are ri,j′>0, rl,j′′>0, rk,j′′′>0?”; the helper expression is
(17)p(0,…,ri,j′,…,rl,j′′,…,rk,j′′′,…,0∣ξ)p(0,…,0∣ξ)=Qϕi,l,k(ri,j′,rl,j′′,rk,j′′′),where factor(18)Q=ϕi(ri,j′)ϕl(rl,j′′)ϕk(rk,j′′′)ϕi,l(ri,j′,rl,j′′)×ϕi,k(ri,j′,rk,j′′′)ϕl,k(rl,j′′,rk,j′′′)is already elicited. Similar expressions follow straightforwardly for higher order interactions.(5)
Calculate the baseline probabilityp(0,…,0∣ξ) of a structure without features by exploiting ∑r∈ΩR‍p(r∣ξ)=1.(6)
Revise the elicited values following the guidelines of Section2.7.(7)
Move to the next chain component.The algorithm presented above may be also applied to chain componentsτ2,τ3,…, that is, with a nonempty set of parents. The elicitation of parameters in a conditional Markov network is performed by iterating the algorithm for the first chain component over each configuration of conditioning parents. It follows that the elicitation burden depends on the number of parents, more precisely, on the cardinality of the Cartesian product where factors are samples spaces of parents.The algorithms defined above lead to a prior distribution on configurations, but the elicitation does not end before the revision of elicited values takes place.


## 2.7. Implementation and Revision Issues

The flexibility achieved by defining predicates entails potential pitfalls that should be considered during an elicitation run with structural features.The revision of elicited beliefs is always needed because the expert is not expected to provide unbiased elicited values, especially with limited training or in complex problem domains. Revision and elaboration of the prior distribution [22, and references therein] are applied to control elicitation bias and other causes of poor elicitation.Overelicitation is an important practice to check for the presence of elicited quantities that do not reflect the expert degree of belief. Nevertheless, in large networks elaboration of elicited probability values is the main resource for checking the quality of elicitation. In the proposed framework, the elicitation through reference features produces proper probability distributions; therefore one elaboration consists of inspecting one or more margins of the probability distributionp(r1,…,rnf∣ξ), for example, the bivariate margin p(ri,rj∣ξ), to look for configurations whose plausibility causes surprise or disbelief in the expert. This operation is particularly meaningful if inspected margins do not straightforwardly relate to elicited odds, for example, taking descriptors belonging to different cliques. Surprise and disbelief ask for the revision of elicited values. If an association between selected margins and bias is suspected, random selection of margins is an option.The stability of elicited values against different order relations on descriptors (reference features) should not be assumed. Nevertheless, in complex problem domains overelicitation made by the repetition of the interview with other ordered partitions not only seems unacceptable as regards the work load, but it could even be cognitively unfeasible, for example, if the original ordered partition selected by the expert is induced by a scientific hypothesis.The core of the proposed approach is based on predicates representing structural features. The expert might believe that some configurations of features are plausible although they are incompatible with DAGs, for example, because they imply the presence of cycles. A positive probability value would be assigned to such a configurationr*, but this is not an instance of elicitation bias if the elicited distribution properly matches expert beliefs. If this is the case, p(z∣ξ)=0 because P[Z=z∣R=r*]=0.The way a predicate is specified determines the granularity of the elicitation and the cardinality of equivalence classes in𝒵. For example, let us consider two nodes vi∈V and vj∈V and the reference feature ℛ1= “Nodes vi, vj are not descendent of other nodes in VR”.Reworking the original proposition we have two simpler predicates:(i)
ℛ1′= “Node vi is not a descendent of other nodes in VR”;(ii)
ℛ1′′= “Node vj is not a descendent of other nodes in VR.”and they can be considered by conjunction. In Table 2, the relation between the original reference feature (right) and the conjoint components (left) is shown: ¬ℛ1 collects three configurations generated by the conjunction of simpler predicates. It follows that simpler predicates are needed if the expert degree of belief changes over collapsed configurations in ℛ1, that is, ¬ℛ1′∩ℛ1′′, ℛ1′∩¬ℛ1′′, ¬ℛ1′∩¬ℛ1′′. While nothing prevents the expert from defining rich predicates, care should be taken to select a granularity suited to properly represent expert beliefs.Table 2
Relation between the original feature (right) and the conjunction of two simpler structural features (left).
Simpler feature1
Simpler feature2
Original feature
ℛ
1
′
ℛ
1
′
′
ℛ
1
¬
ℛ
1
′
ℛ
1
′
′
¬
ℛ
1
ℛ
1
′
¬
ℛ
1
′
′
¬
ℛ
1
¬
ℛ
1
′
¬
ℛ
1
′
′
¬
ℛ
1There are indeed several different ways of formulating a predicate. If two reference sets of features induce the same set of equivalence classes𝒵 then they are operationally equivalent. Nevertheless, from a cognitive standpoint, substantial differences in the ease of elicitation might depend on the way propositions are formulated [22]. Let us consider the reference feature ℛ1= “Nodes vi,vj precede all other nodes.” Does the expert use “precede” as “come before” all other nodes in the order relationship of nodes in V? Given a DAG with vj disconnected from other nodes how to answer? Does the expert use “precede” meaning “are ancestors of all other nodes in V” or to indicate that those nodes do not receive arrows coming from other nodes? This example makes clear that some training is mandatory in order to make an expert effective in the elicitation: a trained expert is expected to choose the right granularity in the elicitation and to define meaningful predicates, that is, statements straightforwardly true/false when applied to any DAG defined on V.Several reference features are jointly considered in actual applications, a number typically far beyond what the expert may simultaneously consider with success. Thus there is the possibility that exclusion and implication relations are not recognized. For example, let us consider two features:ℛ3= “The indegree is three or less” and ℛ25= “vi is a sink node.” After reworking the last feature, the expert reformulates the statement as “Node vi has indegree ten or more.” Clearly the plausibility of ℛ3∧ℛ25 should be null; otherwise the actual interpretation of ℛ3 is “The indegree of all but vi is three or less.”Implication must be recognized to maintain the probabilistic coherence; for example, letℛ3= “Variable Xvi is a direct cause of Xvj” and ℛ5= “Variable Xvi is an ancestral cause of Xvj” be two reference features. Clearly ℛ3 implies ℛ5 and the joint plausibility of both features is bound to the plausibility of ℛ3. In a normalized reference set logical relationships among features are properly handled.


## 3. Results and Discussion

A few seminal approaches to the elicitation of beliefs on structures are reconsidered as special cases in the proposed framework based on structural features. A published case study on breast cancer (BC case-study) [23] will be (partially) exploited at illustrative purposes. The whole set of nodes V includes: age (AGE), the proliferative index-marker Ki67/MIB-1 (PROLN), oestrogen receptors (ER), progesterone receptors (PR), the receptor tyrosine kinase HER2/neu (NEU), and the P53 protein.


### 3.1. Buntine 1991

In the seminal paper of Buntine [24], a prior distribution on the set of DAGs for a fixed set V is defined by assuming a total ordering of nodes in the context ξ. The probability of the parent set for node v is defined by the product of probability for events like “There is an edge y→v,” shortly P["y→v"], extended to each y preceding v in the order relation. The subjective probability value elicited for a network structure z is calculated by marginal independence of parent sets. The original formalization defines a total ordering ≺, so that if y≺xi it may belong to the parent’s set Πi of xi; then a full specification of beliefs for each edge in the directed graph is needed and measured in units of subjective probability; finally the independence of parent sets (Π1,…,Πn) is assumed. The distribution on the set of structures is ([24], modified)
(19)p(z∣≺,ξ)=p(Π1,…,Πn∣≺,ξ)=∏i=1n‍p(Πi∣≺,ξ)p(Πi∣≺,ξ)=∏y∈Πi∧y≺x‍P("y→xi"∣≺,ξ)·(20)∏¬{y∈Πi}∧y≺x‍[1-P("y→xi"∣≺,ξ)].
Given the order relation (v3,v4,v1,v2) on four nodes, n=4, parent sets are Π3=∅, Π4⊆{yv3}, Π1⊆{yv3,yv4}, Π2⊆{yv3,yv4,yv1}.Despite the huge importance of Buntine’s seminal work [24], some limitations should be underlined. Under the causal semantic of BNs, the expert might fail in stating such node order which defines the “causal flow” along nodes. In large regulatory networks of system biology a lot of assignments are expected to be 0.5 because expert beliefs may involve a small subset of arrows. We also expect that some plausibility assignments depend on what is already assigned, for example, due to biological substantive laws, but the need of such conditioning is not accounted for. Finally, we remark that several node orders are compatible with the same sparse DAG on V; therefore the specification of a strict order should not be enforced.The above example may be straightforwardly cast in terms of reference features. The order relation is part of the contextξ, and we define the set of reference features Ri,j= “An arrow is from vi to vj,” with all possible pairs (vi,vj) in which vi precedes vj in the total ordering ≺ equal to (v3,v4,v1,v2). The reference set is ℛ={ℛ3,4,ℛ3,1,ℛ3,2,ℛ4,1,ℛ4,2,ℛ1,2}. In this case a trivial Bayesian network without arrows is equivalent to the prior distribution in (19) and (20) if conditional probability tables are defined by the same values specified in (20) for each feature; that is, features are marginally independent. To see this, note that the context puts a sharp constraint on the space of structures and that the cardinality of the equivalence classes for each configuration r of reference features is equal to one; that is,
(21)p(z∣ξ)=∏i=1nf‍P[Ri=ri[z]∣ξ].The formal approach proposed in this paper allows much more flexibility, for example, by restricting the set of nodes to be ordered and by introducing dependence among arrows. For example, letVB=(v3,v4,v1,v2)⊂V be an ordered subset of nodes taken from BC case study. The reference feature ℛo= “VB is the order on the relevant subset of nodes” induces two equivalence classes, the first made by DAGs on V which do not satisfy ℛo, the second one is made by structures without arrows violating the left-to-right order defined in VB. Besides the sharp constraint obtained by setting P[Ro=1∣ξ]=1, the expert might consider such feature uncertain, thus preferring a degree of belief in the set (0,1). Further features could be defined to build the reference set:
(22)ℛ=ℛo∪{ℛ3,4,ℛ3,1,ℛ3,2,ℛ4,1,ℛ4,2,ℛ1,2}
and the ordered partition
(23)𝒪=({Ro},{R3,4,…,R1,2})
captures weaker causal relationships by features like ℛi,j= “Variable Xvi is a causal ancestor of Xvj.” A CG model made by two components makes possible to elicit conditional beliefs about v4→v2 given the presence of an arrow v3→v4 and given the lack of v3→v1, without assuming a strict order relation on nodes in V.


### 3.2. Heckerman et al. (1995)

In Heckerman et al. [9], a prior network z𝒫 was elicited and compared to a candidate network z by counting the number of different edges, δ, with a high degree of belief assigned to structures closely resembling the prior network. The authors suggested to elicit the hyperparameter 0<k<1 and to define the prior distribution to be proportional to kδ.Among the limitations penalizing the use of this prior we found the following:(i)
The impossibility of specifying the degree of belief if it depends not only on the number of different edgesδ but also on their position and type; the presence/absence/direction of an arrow may have an impact on the belief about other edges.(ii)
The elicitation about a subset ofV is not addressed.(iii)
The causal semantic is natural for this approach because each arrow represents an immediate cause; it seems difficult to mix probabilistic and causal beliefs by counting differences in arrows because a DAG in the probabilistic semantic is just a member of an equivalence class of DAGs representing the same collection of conditional independence relationships.A simple reformulation is oriented to computation and it involves three operators: arrow deletion, insertion, and change of direction. The reference set of atomic features isℛ={ℛj(z):j=0,1,2,…,J}∪{Ra(z):a=J+1} with ℛj(z)= “The application of j operations produces z𝒫,” and where Ra(z)= “The application of a or more operations produces z𝒫.” An undirected graph made by just one clique is associated with potentials represented in Table 3, where J=3. On the right of Table 3, values of the potential function are shown. The normalization constant is σ=(θ0+θ1+θ2+θ3+θa)-1. It is obviously possible to make the two approaches as close as desired by setting θi∝ki with i≤3.Table 3
Potential function for the Heckerman et al. [9] prior.
R
0
R
1
R
2
R
3
R
a
Potential
0
0
0
1
0
θ
3
0
0
1
0
0
θ
2
0
1
0
0
0
θ
1
1
0
0
0
0
θ
0
0
0
0
0
1
θ
a
Otherwise
0An even simpler reformulation exploits the incompatibility of the above features and it is based on just one reference feature:ℛ1(z,a)= “The application of a operations produces z𝒫,” with a∈𝒜1=({0},{1},…,{4,5,…}), with arrow manipulations (insert-delete-change) defined as above. Feature ℛ1 essentially defines a plausible neighborhood with respect to a prior network. Further reference features may be introduced to refine such plausibility, for example, by also considering one causal, ℛ2, and one C.I., ℛ3, relationships. In other words, a candidate network “close” to the prior network could be associated with a high prior probability which is then tuned according to the presence/absence of two other relevant features. A natural ordered partition on descriptors could be 𝒪=({R1},{R2,R3}); thus a two-component chain graph model may support the elicitation.A different kind of reformulation is based on the full-probabilistic semantic in which a prior DAGz𝒫 is just a way to define a collection of C.I. relationships. In this case, it is natural to define a structural feature for each conditional independence relationship if it is relevant according to the expert among those represented by z𝒫. The reference set ℛ={ℛ1,…,ℛnf} in this case is a collection of plausible C.I. statements taken from the prior DAG. Note that, although it is not essential to draw such prior DAG, it may be useful because a general collection of C.I. statements does not necessarily imply the existence of a compatible DAG.If none among the C.I. relations is preeminent the ordered partition of descriptors contains just one element, say𝒪=({R1,…,Rnf}), and a CG model made by one component (undirected graphical model) is suited for the elicitation.


### 3.3. Imoto et al. (2003)

In the seminal paper of Imoto et al. [25], the authors developed a framework for combining microarray data and biological knowledge while learning the structure of a BN representing relationships among genes. The proposed criterion has two components, and the second one is particularly interesting because it captures the a priori biological knowledge.Following their original notation with minor modifications,π(G) is the prior distribution of network G. Then, the interaction energy Ui,j of the edge from (gene) vi to (gene) vj is defined on a sample space which is categorized into I values, say H1,H2,…,HI. For example if gene vi regulates gene vj then Ui,j=H1>0, but if not much is known about such potential regulation then Ui,j=H2>H1. The total energy of a network G is E(G)=∑(vi,vj)∈G‍Ui,j; thus the sum is taken over existing edges in G. The total energy may be rewritten by collecting the parents of each node; thus E(G)=∑vj∈V‍∑vi∈pa(vj)‍Ui,j=∑vj∈V‍Ej. The (prior) probability of network G is modeled by the Gibbs distribution:
(24)π(G)=σ-1exp(-ζE(G)),
where ζ>0 is a hyperparameter and σ is the normalizing constant, also called the partition function σ=∑G∈𝒢V‍exp(-ζE(G)), with 𝒢V being the collection of all DAGs on a fixed set of nodes V.Operationally, the prior information is coded into a square matrixU of size defined by the number of genes, with each ui,j corresponding to ζH1 or ζH2 according to the prior belief. Beliefs in protein-protein interactions are coded by ui,j=uj,i=ζH1. Protein-DNA interactions between the transcription regulator vi and the controlled gene vj are accounted by setting ui,j=ζH1 and uj,i=ζH2. Some genes are controlled by a transcription regulator through a consensus motif in their DNA promoter region. If genes vj1,vj2,…,vjn have the consensus motif and they are regulated by gene vi then uj1,i=⋯=ujn,i=ζH2 and ui,j1=⋯=ui,jn=ζH1.The seminal approach of Imoto et al. [25] suffers of two main limitations. They stated that the biological knowledge should suggest the partitioning of the underlining continuous energy function but it is not clear how, even after invoking the metaphor of energy from physics. In Example 3.2, they tried ζH1=0.5 (but also ζH1=1) and optimized the selection of ζH2=2.5, a procedure not much in line with pure preexperimental Bayesian elicitation. Moreover, the sum in the partition function is taken on the set of DAGs on V, and it becomes quite intractable from 5 nodes on. It follows that their approach is substantially a way to build a score function in a spirit similar to [26], but without providing a flexible support for the calibration of the score function.The prior information considered by these authors can be also expressed in terms of reference features, for example, by considering features likeℛi,j= “Gene vi regulates gene vj,” for all relevant pairs of genes. If preeminent features are absent, an UG model formally resembles the expression based on energy functions but with some major differences:(1)
the parameterization is related to subjective probability through odds values;(2)
the normalization constant is easier to calculate, at least if the number of features is less than the total pair of genes (some genes omitted);(3)
the general calibration constantζ disappears, although a similar tuning has been considered elsewhere [21] to smooth raw elicited odds values while trying to compensate for the well known human tendency towards overstating odds values.


### 3.4. Werhli and Husmeier (2007)

The authors [27], building on the work of Imoto et al. [25], defined a prior information matrix B whose elements Bi,j∈[0,1], with i,j being a pair of integers for nodes vi and vj and the relation vi→vj. If no prior preference about the presence of such arrow is elicited, then Bi,j=0.5; if 0≤Bi,j<0.5 elicited beliefs put more plausibility on the lack of arrow vi→vj; if 0.5<Bi,j≤1 higher plausibility favors the presence of the arrow vi→vj. Note that elements Bi,j in B are not probabilities.The calculation of the prior probability for a candidate DAGG is straightforward if it is represented through an adjacency matrix in which the element Gi,j in row i and column j is 1 if the DAG has an arrow form vi to vj, and it is zero otherwise. At first the energy of the DAG is calculated as
(25)ℰ(G)=∑i,j‍|Bi,j-Gi,j|
with Gi,j and Bi,j being, respectively, the elements of G and B. The probability elicited for G is
(26)π(G∣β,B)=σ-1(β,B)exp(-βℰ(G)),
where β is a hyperparameter regulating the overall strength of the elicited degree of belief and with σ being the partition function which depends on a sum of energy values over the collection of all DAGs on a fixed set of nodes V. It follows that an increase of energy is related to a larger mismatch with the prior B. In the limit of β→0 a noninformative prior is obtained, while for β diverging to infinity peaked priors are defined.The relationship with the seminal approach of Imoto et al. [25] is clear but despite the easier parameterization chosen to define the energy function, several limitations are still present. First of all, the overall strength β is not straightforwardly related to (subjective) probability; therefore at some point the expert has to play numerically with fake examples to get the feeling on reasonable values for such hyperparameter. The calibration of the whole approach is difficult because it depends on the calculation of the normalization constant which is hard to obtain due to the calculation of energy values for all DAGs on a fixed set of nodes V. Despite the shortcut proposed by the authors, the calculation of the normalization constant remains as a bottleneck of the approach.The authors increased the flexibility of their approach in representing prior beliefs by using two matricesB(1) and B(2):
(27)π(G∣β1,B(1),β2,B(2))={σ(β1,B(1),β2,B(2))}-1×exp(-{β1ℰ1(G)+β2ℰ2(G)})
with ℰ1(G), ℰ2(G) being the energy functions, respectively, depending on B(1) and B(2). In the elicitation based on a widely adopted database of metabolic pathways, values of B are defined as the ratio mi,j/Mi,j, with Mi,j being the number of times two genes appear in a pathway and mi,j the number of times that they are linked inside such pathway. Results of a study on 25 genes measured at 73 time points suggested that the procedures using prior information outperformed those without it. Nevertheless, the above mentioned limitations get even worse under such generalization; for example, the explosion in the number of terms entering the partition function must be taken under control by introducing a limitation on the number of arrows entering a node, a technical constraint which may lead to biased elicited beliefs.The expressivity obtained by these authors through the use of twoB matrices may be similarly obtained using reference features with some advantages. Back to the BC case study, a major hypothesis under which the elicitation may be performed could be ℛ1= “ER is a hub gene,” so that given the configuration R1=1 the expert has to express further conditional beliefs about the other markers; for example, ℛ2= “PR regulates P53,” ℛ3= “P53 acts on KI67,” ℛ4= “P53 acts on NEU.” In this context it is natural to consider the order relation on features 𝒪=({R1},{R2,R3,R4}); thus a chain graph model is suited for the elicitation with just node R1 in the first chain component and as many arrows leaving from R1 as needed to capture changes of conditional beliefs due to a switch in the major hypothesis, from R1=1 to R1=0. Note that, using a chain graph on features, it is possible to tune the flexibility exactly of the amount needed, without the unnecessary and blind consideration of all pairs of genes.


### 3.5. Discussion

The expressivity achieved through reference features is wide whether probabilistic or causal information is elicited. Many limitations found in other approaches depend on the consideration of arrows as the key building block of the elicitation. Further restrictions are due to the use of marginal relationships among arrows, so that severe constraints on the expressivity of the approach follow. The elicitation based on reference features has maximum resolution because, for a suitable set of reference features, it is possible to define a prior distribution characterized by a probability value for each DAG defined onV.There are useful side effects in the approach based on structural features. First, the elicitation effort does not depend on the size of the space of DAGs onV, but on the size of the collection of features. Another side effect is related to the cardinality of equivalence classes in 𝒵. Two candidate DAGs z′ and z′′ characterized by two different configurations of structural features may receive a very different prior probability value just because the cardinality of the equivalence classes they belong to is very different, say nr[z′]⋘⋙nr[z′′] (see (11)). Note that the cardinality of equivalence classes must be taken into account to preserve probabilistic coherence.The above description of the elicitation process did not deal with computational issues that are very important for applications. The proposed approach is suited to large networks if the number of DAGs within each equivalence class is available, for example, as a Monte Carlo point estimate. Monte Carlo simulation makes it possible to explore the space of DAGs and it provides evidences about the presence of features which are logically incompatible; thus it may also suggest predicates on which the expert should focus to improve the definition of reference features. Moreover, the elicitation defines a proper prior distribution, so that different MCMC algorithms could be developed to obtain the posterior distribution ofZ. New greedy search algorithms to find plausible structures in ΩZ could be investigated to exploit the elicited reference features, for example, by developing an algorithm that generates candidate DAGs belonging to different equivalence classes in 𝒵 at each iteration of the optimization.The approach to elicitation based on graphical models does not necessitate very esoteric software, although software libraries for platforms commonly adopted in scientific computing would offer the opportunity of performing extensive testing and of investigating human heuristics specifically relevant in this elicitation framework. It is well known that graphical user interfaces facilitate the elicitation, especially if experts are not much trained in Bayesian statistics, and this is a resource which is anyway almost mandatory to face the elicitation in problem domains involving a large number of structural features.


## 3.1. Buntine 1991

In the seminal paper of Buntine [24], a prior distribution on the set of DAGs for a fixed set V is defined by assuming a total ordering of nodes in the context ξ. The probability of the parent set for node v is defined by the product of probability for events like “There is an edge y→v,” shortly P["y→v"], extended to each y preceding v in the order relation. The subjective probability value elicited for a network structure z is calculated by marginal independence of parent sets. The original formalization defines a total ordering ≺, so that if y≺xi it may belong to the parent’s set Πi of xi; then a full specification of beliefs for each edge in the directed graph is needed and measured in units of subjective probability; finally the independence of parent sets (Π1,…,Πn) is assumed. The distribution on the set of structures is ([24], modified)
(19)p(z∣≺,ξ)=p(Π1,…,Πn∣≺,ξ)=∏i=1n‍p(Πi∣≺,ξ)p(Πi∣≺,ξ)=∏y∈Πi∧y≺x‍P("y→xi"∣≺,ξ)·(20)∏¬{y∈Πi}∧y≺x‍[1-P("y→xi"∣≺,ξ)].
Given the order relation (v3,v4,v1,v2) on four nodes, n=4, parent sets are Π3=∅, Π4⊆{yv3}, Π1⊆{yv3,yv4}, Π2⊆{yv3,yv4,yv1}.Despite the huge importance of Buntine’s seminal work [24], some limitations should be underlined. Under the causal semantic of BNs, the expert might fail in stating such node order which defines the “causal flow” along nodes. In large regulatory networks of system biology a lot of assignments are expected to be 0.5 because expert beliefs may involve a small subset of arrows. We also expect that some plausibility assignments depend on what is already assigned, for example, due to biological substantive laws, but the need of such conditioning is not accounted for. Finally, we remark that several node orders are compatible with the same sparse DAG on V; therefore the specification of a strict order should not be enforced.The above example may be straightforwardly cast in terms of reference features. The order relation is part of the contextξ, and we define the set of reference features Ri,j= “An arrow is from vi to vj,” with all possible pairs (vi,vj) in which vi precedes vj in the total ordering ≺ equal to (v3,v4,v1,v2). The reference set is ℛ={ℛ3,4,ℛ3,1,ℛ3,2,ℛ4,1,ℛ4,2,ℛ1,2}. In this case a trivial Bayesian network without arrows is equivalent to the prior distribution in (19) and (20) if conditional probability tables are defined by the same values specified in (20) for each feature; that is, features are marginally independent. To see this, note that the context puts a sharp constraint on the space of structures and that the cardinality of the equivalence classes for each configuration r of reference features is equal to one; that is,
(21)p(z∣ξ)=∏i=1nf‍P[Ri=ri[z]∣ξ].The formal approach proposed in this paper allows much more flexibility, for example, by restricting the set of nodes to be ordered and by introducing dependence among arrows. For example, letVB=(v3,v4,v1,v2)⊂V be an ordered subset of nodes taken from BC case study. The reference feature ℛo= “VB is the order on the relevant subset of nodes” induces two equivalence classes, the first made by DAGs on V which do not satisfy ℛo, the second one is made by structures without arrows violating the left-to-right order defined in VB. Besides the sharp constraint obtained by setting P[Ro=1∣ξ]=1, the expert might consider such feature uncertain, thus preferring a degree of belief in the set (0,1). Further features could be defined to build the reference set:
(22)ℛ=ℛo∪{ℛ3,4,ℛ3,1,ℛ3,2,ℛ4,1,ℛ4,2,ℛ1,2}
and the ordered partition
(23)𝒪=({Ro},{R3,4,…,R1,2})
captures weaker causal relationships by features like ℛi,j= “Variable Xvi is a causal ancestor of Xvj.” A CG model made by two components makes possible to elicit conditional beliefs about v4→v2 given the presence of an arrow v3→v4 and given the lack of v3→v1, without assuming a strict order relation on nodes in V.


## 3.2. Heckerman et al. (1995)

In Heckerman et al. [9], a prior network z𝒫 was elicited and compared to a candidate network z by counting the number of different edges, δ, with a high degree of belief assigned to structures closely resembling the prior network. The authors suggested to elicit the hyperparameter 0<k<1 and to define the prior distribution to be proportional to kδ.Among the limitations penalizing the use of this prior we found the following:(i)
The impossibility of specifying the degree of belief if it depends not only on the number of different edgesδ but also on their position and type; the presence/absence/direction of an arrow may have an impact on the belief about other edges.(ii)
The elicitation about a subset ofV is not addressed.(iii)
The causal semantic is natural for this approach because each arrow represents an immediate cause; it seems difficult to mix probabilistic and causal beliefs by counting differences in arrows because a DAG in the probabilistic semantic is just a member of an equivalence class of DAGs representing the same collection of conditional independence relationships.A simple reformulation is oriented to computation and it involves three operators: arrow deletion, insertion, and change of direction. The reference set of atomic features isℛ={ℛj(z):j=0,1,2,…,J}∪{Ra(z):a=J+1} with ℛj(z)= “The application of j operations produces z𝒫,” and where Ra(z)= “The application of a or more operations produces z𝒫.” An undirected graph made by just one clique is associated with potentials represented in Table 3, where J=3. On the right of Table 3, values of the potential function are shown. The normalization constant is σ=(θ0+θ1+θ2+θ3+θa)-1. It is obviously possible to make the two approaches as close as desired by setting θi∝ki with i≤3.Table 3
Potential function for the Heckerman et al. [9] prior.
R
0
R
1
R
2
R
3
R
a
Potential
0
0
0
1
0
θ
3
0
0
1
0
0
θ
2
0
1
0
0
0
θ
1
1
0
0
0
0
θ
0
0
0
0
0
1
θ
a
Otherwise
0An even simpler reformulation exploits the incompatibility of the above features and it is based on just one reference feature:ℛ1(z,a)= “The application of a operations produces z𝒫,” with a∈𝒜1=({0},{1},…,{4,5,…}), with arrow manipulations (insert-delete-change) defined as above. Feature ℛ1 essentially defines a plausible neighborhood with respect to a prior network. Further reference features may be introduced to refine such plausibility, for example, by also considering one causal, ℛ2, and one C.I., ℛ3, relationships. In other words, a candidate network “close” to the prior network could be associated with a high prior probability which is then tuned according to the presence/absence of two other relevant features. A natural ordered partition on descriptors could be 𝒪=({R1},{R2,R3}); thus a two-component chain graph model may support the elicitation.A different kind of reformulation is based on the full-probabilistic semantic in which a prior DAGz𝒫 is just a way to define a collection of C.I. relationships. In this case, it is natural to define a structural feature for each conditional independence relationship if it is relevant according to the expert among those represented by z𝒫. The reference set ℛ={ℛ1,…,ℛnf} in this case is a collection of plausible C.I. statements taken from the prior DAG. Note that, although it is not essential to draw such prior DAG, it may be useful because a general collection of C.I. statements does not necessarily imply the existence of a compatible DAG.If none among the C.I. relations is preeminent the ordered partition of descriptors contains just one element, say𝒪=({R1,…,Rnf}), and a CG model made by one component (undirected graphical model) is suited for the elicitation.


## 3.3. Imoto et al. (2003)

In the seminal paper of Imoto et al. [25], the authors developed a framework for combining microarray data and biological knowledge while learning the structure of a BN representing relationships among genes. The proposed criterion has two components, and the second one is particularly interesting because it captures the a priori biological knowledge.Following their original notation with minor modifications,π(G) is the prior distribution of network G. Then, the interaction energy Ui,j of the edge from (gene) vi to (gene) vj is defined on a sample space which is categorized into I values, say H1,H2,…,HI. For example if gene vi regulates gene vj then Ui,j=H1>0, but if not much is known about such potential regulation then Ui,j=H2>H1. The total energy of a network G is E(G)=∑(vi,vj)∈G‍Ui,j; thus the sum is taken over existing edges in G. The total energy may be rewritten by collecting the parents of each node; thus E(G)=∑vj∈V‍∑vi∈pa(vj)‍Ui,j=∑vj∈V‍Ej. The (prior) probability of network G is modeled by the Gibbs distribution:
(24)π(G)=σ-1exp(-ζE(G)),
where ζ>0 is a hyperparameter and σ is the normalizing constant, also called the partition function σ=∑G∈𝒢V‍exp(-ζE(G)), with 𝒢V being the collection of all DAGs on a fixed set of nodes V.Operationally, the prior information is coded into a square matrixU of size defined by the number of genes, with each ui,j corresponding to ζH1 or ζH2 according to the prior belief. Beliefs in protein-protein interactions are coded by ui,j=uj,i=ζH1. Protein-DNA interactions between the transcription regulator vi and the controlled gene vj are accounted by setting ui,j=ζH1 and uj,i=ζH2. Some genes are controlled by a transcription regulator through a consensus motif in their DNA promoter region. If genes vj1,vj2,…,vjn have the consensus motif and they are regulated by gene vi then uj1,i=⋯=ujn,i=ζH2 and ui,j1=⋯=ui,jn=ζH1.The seminal approach of Imoto et al. [25] suffers of two main limitations. They stated that the biological knowledge should suggest the partitioning of the underlining continuous energy function but it is not clear how, even after invoking the metaphor of energy from physics. In Example 3.2, they tried ζH1=0.5 (but also ζH1=1) and optimized the selection of ζH2=2.5, a procedure not much in line with pure preexperimental Bayesian elicitation. Moreover, the sum in the partition function is taken on the set of DAGs on V, and it becomes quite intractable from 5 nodes on. It follows that their approach is substantially a way to build a score function in a spirit similar to [26], but without providing a flexible support for the calibration of the score function.The prior information considered by these authors can be also expressed in terms of reference features, for example, by considering features likeℛi,j= “Gene vi regulates gene vj,” for all relevant pairs of genes. If preeminent features are absent, an UG model formally resembles the expression based on energy functions but with some major differences:(1)
the parameterization is related to subjective probability through odds values;(2)
the normalization constant is easier to calculate, at least if the number of features is less than the total pair of genes (some genes omitted);(3)
the general calibration constantζ disappears, although a similar tuning has been considered elsewhere [21] to smooth raw elicited odds values while trying to compensate for the well known human tendency towards overstating odds values.


## 3.4. Werhli and Husmeier (2007)

The authors [27], building on the work of Imoto et al. [25], defined a prior information matrix B whose elements Bi,j∈[0,1], with i,j being a pair of integers for nodes vi and vj and the relation vi→vj. If no prior preference about the presence of such arrow is elicited, then Bi,j=0.5; if 0≤Bi,j<0.5 elicited beliefs put more plausibility on the lack of arrow vi→vj; if 0.5<Bi,j≤1 higher plausibility favors the presence of the arrow vi→vj. Note that elements Bi,j in B are not probabilities.The calculation of the prior probability for a candidate DAGG is straightforward if it is represented through an adjacency matrix in which the element Gi,j in row i and column j is 1 if the DAG has an arrow form vi to vj, and it is zero otherwise. At first the energy of the DAG is calculated as
(25)ℰ(G)=∑i,j‍|Bi,j-Gi,j|
with Gi,j and Bi,j being, respectively, the elements of G and B. The probability elicited for G is
(26)π(G∣β,B)=σ-1(β,B)exp(-βℰ(G)),
where β is a hyperparameter regulating the overall strength of the elicited degree of belief and with σ being the partition function which depends on a sum of energy values over the collection of all DAGs on a fixed set of nodes V. It follows that an increase of energy is related to a larger mismatch with the prior B. In the limit of β→0 a noninformative prior is obtained, while for β diverging to infinity peaked priors are defined.The relationship with the seminal approach of Imoto et al. [25] is clear but despite the easier parameterization chosen to define the energy function, several limitations are still present. First of all, the overall strength β is not straightforwardly related to (subjective) probability; therefore at some point the expert has to play numerically with fake examples to get the feeling on reasonable values for such hyperparameter. The calibration of the whole approach is difficult because it depends on the calculation of the normalization constant which is hard to obtain due to the calculation of energy values for all DAGs on a fixed set of nodes V. Despite the shortcut proposed by the authors, the calculation of the normalization constant remains as a bottleneck of the approach.The authors increased the flexibility of their approach in representing prior beliefs by using two matricesB(1) and B(2):
(27)π(G∣β1,B(1),β2,B(2))={σ(β1,B(1),β2,B(2))}-1×exp(-{β1ℰ1(G)+β2ℰ2(G)})
with ℰ1(G), ℰ2(G) being the energy functions, respectively, depending on B(1) and B(2). In the elicitation based on a widely adopted database of metabolic pathways, values of B are defined as the ratio mi,j/Mi,j, with Mi,j being the number of times two genes appear in a pathway and mi,j the number of times that they are linked inside such pathway. Results of a study on 25 genes measured at 73 time points suggested that the procedures using prior information outperformed those without it. Nevertheless, the above mentioned limitations get even worse under such generalization; for example, the explosion in the number of terms entering the partition function must be taken under control by introducing a limitation on the number of arrows entering a node, a technical constraint which may lead to biased elicited beliefs.The expressivity obtained by these authors through the use of twoB matrices may be similarly obtained using reference features with some advantages. Back to the BC case study, a major hypothesis under which the elicitation may be performed could be ℛ1= “ER is a hub gene,” so that given the configuration R1=1 the expert has to express further conditional beliefs about the other markers; for example, ℛ2= “PR regulates P53,” ℛ3= “P53 acts on KI67,” ℛ4= “P53 acts on NEU.” In this context it is natural to consider the order relation on features 𝒪=({R1},{R2,R3,R4}); thus a chain graph model is suited for the elicitation with just node R1 in the first chain component and as many arrows leaving from R1 as needed to capture changes of conditional beliefs due to a switch in the major hypothesis, from R1=1 to R1=0. Note that, using a chain graph on features, it is possible to tune the flexibility exactly of the amount needed, without the unnecessary and blind consideration of all pairs of genes.


## 3.5. Discussion

The expressivity achieved through reference features is wide whether probabilistic or causal information is elicited. Many limitations found in other approaches depend on the consideration of arrows as the key building block of the elicitation. Further restrictions are due to the use of marginal relationships among arrows, so that severe constraints on the expressivity of the approach follow. The elicitation based on reference features has maximum resolution because, for a suitable set of reference features, it is possible to define a prior distribution characterized by a probability value for each DAG defined onV.There are useful side effects in the approach based on structural features. First, the elicitation effort does not depend on the size of the space of DAGs onV, but on the size of the collection of features. Another side effect is related to the cardinality of equivalence classes in 𝒵. Two candidate DAGs z′ and z′′ characterized by two different configurations of structural features may receive a very different prior probability value just because the cardinality of the equivalence classes they belong to is very different, say nr[z′]⋘⋙nr[z′′] (see (11)). Note that the cardinality of equivalence classes must be taken into account to preserve probabilistic coherence.The above description of the elicitation process did not deal with computational issues that are very important for applications. The proposed approach is suited to large networks if the number of DAGs within each equivalence class is available, for example, as a Monte Carlo point estimate. Monte Carlo simulation makes it possible to explore the space of DAGs and it provides evidences about the presence of features which are logically incompatible; thus it may also suggest predicates on which the expert should focus to improve the definition of reference features. Moreover, the elicitation defines a proper prior distribution, so that different MCMC algorithms could be developed to obtain the posterior distribution ofZ. New greedy search algorithms to find plausible structures in ΩZ could be investigated to exploit the elicited reference features, for example, by developing an algorithm that generates candidate DAGs belonging to different equivalence classes in 𝒵 at each iteration of the optimization.The approach to elicitation based on graphical models does not necessitate very esoteric software, although software libraries for platforms commonly adopted in scientific computing would offer the opportunity of performing extensive testing and of investigating human heuristics specifically relevant in this elicitation framework. It is well known that graphical user interfaces facilitate the elicitation, especially if experts are not much trained in Bayesian statistics, and this is a resource which is anyway almost mandatory to face the elicitation in problem domains involving a large number of structural features.


## 4. Conclusions

Graphical models may be exploited to elicit beliefs about the structure of an unknown BN from experts. The joint plausibility on configurations of structural features is decomposed according to conditional independence relationships that are considered plausible by an expert. The expert may use CG models to elicit structural prior information in quite complex domains without leaving a full Bayesian framework. From the elicited CG model, and eventually by using an auxiliary Monte Carlo simulation to estimate the cardinality of equivalence classes, a (proper) subjective prior distribution on the space of DAGs is built and ready to be used with the likelihood function in order to find BN structures supported both by expert beliefs and by collected observations.No surprise that in complex domains the definition of a prior distribution may be costly. A trade-off should be found by considering the goal of the analysis, how much prior information is available, and the cost and importance of collected data. Here system biology and medicine are expected to be fields in which the proposed approach might be useful, because subjective prior information, besides providing the above mentioned benefits, also tempers the curse of dimensionality caused by structures defined on a high number of variables.



---

*Source: 749150-2014-02-05.xml*