Search code examples
rbayesianbayesian-networksbnlearn

BNLearn: How to merge the estimating parameters of a Gaussian Bayesian network with its conditional structure?


I defined the structure of a Gaussian Bayesian network using the iamb function and then estimated the coefficients of the nodes using bn.fit.


Library

library(bnlearn)

Data

{  C       E       G       N       V       W
48.83   51.48   42.64   54.1    42.96   41.96
48.85   73.43   40.97   60.07   65.29   48.96
67.01   71.1    52.52   51.64   63.22   62.03
37.83   49.33   56.15   49.01   47.75   38.77
55.3    49.27   63.55   54.62   60.57   56.66
56.12   48.72   66.02   43.95   55.54   52.39}

Code

# Definition of mandatory and forbidden nodes - here the white list
wl = data.frame(from = c("E","G","V","W","N"), to = c("V", "V","W","C","C"))

# Definition of the constrained network
network <- iamb(Data, test = "cor", whitelist = wl)

# Estimation of the coefficients according to the structure of the network
est.para <- bn.fit(network, data = Data)

The problem is that est.para is a list and not a GBN that can be plotted, etc.. I would like to know how to merge the network and the estimated parameters?


Solution

  • If you want to have some network plot that displays some extra information besides the connections, you could use strength.plot. Following your example:

    library(Rgraphviz)
    
    strength <- arc.strength(network, Data)
    strength.plot(network, strength, shape = "ellipse")
    

    In case it is absolutely necessary to use the results of the parameters of the GBN est.para, you could use the graphviz.plot arguments to highlight edges and nodes (can be done using edgeRenderInfo and nodeRenderInfo). Just as an example, you could use the parameters to choose the width of the edges:

    library(data.table) 
    
    plot <- graphviz.plot(network, shape = "ellipse")
    
    arc.sizes <- data.table(network$arcs)
    arc.sizes[, edge.name := paste0(arc.sizes$from, "~", arc.sizes$to)]
    arc.sizes[, param := abs(est.para[[to]]$coefficients[[from]]), by = .(from, to)]
    arc.sizes[, lwd := 5*((param - min(param))/(max(param) - min(param)))]
    
    lwd <- as.vector(arc.sizes$lwd)
    names(lwd) <- arc.sizes$edge.name
    edgeRenderInfo(plot) <- list(lwd = lwd)
    
    renderGraph(plot)