Bayesian Networks in the Absence of Temporal Information - Bayesian Networks in R: With Applications in Systems Biology - page 49

The code above implements a setup similar to Sachs et al. ( 2005 ). We learn a net-
work structure from each bootstrap sample with a hill-climbing search and a BDe
posterior density with a very low imaginary sample size. Arcs are considered signif-
icant if they appear in at least 85 % of the networks and in the direction that appears
most frequently. Interestingly, lowering the threshold to 50 % does not change the
results of the analysis, which seems to indicate that its value is not critical in this
case.
Having computed the significance for all possible arcs, we can now build the
averaged network structure as follows:
> avg.boot = averaged.network(boot, threshold = 0.85)
As shown in Fig. 2.8 , the averaged network contains the same arcs as the network
learned from the observational data in Sachs et al. ( 2005 ). The only difference is that
since the probability of the presence of each arc and of its possible directions are
computed separately in boot.strength , we can also determine which direction
has better support from the discretized data.
As an alternative, we can average the results of several hill-climbing searches,
each starting from a different network. Such networks can be generated randomly
from a uniform distribution over the space of connected graphs with the algorithm
proposed by Ide and Cozman ( 2002 ) and implemented in random.graph .This
ensures both that no systematic bias is introduced in the learned networks and that
the search space is covered more thoroughly.
> nodes = names(dsachs)
> start = random.graph(nodes = nodes,
+ method = "ic-dag", num = 500)
> netlist = lapply(start, function(net) {
+ hc(dsachs, score = "bde", iss = 10, start = net)
+})
> rnd = custom.strength(netlist, nodes = nodes)
> rnd[(rnd$strength > 0.85) &
+ (rnd$direction >= 0.5), ]
from to strength direction
11 pmek praf 1 0.518
23 plcg PIP2 1 0.748
43 PIP3 plcg 1 0.542
44 PIP3 PIP2 1 0.790
56 p44.42 pakts473 1 0.682
57 p44.42 PKA 1 0.704
67 pakts473 PKA 1 0.522
89 PKC P38 1 0.540
90 PKC pjnk 1 0.945
100 P38 pjnk 1 0.905
> avg.start = averaged.network(rnd, threshold = 0.85)