R Under development (unstable) (2026-07-17 r90265 ucrt) -- "Unsuffered Consequences" Copyright (C) 2026 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > # Backward-compatibility reference values captured from bootnet 1.8 > # (commit 580f1fe) on R 4.6.1 — these results must not change unless a > # NEWS-documented, intentional break is made. > # > # All data are generated with fixed seeds; every reference value below was > # produced twice from the installed bootnet 1.8 and verified to be bitwise > # identical across runs (including mgm, which is seeded immediately before > # the call). > > library(bootnet) Loading required package: ggplot2 This is bootnet 1.9.1 For questions and issues, please see github.com/SachaEpskamp/bootnet. > > ## --------------------------------------------------------------------------- > ## Data generation (fixed seeds — do not change) > ## --------------------------------------------------------------------------- > > # Gaussian data: 6 variables, 200 rows, Toeplitz correlation matrix 0.6^|i-j| > genGauss <- function(){ + Sigma <- 0.6 ^ abs(outer(1:6, 1:6, "-")) + set.seed(20260716) + d <- as.data.frame(mvtnorm::rmvnorm(200, sigma = Sigma)) + colnames(d) <- paste0("V", 1:6) + d + } > > # Binary data: 5 nodes, chain Ising graph (weight 0.7), thresholds -0.7, n = 500 > genBinary <- function(){ + ig <- matrix(0, 5, 5) + ig[cbind(1:4, 2:5)] <- 0.7 + ig <- ig + t(ig) + set.seed(31415) + d <- as.data.frame(IsingSampler::IsingSampler(500, ig, rep(-0.7, 5))) + colnames(d) <- paste0("I", 1:5) + d + } > > gaussData <- genGauss() > binData <- genBinary() > > ## --------------------------------------------------------------------------- > ## a. EBICglasso > ## --------------------------------------------------------------------------- > > net_ebic <- estimateNetwork(gaussData, default = "EBICglasso", verbose = FALSE) > > ref_ebic_graph <- structure(c(0, 0.5013921896, 0, 0.0458909248, 0, 0, 0.5013921896, + 0, 0.3989122116, 0.0509131063, 0.03146761, 0, 0, 0.3989122116, + 0, 0.3440014768, 0, 0.0105654823, 0.0458909248, 0.0509131063, + 0.3440014768, 0, 0.3829775915, 0.0133346708, 0, 0.03146761, 0, + 0.3829775915, 0, 0.4493225882, 0, 0, 0.0105654823, 0.0133346708, + 0.4493225882, 0), dim = c(6L, 6L)) > > stopifnot(isTRUE(all.equal(round(unname(net_ebic$graph), 10), + ref_ebic_graph, tolerance = 1e-8))) > cat("OK: EBICglasso graph\n") OK: EBICglasso graph > > ## --------------------------------------------------------------------------- > ## b. ggmModSelect > ## --------------------------------------------------------------------------- > > net_gms <- estimateNetwork(gaussData, default = "ggmModSelect", verbose = FALSE) > > ref_gms_graph <- structure(c(0, 0.5588904338, 0, 0, 0, 0, 0.5588904338, 0, 0.4519451477, + 0, 0, 0, 0, 0.4519451477, 0, 0.4105052623, 0, 0, 0, 0, 0.4105052623, + 0, 0.4294993889, 0, 0, 0, 0, 0.4294993889, 0, 0.5029608362, 0, + 0, 0, 0, 0.5029608362, 0), dim = c(6L, 6L)) > > stopifnot(isTRUE(all.equal(round(unname(net_gms$graph), 10), + ref_gms_graph, tolerance = 1e-8))) > cat("OK: ggmModSelect graph\n") OK: ggmModSelect graph > > ## --------------------------------------------------------------------------- > ## c. IsingFit > ## --------------------------------------------------------------------------- > > net_ising <- suppressWarnings( + estimateNetwork(binData, default = "IsingFit", verbose = FALSE)) > > ref_ising_graph <- structure(c(0, 0.8472137518, 0, 0, 0, 0.8472137518, 0, 0.3789165163, + 0, 0, 0, 0.3789165163, 0, 0.6413969983, 0, 0, 0, 0.6413969983, + 0, 0.4259031921, 0, 0, 0, 0.4259031921, 0), dim = c(5L, 5L)) > > ref_ising_intercepts <- c(-0.7827314985, -0.7543547939, -0.4739905321, + -0.7296237443, -0.5371955831) > > stopifnot(isTRUE(all.equal(round(unname(net_ising$graph), 10), + ref_ising_graph, tolerance = 1e-8))) > stopifnot(isTRUE(all.equal(round(unname(unlist(net_ising$intercepts)), 10), + ref_ising_intercepts, tolerance = 1e-8))) > cat("OK: IsingFit graph + intercepts\n") OK: IsingFit graph + intercepts > > ## --------------------------------------------------------------------------- > ## d. mgm (seed set immediately before the call; criterion = EBIC default) > ## --------------------------------------------------------------------------- > > set.seed(42) > net_mgm <- suppressWarnings(suppressMessages( + estimateNetwork(binData, default = "mgm", verbose = FALSE))) > > ref_mgm_graph <- structure(c(0, 0.423607139, 0, 0, 0, 0.423607139, 0, 0.1894588907, + 0, 0, 0, 0.1894588907, 0, 0.3206992763, 0, 0, 0, 0.3206992763, + 0, 0.2129518263, 0, 0, 0, 0.2129518263, 0), dim = c(5L, 5L)) > > stopifnot(isTRUE(all.equal(round(unname(net_mgm$graph), 10), + ref_mgm_graph, tolerance = 1e-8))) > cat("OK: mgm graph\n") OK: mgm graph > > ## --------------------------------------------------------------------------- > ## e. Bootstrap results (nonparametric + person/case-drop) on EBICglasso net > ## --------------------------------------------------------------------------- > > set.seed(1234) > boot_np <- suppressMessages( + bootnet(net_ebic, nBoots = 6, type = "nonparametric", + statistics = c("edge", "strength"), verbose = FALSE)) > > # Sorted bootTable values (sorting makes this robust to row ordering) > ref_np_values <- c(-0.13939805, -0.10966104, -0.06282389, -0.03739441, -0.03564222, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.00301201, 0.0088298, 0.00883062, + 0.01299946, 0.01586979, 0.01590527, 0.01954646, 0.03346298, 0.03604092, + 0.03668636, 0.03819781, 0.03957764, 0.03986875, 0.04030457, 0.04731852, + 0.05130405, 0.0623849, 0.06917367, 0.09438307, 0.10960173, 0.1112265, + 0.11175038, 0.11859931, 0.11885293, 0.1198825, 0.24607124, 0.28089999, + 0.30633193, 0.33044699, 0.34043903, 0.34078203, 0.35364731, 0.37190248, + 0.37300506, 0.37776775, 0.38969664, 0.40522278, 0.4107018, 0.41833464, + 0.41899725, 0.42136412, 0.42919942, 0.43874591, 0.46995338, 0.48533852, + 0.48840078, 0.4921327, 0.49328085, 0.49478216, 0.49496943, 0.51431166, + 0.51918522, 0.52293325, 0.52983258, 0.54763095, 0.54901693, 0.56858632, + 0.58465914, 0.5883564, 0.59526182, 0.59548093, 0.59597349, 0.6803367, + 0.69123116, 0.69573331, 0.69853425, 0.72050956, 0.75085337, 0.75088232, + 0.76085628, 0.76180316, 0.77848275, 0.79152637, 0.86249941, 0.8680938, + 0.86872975, 0.87372246, 0.87974868, 0.90398624, 0.90438782, 0.92202198, + 0.93543972, 0.93566801, 0.94559053, 0.97514951, 0.9868398, 1.0023142, + 1.01596065, 1.06289669, 1.08311673, 1.14571073) > > ref_np_dim <- c(126L, 12L) > ref_np_cols <- c("name", "type", "node1", "node2", "value", "id", "nNode", + "nPerson", "rank_avg", "rank_min", "rank_max", "graph") > > stopifnot(isTRUE(all.equal(round(sort(boot_np$bootTable$value), 8), + ref_np_values, tolerance = 1e-8))) > stopifnot(identical(dim(boot_np$bootTable), ref_np_dim)) > stopifnot(identical(colnames(boot_np$bootTable), ref_np_cols)) > cat("OK: nonparametric bootstrap values + bootTable structure\n") OK: nonparametric bootstrap values + bootTable structure > > set.seed(5678) > boot_cd <- suppressMessages( + bootnet(net_ebic, nBoots = 6, type = "person", caseN = 3, + verbose = FALSE)) > > ref_cs <- c(edge = 0.75, strength = 0.75) > > cs <- round(suppressWarnings(corStability(boot_cd, verbose = FALSE)), 8) > stopifnot(isTRUE(all.equal(cs[names(ref_cs)], ref_cs, tolerance = 1e-8))) > cat("OK: person-drop bootstrap CS-coefficients\n") OK: person-drop bootstrap CS-coefficients > > cat("All backward-compatibility reference checks passed.\n") All backward-compatibility reference checks passed. > > proc.time() user system elapsed 8.04 1.03 9.07