# read in data and true network
data(example)
Y.n=example$Y # data matrix
true.dir=example$true.dir #adjacency matrix of the data generating DAG
true.ske=example$true.ske # skeleton graph of the data generating DAG
## (i) DAG learning by hill climbing: no aggregation
# learn DAG using "BIC" score
temp=score(Y=Y.n, n.boot=0, score.type="BIC")
adj=temp$adj.matrix
# Find the total number of skeleton edges
# and the number of correct skeleton edges of the estimated DAG
tt=adj+t(adj) ##symmetrization
correct.c=sum((tt>0)&(true.ske>0))/2
total.c=sum(tt>0)/2
total.true=sum(true.ske>0)/2
c(total.c, correct.c, total.true) ##326, 96, 109
## (ii) DAG learning by bagging
set.seed(1)
# get the ensemble of DAGs from bootstrap resamples
temp.boot=score(Y.n, n.boot=10, score.type="BIC")
boot.adj=temp.boot$adj.matrix
temp.bag=score_shd(boot.adj, alpha = 1) ##aggregating
adj.bag=temp.bag$adj.matrix
# Find the total number of skeleton edges and the number of
# correct skeleton edges of the estimated DAG
tt=adj.bag+t(adj.bag) ##symmetrization
correct.c=sum((tt>0)&(true.ske>0))/2
total.c=sum(tt>0)/2
total.true=sum(true.ske>0)/2
c(total.c, correct.c, total.true) ##121, 89, 109
#Note: much less number of false edges and slightly less number of correct edges
#compared with the non-bagging result
## try different aggregation parameter alpha
#(i) alpha=0.5: less aggressive, more edges retained
temp.bag=score_shd(boot.adj, alpha = 0.5) ##aggregating
adj.bag05=temp.bag$adj.matrix
# Find the total number of skeleton edges
# and the number of correct skeleton edges of the estimated DAG
tt=adj.bag05+t(adj.bag05) ##symmetrization
correct.c=sum((tt>0)&(true.ske>0))/2
total.c=sum(tt>0)/2
total.true=sum(true.ske>0)/2
c(total.c, correct.c, total.true) ##132, 91, 109
#(ii) alpha=1.5: more aggressive, less edges retained
temp.bag=score_shd(boot.adj, alpha = 1.5) ##aggregating
adj.bag15=temp.bag$adj.matrix
# Find the total number of skeleton edges
# and the number of correct skeleton edges of the estimated DAG
tt=adj.bag15+t(adj.bag15) ##symmetrization
correct.c=sum((tt>0)&(true.ske>0))/2
total.c=sum(tt>0)/2
total.true=sum(true.ske>0)/2
c(total.c, correct.c, total.true) ##115, 87, 109
#(iii) alpha=2: more aggressive, less edges retained
temp.bag=score_shd(boot.adj, alpha = 2) ##aggregating
adj.bag2=temp.bag$adj.matrix
# Find the total number of skeleton edges
# and the number of correct skeleton edges of the estimated DAG
tt=adj.bag2+t(adj.bag2) ##symmetrization
correct.c=sum((tt>0)&(true.ske>0))/2
total.c=sum(tt>0)/2
total.true=sum(true.ske>0)/2
c(total.c, correct.c, total.true) ##86, 67, 109
Run the code above in your browser using DataLab