Analyzing the output
run_bbni() returns a list containing the sampled path of
the Markov chain:
log_posterior- numeric vector of log-posterior values stored after each node-level update.networks- list of sampled transition-function matrices. Each element shares the same format astrue_network.
Convergence: trace plot
A trace plot of the log-posterior over MCMC updates evaluates chain
behavior. Given that we know the specific ground-truth network, we can
further compute the log-posterior of the true network via the package’s
internal Error_LLH() function and graph it as a reference
baseline:
# Evaluate log-posterior value of the true network
# Call the unexported package function to compute a reference log-posterior
# for the ground-truth network.
true_logpost_raw <- BBNI:::Error_LLH(
TRFUM = true_network, GeneData = simulated_data,
SampleSize = sample_size, num.node = num_nodes,
prior_para = prior_para, penalty = penalty
)
# The log-posterior value is the last element of the first returned vector.
true_logpost <- true_logpost_raw[[1]][length(true_logpost_raw[[1]])]
# Plot Markov chain trace against individual node-level updates.
plot(mcmc_results$log_posterior,
type = "l",
xlab = "MCMC Step (one per node update)", ylab = "Log-Posterior",
main = "MCMC Trace Plot",
ylim = range(c(mcmc_results$log_posterior, true_logpost))
)
# Add true-network reference value.
abline(h = true_logpost, col = "firebrick", lty = 2, lwd = 2)
legend("bottomright",
legend = "True network log-posterior",
col = "firebrick", lty = 2, lwd = 2, bty = "n"
)
Note that the horizontal axis has units of node-level updates rather than outer MCMC iterations. The convergence of the log-posterior towards the baseline value is an indication of chain success. With statistical variation being inevitable, a perfect match with the true topology and Boolean functions is not expected from a single posterior sample.
Final sampled network
The final sampled network can be outputted as follows:
final_network <- tail(mcmc_results$networks, 1)[[1]]
final_network
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]
#> [1,] 0 0 0 0 0 0 0 0 0 11 0 0 0 0 0 0 0 0 0 0
#> [2,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#> [3,] 0 0 0 0 0 0 0 11 0 0 0 0 0 0 0 0 0 0 0 0
#> [4,] 0 0 0 0 0 0 12 0 0 0 0 0 0 0 0 0 0 0 0 0
#> [5,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6 0 0 0 0 6
#> [6,] 0 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#> [7,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 11 0 0 0 0
#> [8,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#> [9,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12 0 0 0
#> [10,] 0 0 0 0 0 8 0 0 0 0 0 0 0 0 0 0 0 0 0 8
#> [11,] 0 0 0 0 0 0 3 3 0 0 0 0 0 0 0 0 0 0 0 0
#> [12,] 0 0 0 0 0 4 0 0 0 0 0 0 4 0 0 0 0 0 0 0
#> [13,] 5 0 0 0 0 0 0 0 0 0 0 0 0 0 5 0 0 0 0 0
#> [14,] 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#> [15,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#> [16,] 0 0 0 0 0 0 0 0 0 0 0 2 2 0 0 0 0 0 0 0
#> [17,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
#> [18,] 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 2 0 0 0
#> [19,] 0 0 0 0 0 0 0 0 0 0 0 0 0 12 0 0 0 0 0 0
#> [20,] 0 0 7 0 0 0 0 7 0 0 0 0 0 0 0 0 0 0 0 0Recovery of the true network
As we know the true network, we can quantify the recovery by comparing the final sampled network with the true edge structure and Boolean transition functions:
true_edges <- true_network > 0
final_edges <- final_network > 0
true_positive <- sum(true_edges & final_edges) # correctly recovered edges
false_positive <- sum(!true_edges & final_edges) # spurious edges
false_negative <- sum(true_edges & !final_edges) # missed true edges
cat("True edges in the network: ", sum(true_edges), "\n")
#> True edges in the network: 29
cat("Correctly recovered edges: ", true_positive, "\n")
#> Correctly recovered edges: 23
cat("Spurious (false positive) edges: ", false_positive, "\n")
#> Spurious (false positive) edges: 1
cat("Missed (false negative) edges: ", false_negative, "\n")
#> Missed (false negative) edges: 6
# Among the correctly-recovered edges, how many also got the right
# Boolean function?
correct_function_final <- sum(true_edges & final_edges & (true_network == final_network))
cat("...of which function also correct: ", correct_function_final, "out of", true_positive, "\n")
#> ...of which function also correct: 20 out of 23Once again, some deviation is reasonable because this comparison uses only one sample from the end of the Markov chain. A posterior summary across multiple samples may allow a more stable estimate.
Posterior summarization by Bayesian model averaging
Instead of relying on a single sampled network, we can summarize the posterior distribution by Bayesian model averaging (BMA) over several retained network samples. This allows calculation of posterior edge-inclusion probabilities, which are often more informative than a single end sample. In biological applications, reducing false-positive edges is particularly important, as each predicted interaction may warrant costly experimental validation.
We first discard an initial burn-in period, here set to the first 75%
of the recorded chain. This is to reduce the impact of unpredictability
of the arbitrary starting structure. We then retain one network per
outer iteration (every num.node node-level updates) to
account for the strong statistical dependence of within-iteration
samples. For each directed edge, the posterior inclusion probability is
calculated as the fraction of retained samples in which that edge
appears:
# Set burn-in to 75%, determined by visual inspection of the trace plot above.
burnin_cutoff <- round(0.75 * length(mcmc_results$log_posterior))
# Retain one sampled network per full outer iteration.
keep_idx <- seq(burnin_cutoff + 1, length(mcmc_results$log_posterior), by = num_nodes)
cat("Post-burn-in, pruned samples retained:", length(keep_idx), "\n")
#> Post-burn-in, pruned samples retained: 1000
# Estimate posterior edge-inclusion probabilities across retained samples.
edge_prob <- Reduce(`+`, lapply(mcmc_results$networks[keep_idx], function(m) (m > 0) * 1)) / length(keep_idx)
bma_edges <- edge_prob > 0.5 # classify an edge as present if its posterior probability exceeds 0.5
bma_tp <- sum(true_edges & bma_edges)
bma_fp <- sum(!true_edges & bma_edges)
bma_fn <- sum(true_edges & !bma_edges)
cat("Correctly recovered edges: ", bma_tp, "\n")
#> Correctly recovered edges: 22
cat("Spurious (false positive) edges: ", bma_fp, "\n")
#> Spurious (false positive) edges: 0
cat("Missed (false negative) edges: ", bma_fn, "\n")
#> Missed (false negative) edges: 7
# To evaluate Boolean function accuracy, take the posterior mode (most frequently
# sampled function code) for each matrix coordinate across the retained samples,
# then construct a consensus network topology.
get_mode_nonzero <- function(v) {
nz <- v[v > 0]
if (length(nz) == 0) return(0)
tab <- table(nz)
as.numeric(names(tab)[which.max(tab)])
}
sample_array <- simplify2array(mcmc_results$networks[keep_idx])
bma_function <- apply(sample_array, c(1, 2), get_mode_nonzero)
bma_network <- matrix(0, nrow = num_nodes, ncol = num_nodes)
bma_network[bma_edges] <- bma_function[bma_edges]
correct_function_bma <- sum(true_edges & bma_edges & (true_network == bma_network))
cat("...of which Boolean function also correct: ", correct_function_bma, "out of", bma_tp, "\n")
#> ...of which Boolean function also correct: 19 out of 22
missed_edge_idx <- which(true_edges & !final_edges)
fp_edge_idx <- which(!true_edges & final_edges)
cat("Posterior probability of the missed true edge:", edge_prob[missed_edge_idx], "\n")
#> Posterior probability of the missed true edge: 0.001 0.024 0.211 0.011 0.654 0.016
cat("Posterior probability of the false-positive edge:", edge_prob[fp_edge_idx], "\n")
#> Posterior probability of the false-positive edge: 0.358Visualizing network recovery
The following plot compares the true network, the final sampled network, and the BMA consensus network. Edges that are structurally correct but assigned the wrong Boolean function are marked differently from edges for which both the parent relation and function are correctly recovered:
# Evaluate each matrix entry compared to the true network: correct absence,
# fully correct edge, correct edge with wrong function, false positive,
# or missed edge (false negative).
classify_match <- function(true_net, sample_net) {
true_e <- true_net > 0
samp_e <- sample_net > 0
out <- matrix(0, nrow = nrow(true_net), ncol = ncol(true_net))
out[!true_e & !samp_e] <- 0 # correct absence
out[true_e & samp_e & (true_net == sample_net)] <- 1 # edge + function correct
out[true_e & samp_e & (true_net != sample_net)] <- 2 # edge correct, function wrong
out[!true_e & samp_e] <- 3 # false positive
out[true_e & !samp_e] <- 4 # missed (false negative)
out
}
final_match <- classify_match(true_network, final_network)
bma_match <- classify_match(true_network, bma_network)
match_cols <- c("white", "darkgreen", "orange", "firebrick", "gray70")
match_breaks <- c(-0.5, 0.5, 1.5, 2.5, 3.5, 4.5)
match_labels <- c("Edge + function correct", "Edge correct, function wrong",
"False positive", "Missed edge")
layout(matrix(c(1, 2, 3, 4, 4, 4), nrow = 2, byrow = TRUE), heights = c(4, 1))
par(mar = c(2, 2, 3, 1))
image(true_edges * 1, col = c("white", "steelblue4"), breaks = c(-0.5, 0.5, 1.5),
axes = FALSE, main = "True Network")
image(final_match, col = match_cols, breaks = match_breaks, axes = FALSE, main = "Final Sample")
image(bma_match, col = match_cols, breaks = match_breaks, axes = FALSE, main = "Model-Averaged")
par(mar = c(0, 0, 0, 0))
plot.new()
legend("center", legend = match_labels, fill = match_cols[2:5], horiz = TRUE, bty = "n")