1

I have 2 questions on bayesian network with bnlearn package in R.

library(parallel)
cl = makeCluster(4) 

set.seed(1)
b1 = boot.strength(data = learning.test, R = 5, algorithm = "hc", cluster = cl,  algorithm.args = list( score="bde", iss = 60 , restart=10, perturb = 5 ))
cl = makeCluster(4) 
set.seed(1)
b2 = boot.strength(data = learning.test, R = 5, algorithm = "hc", cluster = cl, algorithm.args = list( score="bde", iss = 60 , restart=10, perturb = 5 ))
all.equal(b1, b2)
[1] "Attributes: < Component “threshold”: Mean relative difference: 0.5 >"
[2] "Component “strength”: Mean relative difference: 1.166667"            
[3] "Component “direction”: Mean relative difference: 0.5294118"

  1. I am using the boot.strength function in bnlearn package to create bootstraps and create multiple bayesian networks to get the arc strength and direction. Each time I run, I get different results from this boot.strength function. I could not find any seed parameter for this function to get reproducible results. Even if I do set.seed(1) before the line with boot.strength function, it gives different results. Kindly help me to get reproducible results with this function. Reproducible code for this is shown above

  2. I first created one network(version1) after applying some thresholds on arc strength and direction in the output file from boot.strength function which was run without any blacklists. Then I created another network (version2) in the same manner using same thresholds for strength and direction but with some arcs blacklisted. Logically I would expect the version 2 network to have lesser no. of arcs than that in version 1 network. In other words, I have version 1 and then I apply some constraints in the form of blacklists and get version 2. Since some arcs are blacklisted in version 2, I would expect version 2 to have lesser no. of arcs. But when I repeated this experiment (having these 2 versions) multiple times, I always find that version 2 has more no. of arcs than version 1. This is very systematic. I found this issue to be faced by multiple people. Any suggestion on what could cause this anomaly would be greatly helpful. Reproducible code for this is shown below. If the below code is run mutliple times by changing thresholds for strength or changing other parameters like restart, perturb or R, you will find that in all such experiments version 2 (i.e. b2) will have more no. of arcs than version 1 (i.e. b1) systematically.

set.seed(1)
b1 = boot.strength(data = learning.test, R = 5, algorithm = "hc", algorithm.args = list( score="bde", iss = 1 , restart=10, perturb = 5 ))

badarc = expand.grid(c("A","F"), c("B","D"))
colnames(badarc) =  c("from", "to")
set.seed(1)
b2 = boot.strength(data = learning.test, R = 5, algorithm = "hc", algorithm.args = list( score="bde", iss = 1 , restart=10, perturb = 5, blacklist = badarc ))

nrow(b1[b1$strength>0.8 & b1$direction>0.5,])  # no. of arcs in version 1
nrow(b2[b2$strength>0.8 & b2$direction>0.5,])  # no. of arcs in version 2
Kumar
  • 43
  • 7
  • re 1. can you show some reproducible code please e.g. `set.seed(1); b1 = boot.strength(learning.test, algorithm = "hc") ; set.seed(1); b2 = boot.strength(learning.test, algorithm = "hc") ; all.equal(b1, b2)` returns `TRUE` for me – user20650 Sep 22 '22 at 18:59
  • "Logically I would expect the version 2 network to have lesser no. of arcs than that in version 1 network" : not necessarily, as if some edges are prohibited then several may be added to recover the associations found in the data – user20650 Sep 22 '22 at 19:03
  • @user20650, I agree with your response to my 2nd question. But at least, some times version 2 should have less # of arcs and sometimes it should have more # of arcs when this experiment is repeated multiple times. But systematically version 2 always has more # of arcs. I have updated my question above with reproducible code for my question 1. I found that cluster parameter is causing the problem since when I removed that parameter, it gives reproducible results – Kumar Sep 23 '22 at 05:57
  • Try `cl = makeCluster(4); parallel::clusterEvalQ(cl, set.seed(1)); b1 = boot.strength(...); stopCluster(cl)` – user20650 Sep 23 '22 at 13:17
  • Thanks a lot for your suggestion. clusterEvalQ did the trick. For my 2nd question, I have updated my question with some reproducible code. If that code is run multiple times by changing thresholds for strength (in last 2 lines of code) or changing parameters like restart, perturb or R, you will find that in all such experiments version 2 (i.e. b2) will have more no. of arcs than version 1 (i.e. b1) systematically. Any suggestion would be very helpful – Kumar Sep 23 '22 at 15:37
  • @user20650, For my 1st question, if I use clusterEvalQ as suggested, will all the bootstraps be randomly chosen. For eg if I have only 4 cores but I am running 50 bootstraps, will all the 50 bootstraps be randomly chosen ? Or will only the 1st 4 be randomly chosen and the following bootstraps will be a repeat of those 1st 4 bootstraps due to setting the seed? If the 2nd question is true, how do I randomly choose all 50 bootstraps but still make it reproducible? – Kumar Jun 23 '23 at 13:20
  • Kumar; you will want the sampling to be reproducible, so you will want to set a seed. I'd think the example i show above may (probably will) have the same samples from each core. which you don't want. When using parallel, there are specific function to set the seed ... i can't remember of hand but quick search should find them. Have a try and see if you can find it but feel free to give me a ping if needed. – user20650 Jun 23 '23 at 13:56
  • actually see https://www.bnlearn.com/examples/pkg-parallel/ – user20650 Jun 23 '23 at 14:00
  • .. which sasys "clusterSetRNGStream() to get reproducible results across the whole cluster...." .So i think `cl = makeCluster(2) ; clusterSetRNGStream(cl, 73810982) ; b2 = boot.strength(data = learning.test, R = 500, algorithm = "hc", cluster = cl ) ; stopCluster(cl)` should be good – user20650 Jun 23 '23 at 14:08
  • Thanks a lot @user20650. Will this single seed value of 73810982 in your eg take care of generating 500 different reproducible seeds for 500 bootstraps in your example? – Kumar Jun 26 '23 at 05:52
  • no. and you don't want to do that. For example; if you were wanting a reproducible uniform sample, you would set the seed once and then sample. e.g. set.seed(1); runif(5).. You would not use a different seed for each of the five samples ... and I believe it can mess with the pseudo-randomness (on rare occasions). Setting it once should be enough for reproducible results. – user20650 Jun 26 '23 at 13:22
  • @user20650, In your previous example of 500 bootstraps, we would want to sample 500 different bootstraps from the same dataset. Will setting a single seed not result in drawing the same sample for 500 bootstraps? i.e. all 500 boostraps would be the same dataset. – Kumar Jun 27 '23 at 06:22
  • no. From above "clusterSetRNGStream() to get reproducible results across the whole cluster...." . And you can test: if it was the same sample on each comute node then the samples would be the same `cl = makeCluster(2); clusterSetRNGStream(cl, 42); clusterEvalQ(cl, runif(10)); stopCluster(cl)`. Same aplies to the sampling of the rows of the data. – user20650 Jun 27 '23 at 09:17
  • Thanks a lot @user20650 for all your guidance. Its generating different samples each time and is reproducible – Kumar Jun 28 '23 at 04:36

0 Answers0