## ----setup, message = FALSE, warning = FALSE---------------------------------- library(sbm) library(ggplot2) library(knitr) theme_set(theme_bw()) ## ----ari---------------------------------------------------------------------- ARI <- aricode::ARI ## ----text-helpers, include = FALSE-------------------------------------------- # Helpers for the comments of the results: the text is generated from the # results, so that it stays consistent with a new simulation. fmt <- function(x, d = 2) format(round(x, d), nsmall = d) quality <- function(a) { if (a > 0.95) "nearly perfectly" else if (a > 0.8) "well" else if (a > 0.5) "partially" else "poorly" } compare <- function(gain, cost) { if (gain > 2 * cost) "much more than" else if (gain > cost) "more than" else "less than" } ## ----precomputed-fits, include = FALSE---------------------------------------- # The fits are long: they are precomputed by dev/precompute_SBM_nodes_covariates.R # and loaded here. If the file is missing, all the fits are run in the vignette. fitsFile <- "SBM_nodes_covariates_fits.rda" runFits <- !file.exists(fitsFile) if (!runFits) load(fitsFile) ## ----simu-param--------------------------------------------------------------- set.seed(3) n <- 120 # nodes Q <- 4 # blocks p <- 3 # covariates (including the intercept) X <- cbind(1, matrix(rnorm((p - 1) * n), ncol = p - 1)) colnames(X) <- c("intercept", "x1", "x2") head(X) ## ----simu-B------------------------------------------------------------------- B <- rbind(c(-0.01, 3, 4, 2), c( 5, 0.4, -2, 0), c(-0.7, 6, 0, -3)) B <- B - B[, Q] # last block as reference B <- B[, order(colSums(abs(B)), decreasing = TRUE)] rownames(B) <- colnames(X) B ## ----simu-sbm----------------------------------------------------------------- alpha <- matrix(rbeta(Q * Q, 1, 8), Q, Q) round(alpha, 2) mySampler <- sampleSimpleSBM( nbNodes = n, connectParam = list(mean = alpha), model = "bernoulli", directed = TRUE, nodesCovariates = X, nodesCovariatesParam = B ) Z <- mySampler$memberships Y <- mySampler$networkData table(Z) ## ----simu-plot---------------------------------------------------------------- plot(mySampler) ## ----simu-noisy--------------------------------------------------------------- epsilon <- 0.1 Noise <- matrix(rbinom(n * n, 1, epsilon), n, n) Ynoisy <- Y * (1 - Noise) + (1 - Y) * Noise ## ----fit-covar---------------------------------------------------------------- #| eval: !expr exists("runFits") && runFits # estimOptions <- list(verbosity = 0, plot = FALSE) # # fitCovar <- estimateSimpleSBM( # Y, model = "bernoulli", directed = TRUE, # nodesCovariates = X, # estimOptions = estimOptions # ) ## ----print-fit-covar---------------------------------------------------------- fitCovar ## ----fit-no-covar------------------------------------------------------------- #| eval: !expr exists("runFits") && runFits # fitNoCovar <- estimateSimpleSBM( # Y, model = "bernoulli", directed = TRUE, # estimOptions = estimOptions # ) ## ----nb-blocks---------------------------------------------------------------- c(noCovar = fitNoCovar$nbBlocks, covar = fitCovar$nbBlocks) ## ----table-clustering--------------------------------------------------------- table(estimated = fitNoCovar$memberships, true = Z) table(estimated = fitCovar$memberships, true = Z) ## ----text-clean, include = FALSE---------------------------------------------- ariNoCovar <- ARI(fitNoCovar$memberships, Z) ariCovar <- ARI(fitCovar$memberships, Z) sameClust <- abs(ariCovar - ariNoCovar) < 0.05 ## ----B-hat-------------------------------------------------------------------- round(fitCovar$nodesCovariatesParam, 2) ## ----align-B------------------------------------------------------------------ matchBlocks <- function(estimated, true) { apply(table(estimated, true), 2, which.max) } perm <- matchBlocks(fitCovar$memberships, Z) Bhat <- fitCovar$nodesCovariatesParam[, perm] Bhat <- Bhat - Bhat[, Q] rownames(Bhat) <- colnames(X) round(Bhat, 2) B ## ----plot-pi------------------------------------------------------------------ softmax <- function(eta) { eta <- eta - apply(eta, 1, max) exp(eta) / rowSums(exp(eta)) } probs <- data.frame( true = c(softmax(X %*% B)), estimated = c(softmax(X %*% Bhat)), block = factor(rep(1:Q, each = n)) ) ggplot(probs) + aes(x = true, y = estimated, color = block) + geom_point(alpha = 0.6) + geom_abline(intercept = 0, slope = 1) + labs(x = expression(pi[q](x[i])), y = expression(hat(pi)[q](x[i]))) ## ----fit-useless-------------------------------------------------------------- #| eval: !expr exists("runFits") && runFits # set.seed(123) # Xuseless <- cbind(1, matrix(rnorm(2 * n), ncol = 2)) # fitUseless <- estimateSimpleSBM( # Y, model = "bernoulli", directed = TRUE, # nodesCovariates = Xuseless, # estimOptions = estimOptions # ) ## ----res-useless-------------------------------------------------------------- resClean <- data.frame( covariates = c("none", "x1 + x2", "two useless covariates"), nbBlocks = c(fitNoCovar$nbBlocks, fitCovar$nbBlocks, fitUseless$nbBlocks), ARI = c(ARI(fitNoCovar$memberships, Z), ARI(fitCovar$memberships, Z), ARI(fitUseless$memberships, Z)), ICL = c(fitNoCovar$ICL, fitCovar$ICL, fitUseless$ICL) ) kable(resClean, digits = 2) ## ----text-useless, include = FALSE-------------------------------------------- dIclCovar <- fitCovar$ICL - fitNoCovar$ICL dIclUseless <- fitUseless$ICL - fitNoCovar$ICL txtSelect <- paste0( if (ariCovar - ariNoCovar < 0.05) "The covariates barely change the clustering. " else "The covariates improve the clustering. ", if (dIclCovar > 0) "The ICL selects" else "However, the ICL does not select", " the model with $x_1$ and $x_2$ (difference of ICL with the model without covariates: ", fmt(dIclCovar, 1), "), ", if (dIclUseless > 0) "and the useless covariates are selected too" else "whereas the useless covariates are not selected", " (difference of ICL: ", fmt(dIclUseless, 1), ")." ) ## ----logpZ-------------------------------------------------------------------- ElogpZ <- c( none = sum(fitNoCovar$probMemberships%*%log(fitNoCovar$blockProp)), "x1 + x2" = sum(fitCovar$probMemberships*log(softmax(X %*% fitCovar$nodesCovariatesParam))) ) ElogpZ ## ----text-logpZ, include = FALSE---------------------------------------------- gainZ <- unname(diff(ElogpZ)) penZ <- 0.5 * (p - 1) * (fitCovar$nbBlocks - 1) * log(n) ## ----fit-noisy---------------------------------------------------------------- #| eval: !expr exists("runFits") && runFits # fitNoCovarNoisy <- estimateSimpleSBM( # Ynoisy, model = "bernoulli", directed = TRUE, # estimOptions = estimOptions # ) # fitCovarNoisy <- estimateSimpleSBM( # Ynoisy, model = "bernoulli", directed = TRUE, # nodesCovariates = X, # estimOptions = estimOptions # ) ## ----ARI-simple--------------------------------------------------------------- resSimple <- data.frame( network = rep(c("original", "noisy"), each = 2), model = rep(c("without covariates", "with covariates"), 2), nbBlocks = c(fitNoCovar$nbBlocks, fitCovar$nbBlocks, fitNoCovarNoisy$nbBlocks, fitCovarNoisy$nbBlocks), ARI = c(ARI(fitNoCovar$memberships, Z), ARI(fitCovar$memberships, Z), ARI(fitNoCovarNoisy$memberships, Z), ARI(fitCovarNoisy$memberships, Z)) ) kable(resSimple, digits = 2) ## ----text-noisy, include = FALSE---------------------------------------------- ariNoCovarNoisy <- ARI(fitNoCovarNoisy$memberships, Z) ariCovarNoisy <- ARI(fitCovarNoisy$memberships, Z) txtNoisy <- paste0( if (sameClust && min(ariNoCovar, ariCovar) > 0.9) "On the original network, the structure is strong enough to be recovered from the network alone, and both models give similar clusterings. " else sprintf("On the original network, the ARI is %s without covariates and %s with covariates. ", fmt(ariNoCovar), fmt(ariCovar)), sprintf("On the noisy network, the standard SBM selects %s%d blocks (ARI = %s)", if (fitNoCovarNoisy$nbBlocks < Q) "only " else "", fitNoCovarNoisy$nbBlocks, fmt(ariNoCovarNoisy)), if (ariNoCovarNoisy < 0.8) ": the network alone does not carry enough information. " else ". ", if (ariCovarNoisy > ariNoCovarNoisy + 0.05) sprintf("The covariates bring the missing information: the SBM with covariates selects %d blocks (ARI = %s).", fitCovarNoisy$nbBlocks, fmt(ariCovarNoisy)) else sprintf("The SBM with covariates selects %d blocks (ARI = %s): here the covariates do not improve the clustering.", fitCovarNoisy$nbBlocks, fmt(ariCovarNoisy)) ) ## ----test-covar--------------------------------------------------------------- #| eval: !expr exists("runFits") && runFits # x3 <- rnorm(n) # covarSets <- list( # "none" = NULL, # "x1" = X[, c("intercept", "x1")], # "x2" = X[, c("intercept", "x2")], # "x1 + x2" = X, # "x1 + x2 + x3" = cbind(X, x3 = x3) # ) # fits <- lapply(covarSets, function(covar) { # if (is.null(covar)) { # estimateSimpleSBM(Ynoisy, model = "bernoulli", directed = TRUE, # estimOptions = estimOptions) # } else { # estimateSimpleSBM(Ynoisy, model = "bernoulli", directed = TRUE, # nodesCovariates = covar, estimOptions = estimOptions) # } # }) ## ----res-test-covar----------------------------------------------------------- resICL <- data.frame( covariates = names(fits), nbBlocks = sapply(fits, function(f) f$nbBlocks), ICL = sapply(fits, function(f) f$ICL), ARI = sapply(fits, function(f) ARI(f$memberships, Z)) ) resICL$deltaICL <- resICL$ICL - resICL$ICL[1] kable(resICL, digits = 2, row.names = FALSE) ## ----text-test-covar, include = FALSE----------------------------------------- covarLabels <- c("none" = "no covariate", "x1" = "$x_1$ only", "x2" = "$x_2$ only", "x1 + x2" = "$x_1$ and $x_2$", "x1 + x2 + x3" = "$x_1$, $x_2$ and $x_3$") dICL <- setNames(resICL$deltaICL, resICL$covariates) best <- resICL$covariates[which.max(resICL$ICL)] alone <- dICL[c("x1", "x2")] > 0 txtTest <- paste0( if (all(alone)) "Here, each of $x_1$ and $x_2$ has an effect on its own" else if (any(alone)) sprintf("Here, only $x_%d$ has an effect on its own", which(alone)) else "Here, neither $x_1$ nor $x_2$ improves the ICL on its own", if (best == "x1 + x2") ", and the best model (highest ICL) is the one with both of them, which is the true model. " else sprintf(", and the best model (highest ICL) is the one with %s, whereas the true model involves $x_1$ and $x_2$. ", covarLabels[best]), if (dICL["x1 + x2 + x3"] < dICL["x1 + x2"]) "Adding the irrelevant covariate $x_3$ decreases the ICL: its small gain in fit does not compensate for the $Q - 1$ additional parameters." else "Adding the irrelevant covariate $x_3$ increases the ICL here: by chance, it slightly helps to predict the memberships of this sample." ) ## ----simu-bipartite----------------------------------------------------------- set.seed(4) n1 <- 120; Q1 <- 4 n2 <- 100; Q2 <- 3 X1 <- cbind(1, matrix(rnorm(2 * n1), ncol = 2)) # row covariates X2 <- cbind(1, rnorm(n2)) # column covariates B1 <- rbind(c(-0.01, 3, 4, 2), c( 5, 0.4, -2, 0), c(-0.7, 6, 0, -3)) B1 <- B1 - B1[, Q1] B1 <- B1[, order(colSums(abs(B1)), decreasing = TRUE)] B2 <- rbind(c(0, 1, 1), c(2, 0, 0.4)) B2 <- B2 - B2[, Q2] B2 <- B2[, order(colSums(abs(B2)), decreasing = TRUE)] alphaBip <- matrix(0, Q1 , Q2) alphaBip[,1] = c(0.2,0.4,0.5,0.05) alphaBip[,2] = c(0.5,0.2,0.05,0.05) alphaBip[,3] = c(0.1,0.2,0.05,0.3) myBipSampler <- sampleBipartiteSBM( nbNodes = c(n1, n2), connectParam = list(mean = alphaBip), model = "bernoulli", dimLabels = c("row", "col"), nodesCovariates = list(X1, X2), nodesCovariatesParam = list(B1, B2) ) Zbip <- myBipSampler$memberships Ybip <- myBipSampler$networkData lapply(Zbip, table) ## ----simu-bipartite-plot------------------------------------------------------ plot(myBipSampler) ## ----simu-bipartite-NA-------------------------------------------------------- epsilonNA <- c(0.1, 0.3, 0.6,0.8) YbipNA <- lapply(epsilonNA, function(eps) { Yna <- Ybip Yna[matrix(rbinom(n1 * n2, 1, eps), n1, n2) == 1] <- NA Yna }) names(YbipNA) <- paste0(100 * epsilonNA, "% NA") sapply(YbipNA, function(Yna) mean(is.na(Yna))) ## ----simu-plot-bipartite-noise------------------------------------------------ plotMyMatrix(YbipNA[[3]],clustering=Zbip) ## ----fit-bipartite------------------------------------------------------------ #| eval: !expr exists("runFits") && runFits # networksBip <- c(list(original = Ybip), YbipNA) # fitBip <- lapply(networksBip, function(Yobs) { # list( # noCovar = estimateBipartiteSBM( # Yobs, model = "bernoulli", dimLabels = c("row", "col"), # estimOptions = estimOptions), # covar = estimateBipartiteSBM( # Yobs, model = "bernoulli", dimLabels = c("row", "col"), # nodesCovariates = list(X1, X2), estimOptions = estimOptions) # ) # }) ## ----print-fit-bipartite------------------------------------------------------ fitBip$original$covar ## ----bipartite-B-hat---------------------------------------------------------- lapply(fitBip$original$covar$nodesCovariatesParam, round, 2) ## ----ARI-bipartite------------------------------------------------------------ resBip <- do.call(rbind, lapply(names(fitBip), function(net) { data.frame( network = net, model = c("without covariates", "with covariates"), nbRowBlocks = sapply(fitBip[[net]], function(f) f$nbBlocks[1]), nbColBlocks = sapply(fitBip[[net]], function(f) f$nbBlocks[2]), ARIrow = sapply(fitBip[[net]], function(f) ARI(f$memberships[[1]], Zbip[[1]])), ARIcol = sapply(fitBip[[net]], function(f) ARI(f$memberships[[2]], Zbip[[2]])), deltaICL = sapply(fitBip[[net]], function(f) f$ICL - fitBip[[net]]$noCovar$ICL) ) })) kable(resBip, digits = 2, row.names = FALSE) ## ----ARI-bipartite-summary, include = FALSE----------------------------------- lab <- names(YbipNA) ## comments of the results, generated from them nbB <- function(net, mod) fitBip[[net]][[mod]]$nbBlocks ariB <- function(net, mod, side) ARI(fitBip[[net]][[mod]]$memberships[[side]], Zbip[[side]]) blocksTxt <- function(nb) sprintf("%d $\\times$ %d block%s", nb[1], nb[2], if (prod(nb) > 1) "s" else "") pct <- function(i) paste0(100 * epsilonNA[i], "%") # original network ariOrig <- c(ariB("original", "noCovar", 1), ariB("original", "noCovar", 2), ariB("original", "covar", 1), ariB("original", "covar", 2)) trueNb <- function(net) all(nbB(net, "noCovar") == c(Q1, Q2)) && all(nbB(net, "covar") == c(Q1, Q2)) improved <- any(ariOrig[3:4] - ariOrig[1:2] > 0.05) dIclOrig <- fitBip$original$covar$ICL - fitBip$original$noCovar$ICL txtBipOrig <- paste0( if (trueNb("original") && min(ariOrig) > 0.95) "On the original network, the structure is strong enough to be recovered from the network alone: both models find the $Q_1 = 4$ row blocks and the $Q_2 = 3$ column blocks with almost perfect clusterings. " else sprintf("On the original network, the model without covariates finds %s (ARI = %s and %s for rows and columns), and the model with covariates %s (ARI = %s and %s). ", blocksTxt(nbB("original", "noCovar")), fmt(ariOrig[1]), fmt(ariOrig[2]), blocksTxt(nbB("original", "covar")), fmt(ariOrig[3]), fmt(ariOrig[4])), if (!improved && dIclOrig > 0) "The covariates do not improve the clustering, but the ICL still selects them, since they predict the memberships." else if (!improved) "The covariates do not improve the clustering, and the ICL does not select them." else if (dIclOrig > 0) "The covariates improve the clustering, and the ICL selects them." else "The covariates improve the clustering, but the ICL does not select them." ) # moderate levels of missing data ariNoCovRow <- sapply(lab, function(l) ariB(l, "noCovar", 1)) ariCovRow <- sapply(lab, function(l) ariB(l, "covar", 1)) txtBipModerate <- if (all(ariNoCovRow[1:2] > 0.9) && min(ariOrig[1:2]) > 0.9) sprintf("The SBM without covariates is robust to moderate levels of missing data: with %s and %s of missing entries, the clusterings are nearly unchanged (row ARI = %s and %s).", pct(1), pct(2), fmt(ariNoCovRow[1]), fmt(ariNoCovRow[2])) else sprintf("With %s and %s of missing entries, the SBM without covariates gives row ARI = %s and %s, and the SBM with covariates %s and %s.", pct(1), pct(2), fmt(ariNoCovRow[1]), fmt(ariNoCovRow[2]), fmt(ariCovRow[1]), fmt(ariCovRow[2])) # high level of missing data txtBip3 <- paste0( sprintf("With %s of missing entries, ", pct(3)), if (ariNoCovRow[3] < 0.9) "the network alone is no longer informative enough: " else "", sprintf("the standard SBM selects %d row blocks (ARI = %s)", nbB(lab[3], "noCovar")[1], fmt(ariNoCovRow[3])), if (ariCovRow[3] > ariNoCovRow[3] + 0.05) sprintf(", whereas the SBM with covariates finds %d row blocks (ARI = %s).", nbB(lab[3], "covar")[1], fmt(ariCovRow[3])) else sprintf(", and the SBM with covariates %d row blocks (ARI = %s).", nbB(lab[3], "covar")[1], fmt(ariCovRow[3])) ) # very high level of missing data nb4 <- nbB(lab[4], "noCovar") txtBip4 <- paste0( sprintf("With %s of missing entries, ", pct(4)), if (all(nb4 == 1)) "the network alone does not reveal any structure (1 $\\times$ 1 block)" else sprintf("the standard SBM finds %s (ARI = %s and %s)", blocksTxt(nb4), fmt(ariNoCovRow[4]), fmt(ariB(lab[4], "noCovar", 2))), sprintf(", while the SBM with covariates finds %s (ARI = %s and %s for rows and columns).", blocksTxt(nbB(lab[4], "covar")), fmt(ariCovRow[4]), fmt(ariB(lab[4], "covar", 2))) ) # gain brought by the covariates and ICL gainRow <- c(ariOrig[3] - ariOrig[1], ariCovRow - ariNoCovRow) dIclBip <- sapply(names(fitBip), function(net) fitBip[[net]]$covar$ICL - fitBip[[net]]$noCovar$ICL) txtBipEnd <- paste0( if (all(diff(gainRow) > -0.02) && gainRow[length(gainRow)] > gainRow[1] + 0.05) "The more entries are missing, the more the information brought by the covariates matters. " else "", if (all(dIclBip > 0)) "For all the networks, the ICL selects the model with covariates." else if (all(dIclBip <= 0)) "For none of the networks does the ICL select the model with covariates." else sprintf("The ICL selects the model with covariates only for the following networks: %s.", paste(names(fitBip)[dIclBip > 0], collapse = ", ")) ) ## ----fit-bipartite-row-only--------------------------------------------------- #| eval: !expr exists("runFits") && runFits # fitBipRowOnly <- estimateBipartiteSBM( # YbipNA[[lab[3]]], model = "bernoulli", dimLabels = c("row", "col"), # nodesCovariates = list(X1, matrix(0, 0, 0)), # estimOptions = estimOptions # ) ## ----nb-blocks-row-only------------------------------------------------------- fitBipRowOnly$nbBlocks ## ----ICL-bipartite------------------------------------------------------------ iclBip <- c(none = fitBip[[lab[3]]]$noCovar$ICL, rowsOnly = fitBipRowOnly$ICL, rowsAndCols = fitBip[[lab[3]]]$covar$ICL) iclBip